Calibration of Local-Stochastic Volatility Models by Optimal Transport
aa r X i v : . [ q -f i n . M F ] S e p Calibration of Local-Stochastic Volatility Models byOptimal Transport
Ivan Guo , Grégoire Loeper , and Shiyi Wang School of Mathematics, Monash University, Australia Centre for Quantitative Finance and Investment Strategies ∗ ,Monash University, AustraliaSeptember 25, 2019 Abstract
In this paper, we study a semi-martingale optimal transport problem and its ap-plication to the calibration of Local-Stochastic Volatility (LSV) models. Rather thanconsidering the classical constraints on marginal distributions at initial and final time,we optimise our cost function given the prices of a finite number of European options.We formulate the problem as a convex optimisation problem, for which we provide adual formulation. Then we solve numerically the dual problem, which involves a fullynon-linear Hamilton–Jacobi–Bellman equation. The method is tested by calibrating aHeston-like LSV model with simulated data and foreign exchange market data.
Since the introduction of the Black–Scholes model, a lot of effort has been put on developingsophisticated volatility models that properly capture the market dynamics. In the space ofequities and currencies, the most widely used models are the Local Volatility (LV) [10] and theStochastic Volatility (SV) models [13, 20]. Introduced as an extension of the Black–Scholesmodel, the LV model can be exactly calibrated to any arbitrage-free implied volatility surface.Despite this feature, the LV model has often been criticised for its unrealistic volatilitydynamics. The SV models tend to be more consistent with the market dynamics, but theystruggle to fit short term market smiles and skews, and being parametric, they do not haveenough degrees of freedom to match all vanilla market prices. A better fit can be obtainedby increasing the number of stochastic factors in the SV models; however, this also increasesthe complexity of calibration and pricing.Local-Stochastic Volatility (LSV) models, introduced in [23], naturally extend and takeadvantage of both approaches. The idea behind LSV models is to incorporate a local, non-parametric, factor into the SV models. Thus, while keeping consistent dynamics, the modelcan match all observed market prices (as long as one restricts to European claims). Thedetermination of this local factor (also called leverage ) is based on Gyöngy’s mimicking the-orem [17]. Research into the numerical calibration of LSV models has been developed in twodifferent directions. One is based on a Monte Carlo approach, with Henry-Labordere [18],followed by Henry-Labordere and Guyon [16] using a so-called McKean’s particle method. ∗ Acknowledgements
The Centre for Quantitative Finance and Investment Strategies has been sup-ported by BNP Paribas.
Preliminaries
Given a Polish space E equipped with its Borel σ -algebra, let C ( E ) be the space of con-tinuous functions on E and C b ( E ) be the space of bounded continuous functions. Denoteby M ( E ) the space of finite signed Borel measures endowed with the weak- ∗ topology. Let M + ( E ) ⊂ M ( E ) denote the subset of nonnegative measures. If E is compact, the topologi-cal dual of C b ( E ) is given by C b ( E ) ∗ = M ( E ) . More generally, if E is non-compact, C b ( E ) ∗ is larger than M ( E ) . Let P ( E ) be the space of Borel probability measures, BV ( E ) be thespace of functions of bounded variation and L ( dµ ) be the space of µ -integrable functions.We also write C b ( E, R d ) , M ( E, R d ) , BV ( E, R d ) and L ( dµ, R d ) as the vector-valued ver-sions of their corresponding spaces. For convenience, let Λ = [0 , T ] × R d and X = R × R d × S d .We use the notation h· , ·i to denote the duality bracket between C b (Λ , X ) and C b (Λ , X ) ∗ .Let Ω := C ([0 , T ] , R d ) , T > be the canonical space with the canonical process X andthe canonical filtration F = ( F t ) ≤ t ≤ T generated by X . We denote by P the collection of allprobability measures P on (Ω , F T ) under which X ∈ Ω is an ( F , P ) -semi-martingale givenby X t = X + A t + M t , t ∈ [0 , T ] , P -a.s.,where M is an ( F , P ) -martingale with quadratic variation h X t i = h M t i = B t , and theprocesses A and B are P -a.s. absolutely continuous with respect to t . Denote by S d the setof d × d symmetric matrices and S d + ⊂ S d the set of positive semidefinite matrices. For anymatrices A, B ∈ S d , we write A : B := tr( A ⊺ B ) for their scalar product. We say ( α P , β P ) are the characteristics of P if α P t = dA P t dt , β P t = dB P t dt , where ( α P , β P ) takes values in the space R d × S d + . Note that ( α P , β P ) is F -adapted anddetermined up to d P × dt , almost everywhere. Let P ⊂ P be a subset of probabilitymeasures P under which the characteristics ( α P , β P ) are P -integrable on the interval [0 , T ] .In other words, E P Z T | α P | + | β P | dt ! < + ∞ , where | · | is the L -norm.Given a vector τ := ( t , . . . , t m ) ∈ (0 , T ] m , denote by G a vector of m functions suchthat each function G i ∈ C b ( R ) for all i = 1 , . . . , m . Given a Dirac measure µ = δ x and avector c ∈ R m , we define P ( µ , τ, c, G ) ⊂ P as follows: P ( µ , τ, c, G ) := { P : P ∈ P , P ◦ X − = µ and E P [ G i ( X t i )] = c i , i = 1 , . . . , m } . For brevity, we also write E P t,x := E P [ · | X t = x ] .For technical reasons, we restrict ourselves to functions G i in C b ( R ) . In the contextof volatility models calibration, G i are discounted European payoffs. Although the calloption payoff functions are not technically in C b ( R ) , we only work with them in a truncated(compact) space in practice. Alternatively, one may consider only put options using put-callparity. Now let us introduce the formulation of the semi-martingale optimal transport problem underdiscrete constraints. 3efine the cost function F : Λ × R d × S d → R ∪{ + ∞} where F ( t, x, α, β ) = + ∞ if β / ∈ S d + ,and F ( t, x, α, β ) is nonnegative, proper, lower semi-continuous and convex in ( α, β ) . We alsolet F be coercive in the sense that there exist constant p > and C > such that | α | p + | β | p ≤ C (1 + F ( t, x, α, β )) , ∀ ( t, x ) ∈ Λ . Its convex conjugate F ∗ : Λ × R d × S d → R ∪ { + ∞} for ( α, β ) is then defined as F ∗ ( t, x, a, b ) := sup α ∈ R d ,β ∈ S d + { α · a + β : b − F ( t, x, α, β ) } . For simplicity, we write F ( α, β ) := F ( t, x, α, β ) and F ∗ ( a, b ) := F ∗ ( t, x, a, b ) if there is noambiguity.We are interested in the following minimisation problem: Definition 3.1 (Problem 1) . Given µ , τ, c and G , we want to find inf P ∈P ( µ ,τ,c,G ) E P Z T F ( α P t , β P t ) dt. The problem is said to be admissible if P ( µ , τ, c, G ) is nonempty and the infimum above isfinite. The following lemma is an immediate consequence of Itô’s formula and [33, Theorem 2.5](see [12, Theorem 2.6] for a shorter proof in the case of bounded coefficients).
Lemma 3.2.
Let P ∈ P and ρ P be the law of P . Defining α ( t, x ) = E P t,x α P t and β ( t, x ) = E P t,x β P t , then ρ P is a weak solution to the Fokker–Planck equation: ∂ t ρ t ( x ) + ∇ x · ( ρ t ( x ) α ( t, x )) − X i,j ∂ ij ( ρ t ( x ) β ij ( t, x )) = 0 in Λ ,ρ = µ in R d . (1) Moreover, there exists a Markov process X such that, under another probability measure P ′ ∈ P , X has law ρ P ′ = ρ P and solves (cid:26) dX t = α ( t, X t ) dt + σ ( t, X t ) dW P ′ t ,X = x , (2) where σ : Λ → R d × d is defined such that β = σσ ⊺ and W P ′ is a P ′ -Brownian motion. Definition 3.3.
Define P loc to be a subset of P such that, under P ′ ∈ P loc , X is a Markovprocess that has law ρ P ′ = ρ P and takes the form of (2) with coefficients ( α ( t, x ) , β ( t, x )) =( E P t,x α P t , E P t,x β P t ) for some P ∈ P . In the next proposition we write ( E P t,x α P t , E P t,x β P t ) as ( α P ′ loc , β P ′ loc ) for P ∈ P with ρ P ′ = ρ P .The superscript P ′ indicates the dependence of P ′ via ρ P ′ = ρ P . Then, we establish thefollowing result: Proposition 3.4 (Localisation) . Let α P ′ loc = E P t,x α P t and β P ′ loc = E P t,x β P t . Then, inf P ∈P E P Z T F ( α P t , β P t ) dt = inf P ′ ∈P loc E P ′ Z T F ( α P ′ loc , β P ′ loc ) dt. (3)4 roof. By applying Jensen’s inequality together with the tower property of conditional ex-pectation, we have E P Z T F ( α P t , β P t ) dt = E P Z T (cid:0) E P t,x F ( α P t , β P t ) (cid:1) dt ≥ E P Z T F ( E P t,x α P t , E P t,x β P t ) dt. (4)This shows that by localising ( α P , β P ) , one lowers the transportation cost. Denote by ρ P thelaw of P . Applying Lemma 3.2, there exists P ′ ∈ P and a Markov process that has law ρ P ′ = ρ P and takes the form of (2) with coefficients ( α P ′ loc , β P ′ loc ) = ( E P t,x α P t , E P t,x β P t ) . Thus, E P Z T F ( α P t , β P t ) dt ≥ E P ′ Z T F ( α P ′ loc , β P ′ loc ) dt, (5)and E P is replaced by E P ′ on the right-hand side because ρ P = ρ P ′ . Since P loc ⊂ P , takingthe infimum over P ∈ P on the left-hand side of (5) and over P ′ ∈ P loc on the right-handside of (5), we obtain the required result.Proposition 3.4 shows that it suffices to consider only ( α P ′ loc , β P ′ loc ) . Thus Problem 1 canbe studied via PDE methods. Following the Benamou–Brenier formulation of the classicaloptimal transport [3], we introduce the following problem: Proposition 3.5 (Problem 2) . Given µ , τ, c and G , we want to minimise V = inf ρ,α,β Z T Z R d F ( α ( t, x ) , β ( t, x )) dρ ( t, x ) dt, (6) among all ( ρ, α, β ) ∈ C ([0 , T ] , P ( R d ) − w ∗ ) × L ( dρ t dt, R d ) × L ( dρ t dt, S d ) satisfying (in thedistributional sense) ∂ t ρ ( t, x ) + ∇ x · ( ρ ( t, x ) α ( t, x )) − X i,j ∂ ij ( ρ ( t, x ) β ij ( t, x )) = 0 , (7) Z R d G i ( x ) dρ ( t i , x ) = c i , ∀ i = 1 , . . . , m, and ρ (0 , · ) = µ . (8)The interchange of integrals in (6) is justified by Fubini’s theorem. For the weak continuityof measure ρ in time, the reader can refer to [26, Theorem 3]. This section examines the proof of the duality result by closely following [26, Section 3.2] (seealso [5, 21]).
Theorem 3.6.
Assumed that V is finite. Then, V = sup φ,λ m X i =1 λ i c i − Z R m φ (0 , x ) dµ , (9) where the supremum is taken over all ( φ, λ ) ∈ BV ([0 , T ] , C b ( R d )) × R m satisfying ∂ t φ + m X i =1 λ i G i δ t i + F ∗ ( ∇ x φ, ∇ x φ ) ≤ in [0 , T ) × R d , (10) and φ ( T, · ) = 0 . roof. Define measures A := ρα and B := ρβ , then A and B are absolutely continuouswith respect to ρ . Assuming that there exists an admissible solution (¯ ρ, ¯ A = ¯ ρ ¯ α, ¯ B = ¯ ρ ¯ β ) ofProblem 2 satisfying (7, 8), the constraints can be formulated in the following weak form: ∀ φ ∈ C ∞ c (Λ) , Z Λ ∂ t φ ( dρ − d ¯ ρ ) + ∇ x φ · ( d A − d ¯ A ) + 12 ∇ x φ : ( d B − d ¯ B ) = 0 , (11) ∀ λ ∈ R m , Z Λ m X i =1 λ i G i δ t i ( dρ − d ¯ ρ ) = 0 . (12)Thus Problem 2 can be reformulated as the following saddle point problem: V = inf ρ, A , B sup φ,λ (cid:26) Z Λ F (cid:18) d A dρ , d B dρ (cid:19) dρ − m X i =1 λ i G i δ t i ( dρ − d ¯ ρ ) − ∂ t φ ( dρ − d ¯ ρ ) − ∇ x φ · ( d A − d ¯ A ) − ∇ x φ : ( d B − ¯ B ) (cid:27) . Adopting the terminology of [21], we say the triple ( r, a, b ) is represented by ( φ, λ ) if itsatisfies r + ∂ t φ + m X i =1 λ i G i δ t i = 0 ,a + ∇ x φ = 0 ,b + 12 ∇ x φ = 0 . If we choose ( r, a, b ) from C b (Λ , X ) , then the multiplier φ has possible jump discontinuitiesat t = t i due to the presence of the Dirac delta functions. Now, define functionals Φ : C b (Λ , X ) → R ∪ { + ∞} and Ψ : C b (Λ , X ) → R ∪ { + ∞} as follows: Φ( r, a, b ) = (cid:26) if r + F ∗ ( a, b ) ≤ , + ∞ otherwise, Ψ( r, a, b ) = Z Λ r d ¯ ρ + a · d ¯ A + b : d ¯ B if ( r, a, b ) is represented by ( φ, λ ) ∈ BV ([0 , T ] , C b ( R d )) × R m , + ∞ otherwise.Denote by Φ ∗ and Ψ ∗ the convex conjugates of Φ and Ψ , respectively. For Φ , its convexconjugate Φ ∗ : C b (Λ , X ) ∗ → R ∪ { + ∞} is given by Φ ∗ ( ρ, A , B ) = sup ( r,a,b ) ∈ C b (Λ , X ) {h ( r, a, b ) , ( ρ, A , B ) i ; r + F ∗ ( a, b ) ≤ } . As shown in Lemma A.1, if we restrict Φ ∗ to M (Λ , X ) , then Φ ∗ ( ρ, A , B ) = Z Λ F (cid:18) d A dρ , d B dρ (cid:19) dρ if ρ ∈ M + (Λ , X ) and ( A , B ) ≪ ρ, + ∞ otherwise.Next, its convex conjugate Ψ ∗ : C b (Λ , X ) ∗ → R ∪ { + ∞} is given by Ψ ∗ ( ρ, A , B ) = sup r,a,b h ( r, a, b ) , ( ρ − ¯ ρ, A − ¯ A , B − ¯ B ) i , where the supremum is performed over all triples ( r, a, b ) ∈ C b (Λ , X ) represented by ( φ, λ ) in BV ([0 , T ] , C b ( R d )) × R m . In terms of ( φ, λ ) , Ψ ∗ ( ρ, A , B ) = sup φ,λ h ( − ∂ t φ − m X i =1 λ i G i δ t i , −∇ x φ, − ∇ x φ ) , ( ρ − ¯ ρ, A − ¯ A , B − ¯ B ) i . V can be expressed as V = inf ( ρ, A , B ) ∈M (Λ , X ) (Φ ∗ + Ψ ∗ )( ρ, A , B ) = inf ( ρ, A , B ) ∈ C b (Λ , X ) ∗ (Φ ∗ + Ψ ∗ )( ρ, A , B ) , where the second equality is proved in Lemma A.2.Let O m × n denote the null matrix of size m × n . Consider the point ( r, a, b ) = ( − , O d × , O d × d ) which can be represented by ( φ, λ ) = ( − t, O m × ) . As F is nonnegative, at ( − , O d × , O d × d ) we have − F ∗ ( O d × , O d × d ) = − − inf α ∈ R d ,β ∈ S d F ( α, β ) < . This shows that Φ( − , O d × , O d × d ) = 0 , Ψ( − , O d × , O d × d ) = − T. Thus, at ( − , O d × , O d × d ) , Φ is continuous with respect to the uniform norm (since F ∗ iscontinuous in dom( F ∗ ) ), and Ψ is finite. Furthermore, as the convex functionals Φ and Ψ take values in ( −∞ , + ∞ ] , all of the required conditions are fulfilled to apply the Fenchel–Rockafellar duality theorem (see [6, Chapter 1]). We then obtain V = inf ( ρ, A , B ) ∈ C b (Λ , X ) ∗ { Φ ∗ ( ρ, A , B ) + Ψ ∗ ( ρ, A , B ) } = sup ( r,a,b ) ∈ C b (Λ , X ) {− Φ( − r, − a, − b ) − Ψ( r, a, b ) } , and the infimum is in fact attained. Consequently, V = sup r,a,b (cid:26) Z Λ − r d ¯ ρ − a · d ¯ A − b : ¯ B ; − r + F ∗ ( − a, − b ) ≤ (cid:27) , where the supremum is restricted to all ( r, a, b ) represented by ( φ, λ ) ∈ BV ([0 , T ] , C b ( R d )) × R m . Again, in terms of ( φ, λ ) , this is equivalent to V = sup φ,λ Z Λ ( ∂ t φ + m X i =1 λ i G i δ t i ) d ¯ ρ + ∇ x φ · d ¯ A + 12 ∇ x φ : d ¯ B , where ( φ, λ ) ∈ BV ([0 , T ] , C b ( R d )) × R m subject to (10). Recalling that the admissiblesolution (¯ ρ, ¯ A , ¯ B ) satisfies (7, 8) and fixing φ ( t = T ) = 0 , we obtain the required result byintegrating by parts. Adopting the concept of viscosity solutions, it can be shown that the supremum of theobjective with respect to φ is achieved by the viscosity solution of a HJB equation. Inorder to deal with the jump discontinuities in φ , assuming that T is the longest maturity( T = max t k ), we define the viscosity solution in the following way. Definition 3.7.
Denote by set( τ ) the set of entries of vector τ and by K the cardinality of set( τ ) . Let t = 0 , we define disjoint intervals I k := [ t k − , t k ) such that K [ k =1 I k = [0 , T ) , where t k − < t k and t k ∈ set( τ ) for all k = 1 , . . . , K . efinition 3.8 (Viscosity solution) . For any λ ∈ R m , we say φ is a viscosity subsolution(resp,. supersolution) of ∂ t φ + m X i =1 λ i G i δ t i + F ∗ ( ∇ x φ, ∇ x φ ) = 0 in [0 , T ) × R d ,φ T = 0 in R d , (13) if φ is a classical (continuous) viscosity subsolution (resp,. supersolution) of (13) in I k × R d for all k = 1 , . . . , K , and has jump discontinuities: φ ( t, x ) = φ ( t − , x ) − m X i =1 λ i G i ( x ) ( t = t i ) ∀ ( t, x ) ∈ τ × R d . Then, φ is called a viscosity solution of (13) if φ is both a viscosity subsolution and a viscositysupersolution of (13). Remark 3.9 (Comparison principle) . The comparison principle still holds for viscosity so-lutions of (13). Let u and v be a viscosity subsolution and a viscosity supersolution of theequation (13), respectively. At the terminal time T , u ( T, · ) ≤ v ( T, · ) . Since t K = T is in set( τ ) and u, v have the same jump size at { T } × R d , we get u ( T − , · ) ≤ v ( T − , · ) . Next, inthe interval I K = [ t K − , t K ) , by the classical comparison principle, we get u ≤ v on I K .Applying this argument for all intervals I k for k = 1 , . . . , K , we conclude that u ( t, x ) ≤ v ( t, x ) ∀ ( t, x ) ∈ [0 , T ] × R d . Also, u (0 , · ) ≤ v (0 , · ) . Remark 3.10 (Existence and uniqueness) . As a consequence of the comparison principle,there exists a unique viscosity solution of (13). The uniqueness is a direct consequence of thecomparison principle. The existence can be obtained by the Perron’s method (see [7]) underwhich the comparison principle is a key argument.
Now we show that, with any λ ∈ R m , the optimal φ is the viscosity solution of the HJBequation (13). We prove the result by the smoothing technique used in [4]. The proof issimilar to [4, Theorem 2.4], which we sketch here for completeness. Corollary 3.11.
The dual formulation is equivalent to V = sup λ ∈ R m m X i =1 λ i c i − Z R d φ (0 , x ) dµ , (14) where φ is the viscosity solution to the HJB equation (13). Moreover, if the supremum in φ is attained by some ¯ φ , the optimal ( α, β ) is given by ( α, β ) = ∇ F ∗ ( ∇ x ¯ φ, ∇ x ¯ φ ) . (15) Proof.
Denote by ϕ a viscosity solution of the equation (13) with any λ ∈ R m . From Remark3.10, we know that such ϕ exists and is unique. The proof is divided into two steps. Step 1.
Assuming that there exists a sequence of supersolutions of (13) in BV ([0 , T ] , C b ( R d )) converging to ϕ pointwise, we can show that ϕ achieves the supremum with respect to φ inthe objective of the dual (9). Let φ ∈ BV ([0 , T ] , C b ( R d )) be any solution satisfies (10), and φ is also a (viscosity) supersolution of (13). By Remark 3.9, we have ϕ (0 , x ) ≤ φ (0 , x ) forall x ∈ R d , hence m X i =1 λ i c i − Z R d φ (0 , x ) dµ ≤ m X i =1 λ i c i − Z R d ϕ (0 , x ) dµ . (16)8he equality can be achieved in (16) by taking the supremum with respect to φ on the left-hand side of (16). Step 2.
Now, we shall construct the sequence of supersolutions required in Step 1. We firstsmooth the viscosity solution ϕ in space by a standard regularising kernel. By applying theresult of [4] which relies critically on the fact that F ∗ ( a, b ) is convex in ( a, b ) , it can beshown that ϕ ε are supersolutions of equation (13). If we send ε to , the supersolutions ϕ ε converges to the viscosity solution ϕ pointwise. The desired sequence is then constructed.The second part of the theorem is easy to check by the definition of convex conjugate.Therefore, the claim is proven. In this section, we illustrate our method by calibrating a Heston-like LSV model. This methodcould also be easily extended to other LSV models. We consider the model with followingdynamics under the risk-neutral measure: dZ t = ( r ( t ) − q ( t ) − σ ( t, Z t , V t )) dt + σ ( t, Z t , V t ) dW Zt ,dV t = κ ( θ − V t ) dt + ξ √ V t dW Vt ,dW Zt dW Vt = η ( t, Z t , V t ) dt, (17)where Z t is the logarithm of the stock price at time t . The interpretations of r and q differbetween financial markets. In the equity market, r is the risk-free rate and q is the dividendyield. In the FX market, r is the domestic interest rate and q is the foreign interest rate.The parameters κ, θ, ξ have the same interpretation as in the Heston model. In our method,we assume these parameters are obtained by calibrating a pure Heston model. In contrastto the LSV models in other papers, we consider a local-stochastic volatility σ > and alocal-stochastic correlation η ∈ [ − , whose values depend on ( t, Z t , V t ) . Our objective isto calibrate σ ( t, Z, V ) and η ( t, Z, V ) so that model prices exactly match market prices. Remark 4.1.
If the volatility σ ( t, Z, V ) ≡ √ V and the correlation η ( t, Z, V ) is a constant,the LSV model reduces to a pure Heston model. Furthermore, if σ ( t, Z, V ) is independent ofthe variable V , the model is equivalent to a local volatility model. Let the canonical process X t be the vector of processes ( Z t , V t ) with an initial distribution µ = ( δ Z , δ V ) . We want to find a probability measure P ∈ P ( µ , τ, c, G ) characterised by ( α, β ) where α = (cid:20) ( r − q − σ / κ ( θ − V ) (cid:21) and β = (cid:20) σ ηξ √ V σηξ √ V σ ξ V (cid:21) , (18)for the F -adapted processes σ and η . On other words, we want the dynamics of X to be consistent with (17). Given m European options with prices c ∈ R m + , maturities τ = ( t , . . . , t m ) ∈ (0 , T ] m and discounted payoffs G = ( G , . . . , G m ) where G i : R + → R (e.g., G i ( Z ) = e − R ti r ( s ) ds ( e Z − K ) + for a European call of strike K and maturity t i ), thecalibration problem can be formulated as inf P ∈P ( µ ,τ,c,G ) E P Z T F ( t, X t , α P t , β P t ) dt, (19)where F is a suitable convex cost function that forces ( α, β ) to take the form of (18).One possible way to choose the cost function F is based on the idea of minimising the dif-ference between each element of the covariance matrix β and a reference value while keeping β to be in S . However, it is often impossible to find an explicit formula to approximate F ∗ .9hus numerical optimisation is needed, which makes the method computationally expensive.To overcome this issue, we choose the correlation η ( t, Z, V ) = √ Vσ ¯ η, where ¯ η is a constant correlation obtained (along with κ, θ, ξ ) by calibrating a pure Hestonmodel. In this case, β is positive semidefinite if and only if σ ≥ ¯ η V .Figure 1: The plot of H ( x, ¯ x, s ) for given ¯ x and s. Definition 4.2.
Define function H : R × R + × R → R such that H ( x, ¯ x, s ) := ( a ( x − s ¯ x − s ) p + b ( x − s ¯ x − s ) − p + c if x > s and ¯ x > s, + ∞ otherwise.The parameter p is a constant greater than 1, and a, b, c are constants determined to minimisethe function at x = ¯ x with min H = 0 . Given ¯ x and s satisfying ¯ x > s , the function H is convex in x and minimised at ¯ x . It isfinite only when x > s . A plot of H is given in Figure 1. Then, we define the cost functionas follows. Definition 4.3.
The cost function F : R × R × R × R × S → R is defined as F ( t, Z, V, α, β ) := (cid:26) H ( β , V, ¯ η V ) if ( α, β ) ∈ Γ , + ∞ otherwise,where the convex set Γ is defined as Γ := { ( α, β ) ∈ R × S | α = r − q − β / , α = κ ( θ − V ) , β = β = ¯ ηξV, β = ξ V } . Remark 4.4.
The function H penalises deviations of the LSV model from a pure Hestonmodel by choosing ¯ x = V (see Remark 4.1). This approach seeks to retain the attractivefeatures of the Heston model while still matching all the market prices. We also set s = ¯ η V to ensure that β remains positive semidefinite, so the correlation η is in [ − , . The set Γ forces the dynamics of the semi-martingale X t to be consistent with the LSV model (17). Inparticular, it remains risk neutral. V = sup λ ∈ R m J ( λ ) := sup λ ∈ R m m X i =1 λ i c i − φ (0 , Z , V ) , (20)where φ is the viscosity solution to the HJB equation ∂ t φ + m X i =1 λ i G i δ t i + sup β (cid:26) ( r − q − β ) ∂ Z φ + κ ( θ − V ) ∂ V φ + ¯ ηξV ∂ ZV φ + 12 β ∂ ZZ φ + 12 ξ V ∂
V V φ − H ( β , V, ¯ η V ) (cid:27) = 0 , (21)with a terminal condition φ ( T, · ) = 0 .The dual formulation is easier to solve than its primal. Given any λ ∈ R m , we cancalculate φ (0 , Z , V ) by numerically solving the HJB equation (21). The optimal λ can befound by applying a standard optimisation algorithm. The convergence of the algorithm canbe improved by providing the gradient of the objective. Lemma 4.5.
The gradient of the objective with respect to λ i can be formulated as: ∂ λ i J ( λ ) = c i − E P , ( Z ,V ) G i ( Z t i ) ∀ i = 1 , . . . , m. (22) Proof.
Since λ appears in (21), φ (0 , Z , V ) is also a function of λ . Differentiating (21) withrespect to λ i for any i = 1 , . . . , m and writing φ ′ = ∂ λ i φ , we obtain the following PDE ∂ t φ ′ + ( r − q − σ ) ∂ Z φ ′ + κ ( θ − V ) ∂ V φ ′ + ¯ ηξV ∂ ZV φ ′ + 12 σ ∂ ZZ φ ′ + 12 ξ V ∂
V V φ ′ = − G i δ t i , (23)where σ is the optimal β that achieves the supremum in (21). With the terminal condition φ ′ ( T, · ) = 0 , we solve (23) by the Feynman-Kac formula (see [25, Theorem 7.6]). Thus, φ ′ (0 , Z , V ) = E P , ( Z ,V ) G i ( Z t i ) . Differentiating J ( λ ) with respect to λ i completes the proof. Remark 4.6.
Note that E P , ( Z ,V ) G i ( Z t i ) is the model price of the i -th European option.The gradient can be thus interpreted as the difference between market price and model price.Also, solving (23) with terminal condition φ ′ ( T, · ) = 0 is equivalent to solving ∂ t φ ′ + ( r − q − σ ) ∂ Z φ ′ + κ ( θ − V ) ∂ V φ ′ + ¯ ηξV ∂ ZV φ ′ + 12 σ ∂ ZZ φ ′ + 12 ξ V ∂
V V φ ′ = 0 , (24) with terminal condition φ ′ ( t i , · ) = G i . Once the optimal ( φ, λ ) has been found, the optimal volatility σ can be recovered by(15), which is equivalent to solving ( ∂ ZZ φ − ∂ Z φ ) / ∂ σ H ( σ , V, ¯ η V ) , (25)for which a closed-form solution is available. 11 Numerical examples
In this section, we present a numerical method for solving the dual formulation. Let us firstrestrict the problem to a computational domain Q ⊆ R , which should be large enough sothat the density vanishes at the boundary. We consider a uniform mesh of the domain Q anda uniform mesh of the time interval [0 , T ] with step size ∆ t = T /N T with N T ∈ N . Thus t k = k ∆ t, ∀ k = 0 , . . . , N T . If we have the optimal σ t k and any λ ∈ R m , the HJB equation (21) becomes a linearPDE which can be numerically solved by an ADI finite difference method (see [22] for moredetails). To approximate the optimal σ t k , at each time step t = t k , we solve (25) with φ = φ t k +1 via a closed-form formula. If t k +1 equals to the maturity of any input option,we incorporate the jump discontinuity by adding the jump size P mi =1 λ i G i ( t i = t k +1 ) into φ t k +1 .Setting an initial λ , the objective V is maximised over λ ∈ R m by an optimisationalgorithm. At each optimisation iteration, we calculate φ by the above arguments. Theprocess is aided by numerically computing the gradient (22). We measure the optimality bythe uniform norm of the gradient. Setting a threshold ǫ , the algorithm terminates when thefollowing stopping criterion is reached: k ∂ λ J ( λ ) k ∞ ≤ ǫ . The method is summarised in Algorithm 1.
With simulated data, we provide two numerical examples to demonstrate the calibrationresults. In both examples, the risk-free rate is set to a constant r = 0 . and the dividendyield is set to q = 0 . Let Z = ln 100 and V = 0 . for both models, we set the problemin the computational domain Q = [ Z − √ V , Z + 4 √ V ] × [0 , . with the time interval [0 , . For the spatial discretisation on Q , we use 51 points in Z-direction and 51 points inV-direction. For the time interval, we choose N T = 100 . The LSV model is calibrated toa set of European call options generated by the Heston model with given parameters. The Algorithm 1:
LSV calibration
Data:
Market prices of European option
Result:
A calibrated LSV model that matches all market prices Set an initial λ ; do for k = N T − , . . . , do if t k +1 is equal to the maturity of any input option. then φ t k +1 ← φ t k +1 + P mi =1 λ i G i ( t i = t k +1 ) ; end Approximate the optimal σ t k by solving (25) with φ t k +1 ; Solve the HJB equation (21) by the ADI method at t = t k ; end Solve (24) to calculate the model prices by the ADI method; Calculate the gradient ∂ λ J ( λ ) by (22); Update λ by an optimisation algorithm; while k ∂ λ J ( λ ) k ∞ > ǫ ; 12igure 2: The volatility function σ ( t, Z, V ) for Example 1.option prices are calculated at maturities t i ∈ { . , . , . , . , . } for strikes of all pointsin [ Z − . √ V , Z + 1 . √ V ] . The parameter p = 4 is chosen for the cost function. Forsolving the HJB equation (21) and the backward pricing PDE (24), we use the ADI methodof Douglas and Rachford [9] (also see [22]). We use the parameters κ = 0 . , θ = 0 . , ξ = 0 . and ¯ η = − . for both the LSV model andthe Heston model. This example represents a trivial case, since we know that the optimalanalytical solution is σ ( t, Z, V ) = V and η ( t, Z, V ) = ¯ η if we use the same set of parametersfor both models. With a threshold ǫ = 10 − , we obtain the expected results (see Figure 2). In this example, we give different parameters to the LSV model and the Heston model (seeTable 2). As in Remark 4.1, the LSV model can be converted to a LV model, so it hasthe ability to be exactly calibrated to all option prices generated by Heston model with anyparameters. We obtain the result with threshold ǫ = 0 . . Figures 4 and 3 visualise thevolatility function σ ( t, Z, V ) and the correlation function η ( t, Z, V ) . Figure 5 shows theimplied volatility (IV) of the input and the model generated options. We can see that all themodel IVs match their corresponding input IVs. The values of a subset of implied IVs aregiven in Table 1. 13igure 3: The correlation function η ( t, Z, V ) in Example 2. In this section, we calibrate the LSV model to the FX options data provided in [32]. Thedata are listed in Appendix A.4. The parameters we use are shown in Table 3. In this case, κθ/ξ = 0 . ≪ and the Feller condition is strongly violated.For the numerical settings, the computational domain is set to Q = [ − . , . × [0 , .For improving the accuracy while still keeping a reasonable computation time, we employnon-uniform meshes in both Z- and V-directions and put more points around ( Z , V ) (see[22, Section 2.2] for more details). For the time interval [0 , , we give 30 time steps with anequal step size between any two consecutive maturities.In the simulated data example, F ∗ and the optimal σ t k are simply approximated byusing φ t k +1 for a clear illustration. We therefore introduce an iterative algorithm to improvethe convergence. It is described in Appendix A.3.Setting the threshold ǫ = 6 × − , we obtain an exact calibration. The maximumdifference between the model IVs and market IVs is less than 1 basis point. Figure 6 showsthe IVs of short-maturity options (1 month and 3 months) for the market data, the Hestonmodel and the LSV model. Figure 7 shows the IVs of long-maturity options (2 years and 5years). 14 aturity Log-strike Input IV Model IV Error4.3492 0.2396 0.2396 5.05E-054.4452 0.2291 0.2291 2.81E-06T = 0.2 4.5732 0.2199 0.2199 4.08E-064.7012 0.2138 0.2138 5.02E-064.8292 0.2123 0.2124 8.48E-064.3492 0.2488 0.2488 3.94E-064.4452 0.2422 0.2422 1.16E-06T = 0.4 4.5732 0.2359 0.2359 1.65E-074.7012 0.2303 0.2303 9.82E-074.8292 0.2257 0.2257 4.55E-064.3492 0.2576 0.2576 2.84E-064.4452 0.2523 0.2523 1.60E-06T = 0.6 4.5732 0.2471 0.2471 7.56E-074.7012 0.2423 0.2423 3.35E-074.8292 0.2378 0.2378 3.24E-094.3492 0.2646 0.2646 3.45E-064.4452 0.2600 0.2600 8.69E-07T = 0.8 4.5732 0.2555 0.2555 7.62E-074.7012 0.2512 0.2512 1.42E-074.8292 0.2472 0.2472 4.91E-094.3492 0.2699 0.2699 8.28E-074.4452 0.2659 0.2659 1.83E-06T = 1.0 4.5732 0.2620 0.2620 1.17E-064.7012 0.2581 0.2581 1.54E-064.8292 0.2544 0.2544 9.24E-07 Table 1: A subset of the input IVs and model IVs in Example 2. κ θ ξ ¯ η Heston 2.0 0.09 0.10 -0.6LSV 0.5 0.04 0.16 -0.4Table 2: The parameters used in the Heston model and the LSV modelParameter κ θ ξ ¯ η Z V Value 0.8721 0.0276 0.5338 -0.3566 0.2287 0.012Table 3: The parameters used in the market data example
A Appendix
A.1 Lemma A.1
Lemma A.1.
Define
Φ : C b (Λ , X ) → R ∪ { + ∞} by Φ( r, a, b ) = (cid:26) if r + F ∗ ( a, b ) ≤ , + ∞ otherwise. σ ( t, Z, V ) in Example 2. If we restrict the domain of its convex conjugate Φ ∗ : C b (Λ , X ) ∗ → R ∪ { + ∞} to M (Λ , X ) ,then Φ ∗ ( ρ, A , B ) = Z Λ F (cid:18) d A dρ , d B dρ (cid:19) dρ if ρ ∈ M + (Λ) and ( A , B ) ≪ ρ, + ∞ otherwise.Proof. For any ( ρ, A , B ) ∈ C b (Λ , X ) ∗ , the convex conjugate Φ ∗ is given by Φ ∗ ( ρ, A , B ) = sup ( r,a,b ) ∈ C b (Λ , X ) {h ( r, a, b ) , ( ρ, A , B ) i ; r + F ∗ ( a, b ) ≤ } . If we restrict the domain of Φ ∗ to M (Λ , X ) ⊂ C b (Λ , X ) ∗ , then Φ ∗ ( ρ, A , B ) = sup ( r,a,b ) ∈ C b (Λ , X ) (cid:26) Z Λ r dρ + a · d A + b : d B ; r + F ∗ ( a, b ) ≤ (cid:27) . To show that ρ is a nonnegative measure, we assume that ρ ( E ) < for some measurableset E ⊂ Λ . By the density of C b in L , there exists a sequence of functions in C b ( R ) thatconverges to E . Let r = − λ E for some positive λ , a = O d × and b = O d × d where O m × n denotes the null matrix of size m × n . It is clear that the constraint r + F ∗ ( a, b ) ≤ issatisfied at ( − λ E , O d × , O d × d ) as F is nonnegative. If we send λ to infinity, the function Φ ∗ becomes unbounded.To show ( A , B ) ≪ ρ , we assume that there exists a measurable set E such that ( A , B )( E ) =0 but ρ ( E ) = 0 . Again, we construct a sequence of continuous functions in C b (Λ) convergingto E in L ( dρ ) . Then let r = − F ∗ ( a, b ) , a = λ E I d × and b = λ E I d × d where I m × n denotes the identity matrix of size m × n . The function Φ ∗ goes to infinity if we send λ to + ∞ or −∞ , depending on the sign of A and B on E .16 .35 4.4 4.45 4.5 4.55 4.6 4.65 4.7 4.75 4.8 4.85 Log Strike I m p li ed V o l a t ili t y ( I V ) Implied Volatility
Maturity T = 1.0Maturity T = 0.8Maturity T = 0.6Maturity T = 0.4Maturity T = 0.2Input IV
Figure 5: The input and the model generated IV for Example 2.Now, since the integrand is linear in ( r, a, b ) , if Φ ∗ is finite, the supremum must occur atthe boundary. Thus, assuming that ρ ∈ M + (Λ , X ) and ( A , B ) ≪ ρ , we have Φ ∗ ( ρ, A , B ) = sup r + F ∗ ( a,b )=0 Z Λ (cid:18) r + a · d A dρ + b : d B dρ (cid:19) dρ = sup a,b Z Λ (cid:18) a · d A dρ + b : d B dρ − F ∗ ( a, b ) (cid:19) dρ = Z Λ sup a,b (cid:18) a · d A dρ + b : d B dρ − F ∗ ( a, b ) (cid:19) dρ = Z Λ F (cid:18) d A dρ , d B dρ (cid:19) dρ. The interchange of the supremum and integral is made by choosing a sequence of continuousfunctions ( a, b ) n converging to ∇ F ( d A dρ , d B dρ ) , then applying the dominated convergence theo-rem and the fact that F ∗ is continuous in its effective domain. The last equality holds sincethe cost function F equals to its biconjugate F ∗∗ according to the Fenchel–Moreau theorem(see [6, Theorem 1.11]). The proof is completed. A.2 Lemma A.2
In this section, we prove that, the duality between spaces C b and M can be extended to thenon-compact space [0 , T ] × R d in this particular case. A similar argument for the Kantorovichduality of the classical optimal transport was made in [34, Appendix 1.3]. Lemma A.2.
Denote by K o the set of ( r, a, b ) in C b (Λ , X ) that can be represented bysome ( φ, λ ) in BV ([0 , T ] , C b ( R d )) × R m (see the proof of Theorem 3.5 for the definition of Strike I m p li ed V o l a t ili t y ( I V ) IV for T = 1m
LSVHestonMarket
Strike I m p li ed V o l a t ili t y ( I V ) IV for T = 3m
LSVHestonMarket
Figure 6: IV for options with maturities 1 month and 3 months
Strike I m p li ed V o l a t ili t y ( I V ) IV for T = 2Y
LSVHestonMarket
Strike I m p li ed V o l a t ili t y ( I V ) IV for T = 5Y
LSVHestonMarket
Figure 7: IV for options with maturities 2 years and 5 years ’represented’). Let Φ ∗ : C b (Λ , X ) ∗ → R ∪ { + ∞} and Ψ ∗ : C b (Λ , X ) ∗ → R ∪ { + ∞} be definedby Φ ∗ ( ρ, A , B ) = sup ( r,a,b ) ∈ C b (Λ , X ) {h ( r, a, b ) , ( ρ, A , B ) i ; r + F ∗ ( a, b ) ≤ } , Ψ ∗ ( ρ, A , B ) = sup ( r,a,b ) ∈ K o h ( r, a, b ) , ( ρ − ¯ ρ, A − ¯ A , B − ¯ B ) i . Then, inf ( ρ, A , B ) ∈ C b (Λ , X ) ∗ (Φ ∗ + Ψ ∗ )( ρ, A , B ) = inf ( ρ, A , B ) ∈M (Λ , X ) (Φ ∗ + Ψ ∗ )( ρ, A , B ) . (26) Proof.
Let C (Λ , X ) be the space of continuous functions on Λ valued in X that vanish atinfinity. We decompose ( ρ, A , B ) = (˜ ρ, ˜ A , ˜ B ) + ( δρ, δ A , δ B ) such that (˜ ρ, ˜ A , ˜ B ) ∈ M (Λ , X ) and h ( φ ρ , φ A , φ B ) , ( δρ, δ A , δ B ) i = 0 for any ( φ ρ , φ A , φ B ) ∈ C (Λ , X ) . Since, M (Λ , X ) is asubset of C b (Λ , X ) ∗ , it follows that inf ( ρ, A , B ) ∈ C b (Λ , X ) ∗ (Φ ∗ + Ψ ∗ )( ρ, A , B ) ≤ inf ( ρ, A , B ) ∈M (Λ , X ) (Φ ∗ + Ψ ∗ )( ρ, A , B ) . Next, we show that the converse of the above inequality is also valid. The reader can refer to [34, Appendix 1.3] for the existence of such a decomposition. Φ ∗ , it follows by a standard approximation argument that Φ ∗ ( ρ, A , B ) ≥ sup ( r,a,b ) ∈ C (Λ , X ) {h ( r, a, b ) , ( ρ, A , B ) i ; r + F ∗ ( a, b ) ≤ } = sup ( r,a,b ) ∈ C (Λ , X ) (cid:26) Z Λ r d ˜ ρ + a · d ˜ A + b : d ˜ B ; r + F ∗ ( a, b ) ≤ (cid:27) = Φ ∗ (˜ ρ, ˜ A , ˜ B ) . For Ψ ∗ , if we restrict its domain to (˜ ρ, ˜ A , ˜ B ) ∈ M (Λ , X ) , then Ψ ∗ = 0 if (˜ ρ, ˜ A , ˜ B ) satisfiesthe constraints (7) and (8) or Ψ ∗ = + ∞ otherwise. Whenever Ψ ∗ is finite, Z Λ r d (˜ ρ − ¯ ρ ) + a · d ( ˜ A − ¯ A ) + b : d ( ˜ B − ¯ B ) = 0 ∀ ( r, a, b ) ∈ K o . (27)The equation (27) holds in particular for ( r, a, b ) in the subset K o ∩ C (Λ , X ) . Thus, Ψ ∗ (˜ ρ, ˜ A , ˜ B ) = sup ( r,a,b ) ∈ K o Z Λ r d (˜ ρ − ¯ ρ ) + a · d ( ˜ A − ¯ A ) + b : d ( ˜ B − ¯ B )= sup ( r,a,b ) ∈ K o ∩ C (Λ , X ) Z Λ r d (˜ ρ − ¯ ρ ) + a · d ( ˜ A − ¯ A ) + b : d ( ˜ B − ¯ B )= sup ( r,a,b ) ∈ K o ∩ C (Λ , X ) h ( r, a, b ) , ( ρ − ¯ ρ, A − ¯ A , B − ¯ B ) i≤ Ψ ∗ ( ρ, A , B ) where the last inequality holds since K o ∩ C (Λ , X ) is a subset of K o .Therefore, inf ( ρ, A , B ) ∈ C b (Λ , X ) ∗ (Φ ∗ + Ψ ∗ )( ρ, A , B ) ≥ inf ( ρ, A , B ) ∈M (Λ , X ) (Φ ∗ + Ψ ∗ )( ρ, A , B ) . This completes the proof.
A.3 Algorithm 2
To improve the convergence, we add an iteration step to handle the nonlinearity in the HJBequation. Let φ oldt k = φ t k +1 , we approximate σ t k using φ oldt k . Then, we solve the HJB equation(21) with σ t k at t = t k to get φ newt k . If k φ newt k − φ oldt k k is greater than a threshold ǫ , welet φ oldt k = φ newt k and repeat above steps. Otherwise, we let φ t k = φ newt k and continue. Themodified algorithm is described in Algorithm 2.19 lgorithm 2: LSV calibration with nonlinear iterations
Data:
Market prices of European option
Result:
A calibrated LSV model that matches all market prices Set an initial λ ; do for k = N T − , . . . , do if t k +1 is equal to the maturity of any input option. then φ t k +1 ← φ t k +1 + P mi =1 λ i G i ( t i = t k +1 ) ; end Let φ oldt k = φ t k +1 ; do Approximate the optimal σ t k by solving (25) with φ oldt k ; Solve the HJB equation (21) by the ADI method at t = t k , and set thesolution to φ newt k ; while k φ newt k − φ oldt k k > ǫ ; Let φ t k = φ newt k ; end Solve (24) to calculate the model prices by the ADI method; Calculate the gradient ∂ λ J ( λ ) by (22); Update λ by an optimisation algorithm; while k ∂ λ J ( λ ) k ∞ > ǫ ; A.4 FX option data
The market data used for calibration is EUR/USD option data quoted on 23 August 2012.The spot price S = 1 . USD per EUR. Maturities, strikes and implied volatilities areshown in Table 4. At each maturity, the options correspond to 10-delta calls, 25-delta calls,50-delta calls, 25-delta puts and 10-delta puts. The interest yields are given in Table 5. Notethat the yields data are yield to maturity. Thus we need to convert them to spot rates forour method. 20aturity Option type Strike IV Maturity Option type Strike IVCall 1.3006 0.0905 Call 1.4563 0.1069Call 1.2800 0.0898 Call 1.3627 0.10521m Call 1.2578 0.0915 1Y Call 1.2715 0.1118Put 1.2344 0.0966 Put 1.1701 0.1278Put 1.2110 0.1027 Put 1.0565 0.1491Call 1.3191 0.0897 Call 1.5691 0.1100Call 1.2901 0.0896 Call 1.4265 0.10962m Call 1.2588 0.0933 2Y Call 1.2889 0.1168Put 1.2243 0.1014 Put 1.1421 0.1328Put 1.1882 0.1109 Put 0.9863 0.1540Call 1.3355 0.0912 Call 1.6683 0.1109Call 1.2987 0.0908 Call 1.4860 0.11223m Call 1.2598 0.0955 3Y Call 1.3113 0.1200Put 1.2160 0.1058 Put 1.1308 0.1352Put 1.1684 0.1185 Put 0.9468 0.1547Call 1.3775 0.0960 Call 1.7507 0.1104Call 1.3213 0.0953 Call 1.5351 0.11276m Call 1.2633 0.1013 4Y Call 1.3306 0.1210Put 1.1973 0.1145 Put 1.1226 0.1365Put 1.1236 0.1316 Put 0.9152 0.1554Call 1.4068 0.1013 Call 1.8355 0.1111Call 1.3329 0.1005 Call 1.5835 0.11379m Call 1.2583 0.1068 5Y Call 1.3505 0.1220Put 1.1745 0.1215 Put 1.1180 0.1379Put 1.0805 0.1407 Put 0.8887 0.1571Table 4: The EUR/USD option data quoted on 23 August 2012.Maturity 1m 2m 3m 6m 9m 1Y 2Y 3Y 4Y 5YDomestic yield 0.41 0.51 0.66 0.95 1.19 1.16 0.60 0.72 0.72 0.72Foreign yield 0.04 0.11 0.23 0.47 1.62 0.64 0.03 0.03 0.03 0.03Table 5: The domestic and foreign yields (in %) quoted on 23 August 2012.
References [1]
Abergel, F., and Tachet, R.
A nonlinear partial integro-differential equation frommathematical finance.
Discrete Contin. Dyn. Syst. 27 , 3 (2010), 907–917.[2]
Avellaneda, M., Friedman, C., Holmes, R., and Samperi, D.
Calibrating volatil-ity surfaces via relative-entropy minimization.
Appl. Math. Finance 4 , 1 (1997), 37–64.[3]
Benamou, J.-D., and Brenier, Y.
A computational fluid mechanics solution to theMonge–Kantorovich mass transfer problem.
Numer. Math. 84 , 3 (2000), 375–393.[4]
Bouchard, B., Loeper, G., and Zou, Y.
Hedging of covered options with linearmarket impact and gamma constraint.
SIAM J. Control Optim. 55 , 5 (2017), 3319–3348.[5]
Brenier, Y.
Minimal geodesics on groups of volume-preserving maps and generalizedsolutions of the Euler equations.
Comm. Pure Appl. Math. 52 , 4 (1999), 411–452.216]
Brezis, H.
Functional analysis, Sobolev spaces and partial differential equations . Uni-versitext. Springer, 2011.[7]
Crandall, M. G., Ishii, H., and Lions, P.-L.
User’s guide to viscosity solutions ofsecond order partial differential equations.
Bull. Amer. Math. Soc. (N.S.) 27 , 1 (1992),1–67.[8]
Dolinsky, Y., and Soner, H. M.
Martingale optimal transport and robust hedgingin continuous time.
Probab. Theory Related Fields 160 , 1-2 (2014), 391–427.[9]
Douglas, Jr., J., and Rachford, Jr., H. H.
On the numerical solution of heatconduction problems in two and three space variables.
Trans. Amer. Math. Soc. 82 (1956), 421–439.[10]
Dupire, B.
Pricing with a smile.
Risk Magazine (1994), 18–20.[11]
Engelmann, B., Koster, F., and Oeltz, D.
Calibration of the Heston stochas-tic local volatility model: A finite volume scheme. https://ssrn.com/abstract=1823769 (2011).[12]
Figalli, A.
Existence and uniqueness of martingale solutions for SDEs with rough ordegenerate coefficients.
J. Funct. Anal. 254 , 1 (2008), 109–153.[13]
Gatheral, J.
The volatility surface: a practitioner’s guide , vol. 357. John Wiley &Sons, 2011.[14]
Guo, I., and Loeper, G.
Path dependent optimal transport and model calibration onexotic derivatives. arXiv preprint arXiv:1812.03526 (2018).[15]
Guo, I., Loeper, G., and Wang, S.
Local volatility calibration by optimal transport. (2019), 51–64.[16]
Guyon, J., and Henry-Labordere, P.
The smile calibration problem solved. https://ssrn.com/abstract=1885032 (2011).[17]
Gyöngy, I.
Mimicking the one-dimensional marginal distributions of processes havingan Itô differential.
Probab. Theory Relat. Fields 71 , 4 (1986), 501–516.[18]
Henry-Labordere, P.
Calibration of local stochastic volatility models to marketsmiles: A monte-carlo approach.
Risk Magazine (2009).[19]
Henry-Labordère, P., and Touzi, N.
An explicit martingale version of the one-dimensional brenier theorem.
Finance Stoch. 20 , 3 (2016), 635–668.[20]
Heston, S. L.
A closed-form solution for options with stochastic volatility with appli-cations to bond and currency options.
Rev. Financ. Stud. 6 , 2 (1993), 327–343.[21]
Huesmann, M., and Trevisan, D.
A Bnamou–Brenier formulation of martingaleoptimal transport. arXiv preprint arXiv:1707.01493 (2017).[22]
In ’t Hout, K. J., and Foulon, S.
ADI finite difference schemes for option pricingin the Heston model with correlation.
Int. J. Numer. Anal. Model. 7 , 2 (2010), 303–320.[23]
Jex, M., Henderson, R., and Wang, D.
Pricing exotics under the smile.
RiskMagazine (1999), 72–75.[24]
Kantorovich, L. V.
On a problem of Monge (in Russian).
Zap. Nauchn. Sem. S.-Peterburg. Otdel. Mat. Inst. Steklov. (POMI) 3 (1948), 255–226.2225]
Karatzas, I., and Shreve, S. E.
Brownian motion and stochastic calculus , second ed.,vol. 113 of
Graduate Texts in Mathematics . Springer-Verlag, 1991.[26]
Loeper, G.
The reconstruction problem for the Euler-Poisson system in cosmology.
Arch. Ration. Mech. Anal. 179 , 2 (2006), 153–216.[27]
Monge, G.
Mémoire sur la théorie des déblais et des remblais.
Histoire de l’AcadémieRoyale des Sciences de Paris (1781).[28]
Pal, S., Wong, T.-K. L., et al.
Exponentially concave functions and a new infor-mation geometry.
Ann. Probab. 46 , 2 (2018), 1070–1113.[29]
Ren, Y., Madan, D., and Qian, M. Q.
Calibrating and pricing with embedded localvolatility models.
Risk Magazine 20 , 9 (2007), 138.[30]
Saporito, Y. F., Yang, X., and Zubelli, J. P.
The calibration of stochastic local-volatility models: an inverse problem perspective.
Comput. Math. Appl. 77 , 12 (2019),3054–3067.[31]
Tan, X., and Touzi, N.
Optimal transportation under controlled stochastic dynamics.
Ann. Probab. 41 , 5 (2013), 3201–3240.[32]
Tian, Y., Zhu, Z., Lee, G., Klebaner, F., and Hamza, K.
Calibrating and pricingwith a stochastic-local volatility model.
Journal of Derivatives 22 , 3 (2015), 21.[33]
Trevisan, D.
Well-posedness of multidimensional diffusion processes with weakly dif-ferentiable coefficients.
Electron. J. Probab. 21 (2016), Paper No. 22, 41.[34]
Villani, C.
Topics in optimal transportation , vol. 58 of
Graduate Studies in Mathe-matics . American Mathematical Society, Providence, RI, 2003.[35]
Wyns, M., and Du Toit, J.
A finite volume–alternating direction implicit approachfor the calibration of stochastic local volatility models.