{'fraction_non_alphanumeric': 0.05856133251224304, 'fraction_numerical': 0.01656842907832615, 'mean_word_length': 3.7595254710397765, 'pattern_counts': {'":': 0, '<': 12, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 14, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Calibration of Local-Stochastic Volatility Models by Optimal Transport

In this paper, we study a semi-martingale optimal transport problem and its application to the calibration of Local-Stochastic Volatility (LSV) models. Rather than considering 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 a PDE formulation along with its dual counterpart. Then we solve numerically the dual problem, which involves a fully non-linear Hamilton-Jacobi-Bellman equation. The method is tested by calibrating a Heston-like LSV model with simulated data and foreign exchange market data.

to optimal transport. The calibration of these new SV models is achieved by modifying only the drift and leaving the volatility of volatility unchanged. Apart from continuous models, as an extension of De March and Henry-Labordere (2019), Guyon (2020) constructed a discrete-time model which solves the challenging joint calibration problem of SPX and VIX options. Recently in Guo et al. (2020), we propose a continuous-time model to solve also the joint calibration problem by extending the approach of this paper.
In this paper, we further extend the approach of Guo et al. (2019) and Guo and Loeper (2018) to the calibration of LSV models. The calibration problem is formulated as a semi-martingale optimal transport problem. Unlike Tan and Touzi (2013), we consider a finite number of discrete constraints given by the prices of European claims. As a consequence of Jensen's inequality, we show that an optimal diffusion process can be chosen to be Markovian in the state variables given by the initial SV model. This result leads to a PDE formulation. By following the duality theory of optimal transport introduced in Brenier (1999) and a smoothing argument used in Bouchard et al. (2017), we establish a dual formulation. We also provide a numerical method to solve a fully non-linear Hamilton-Jacobi-Bellman (HJB) equation arising in the dual formulation. Finally, numerical examples show that the model can be fully calibrated to the European options with both simulated data and FX market data.
Despite its accuracy, our method is quite demanding in terms of computational power. The gradient descent demands at each step to solve one non-linear 2 d PDE, and the computation of the gradient requires one (linear) 2 d PDE per instrument. The most costly operation of numerically solving a linear PDE is inverting a large sparse matrix. However, this operation only needs to be carried once per time step because the computations of all components of the gradient are computed by solving the same linear PDE but with different terminal conditions. Alternatively, all gradients can be efficiently computed in one Monte Carlo simulation, which is a choice we did not make here for the sake of accuracy. Going into higher dimensions (a multi-factor stochastic volatility model for example) would require to increase the dimension of the PDEs, which is problematic as soon as d ≥ 3 . When the goal is only to solve the usual LSV calibration problem, i.e., to find the leverage function, other methods achieve the result faster. For example, for a one-factor model, the PDE method of Ren et al. (2007) only requires solving a two-dimensional non-linear PDE once, and the particle method of Guyon and Henry-Labordère (2012) is even faster and can be applied to high dimensional cases (e.g., calibrating an LSV model with multiple stochastic factors). On the other hand, with the technique developed in this paper, one can fit pathdependent products (Guo and Loeper, 2018), SPX and VIX options (Guo et al., 2020) and here LSV models. Therefore the interest of our method is clearly its broad range of applications, at the cost of a relatively heavy computational cost. We also believe that with the recent developments of numerical methods for solving non-linear PDEs in high dimensions (see e.g., E et al., 2020), our method can be greatly improved in terms of computational speed, and become applicable in high dimensions. Also notice that, being based on gradient descent, for a slight update of the market data, only a few gradient iterations should be needed to update the model. Finally, our method provides a rigorous existence result of an LSV type model. Previous works by Abergel and Tachet (2010) only provide an existence result for small times (see also Jourdain and Zhou (2020); Lacker et al. (2020)).
The paper is organised as follows: In Section 2, we introduce some preliminary definitions. In Section 3, we show the connection between the semi-martingale optimal transport problem and a PDE formulation. Duality results are then established for the PDE formulation. In Section 4, we demonstrate the calibration method using a Heston-like LSV model. Numerical method and results with both simulated data and FX market data are provided in Section 5.

Preliminaries
Given a Polish space E equipped with its Borel σ -algebra, let C(E) be the space of continuous functions on E and C b (E) be the space of bounded continuous functions. Denote by 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 topological dual of C b (E) is given by Let P(E) be the space of Borel probability measures, BV (E) be the space of functions of bounded variation and L 1 (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 1 (dµ, R d ) as the vectorvalued versions of their corresponding spaces. If µ t (x) = µ(t, x) is a measure defined on [0, T ] × R d , we will write dµ or dµ t dt in short for µ(t, dx)dt. Denote by S d the set of d × d symmetric matrices and S d + ⊂ S d the set of positive semidefinite matrices. For any matrices A, B ∈ S d , we write A : B := tr(A ⊺ B) for their scalar product. For convenience, let Λ = [0, T ] × R d and X = R × R d × S d . We use the notation ·, · to denote the duality bracket between C b (Λ, X ) and C b (Λ, X ) * .
Let Ω := C([0, T ], R d ), T > 0 be the canonical space with the canonical process X and the canonical filtration F = (F t ) 0≤t≤T generated by X . We denote by P the collection of all probability measures P on (Ω, F T ) under which X ∈ Ω is an (F, P)-semi-martingale given by where M is an (F, P)-martingale with quadratic variation X t = M t = B t , and the processes A and B are P-a.s. absolutely continuous with respect to t. We say P is characterised by (α P , β P ) if where (α P , β P ) take values in the space R d × S d + . Note that (α P , β P ) is F-adapted and determined up to dP × dt, almost everywhere. Let P 1 ⊂ P be the subset of probability measures P under which the characteristics (α P , β P ) are P-integrable on the interval [0, T ]. In other words, where | · | is the L 1 -norm.
For technical reasons, we restrict ourselves to functions G i in C b (R d ). In the context of volatility models calibration, G i are discounted European payoffs. Although the call option payoff functions are not technically in C b (R d ), we only work with them in a truncated (compact) space in practice. Alternatively, one may consider only put options using put-call parity. It is possible to relax the assumption G i ∈ C b (R d ), but it would require a different set up in topological spaces.

Formulations
In this section, we first formulate the semi-martingale optimal transport problem under discrete constraints. Then a PDE formulation is introduced along with its dual counterpart.
Define the cost function F : and F (t, x, α, β) is nonnegative, proper, lower semi-continuous, strongly convex and coercive in (α, β) and uniformly in (t, x). By F being strongly convex in (α, β) we mean that there exists a constant C > 0 such that for all t, x, α, β, α ′ , β ′ and any subderivative ∇F , where ∇ is performed over (α, β), if where · denotes the Euclidean norm on R d and S d . By F being coercive in (α, β) we mean that there exist constants p > 1 and C > 0 for all t, x, α, β such that We remark that S d in (1) can be replaced by S d + due to the assumption that F (t, x, α, β) is finite only if β ∈ S d + . For simplicity, we write F (α, β) := F (t, x, α, β) and F * (a, b) := F * (t, x, a, b) if there is no ambiguity. Note that our definition of strongly convex does not require F to be differentiable, since only subderivatives are used. Nevertheless, it implies that F is strictly convex and thus F * is differentiable. In addition, the coercivity of F implies that F * is finite.
Adopting the convention inf ∅ = +∞, we are interested in the following minimisation problem: Problem 1. Given µ 0 , τ, c and G, we want to find The problem is said to be admissible if P(µ 0 , τ, c, G) is nonempty and the infimum above is finite.
It is well known that the marginal distributions of diffusion processes at fixed times solve the Fokker-Planck equation in the weak sense. The converse result was given by Figalli (2008) and Trevisan (2016). For brevity, we write E P t,x := E P [ · | X t = x]. As an immediate consequence of Itô's formula and Theorem 2.5 in Trevisan (2016), we introduce the following lemma.
Lemma 3.1. Let P ∈ P 1 and ρ P t = ρ P (t, ·) = P • X −1 t be the marginal distribution of X t under P, t ≤ T . Then ρ P is a weak solution to the Fokker-Planck equation: Moreover, there exists another probability measure P ′ ∈ P 1 , characterised by (α P ′ , β P ′ ), under which X has the same marginals, ρ P ′ = ρ P , and is a Markov process solving where W P ′ is a P ′ -Brownian motion, α P ′ (t, X t ) = E P t,Xt α P t and β P ′ (t, X t ) = E P t,Xt β P t .
The above lemma provides a solution to study semi-martingales via Markov processes in the form of (3). It is worth noting that the idea of using diffusion processes to mimic an Itô process by matching their marginals at fixed times traces back to the classical mimicking theorem of Gyöngy (1986). The uniform ellipticity condition of Gyöngy's mimicking theorem was later relaxed by Brunick and Shreve (2013). In fact, if X is an Itô process under P, Lemma 3.1 can be seen as a reformulation of Brunick and Shreve (2013, Corollary 3.7) which was constructed by a completely different approach. The Markov processes X in (3) are also called Markovian projections in the literature. Note that, in Brunick and Shreve (2013), even though the main results are given for Itô processes, the authors first provide more general results for semi-martingales (see Brunick and Shreve (2013, Theorem 7.1)) and then prove the main results for Itô processes by the Itô representation theorem. Therefore, Lemma 3.1 can also be proved by the results of Brunick and Shreve (2013).
If P(µ 0 , τ, c, G) is not empty, by Lemma 3.3, P loc (µ 0 , τ, c, G) is not empty. For any P ∈ P(µ 0 , τ, c, G), let P ′ ∈ P loc (µ 0 , τ, c, G) be a probability measure such that X has the same marginals under P and P ′ . Applying Jensen's inequality together with the tower property of conditional expectation, we have The last E P is replaced by E P ′ because the marginal of X is the same under P and P ′ . Since P loc (µ 0 , τ, c, G) ⊂ P(µ 0 , τ, c, G), taking infimum over all P ∈ P(µ 0 , τ, c, G) on the left-hand side and over all P ′ ∈ P loc (µ 0 , τ, c, G) on the right-hand side of (5), we obtain the required result.
Proposition 3.4 shows that it suffices to consider only the probability measures in P loc (µ 0 , τ, c, G).
Thus, by the connections established in Lemma 3.1, Problem 1 can be studied via PDE methods. Following the Benamou-Brenier formulation of the classical optimal transport from Benamou and Brenier (2000), we introduce the following problem: Problem 2 (PDE formulation). Given µ 0 , τ, c and G, we want to solve The interchange of integrals in (6) is justified by Fubini's theorem as F is nonnegative. For the weak continuity of measure ρ in time, the reader can refer to Loeper (2006, Theorem 3).
Based on the results of Sections 3.2 and 3.3 below, we shall introduce a dual formulation of Problem 2. In the proposition below, C 2 b (R d ) is the space of twice continuously differentiable functions with bounded partial derivatives up to order 2, and it is equipped with the norm given by the supremum of all partial derivatives up to order 2. The subscript of φ λ indicates the implicit dependence of φ on λ via the HJB equation. The definition of the viscosity solution to (10) and the proof will be given in Section 3.3.
Proposition 3.5 (Dual formulation). If Problem 1 is admissible, then where φ is the viscosity solution to the HJB equation (7) and (8) (in the distributional sense), then the infimum of Problem 2 is attained. If the supremum is attained by some λ * ∈ R m and (ρ, α, β) is an optimal solution of Problem 2, then (α, β) is given by Before ending this section, it is worth commenting on the admissibility of Problem 1. We have chosen to impose the admissibility assumption in order to simplify our presentation and arguments. With some modifications, it is possible to remove this assumption from the primal problem and still obtain duality. In particular, both sides of the duality would be infinite if the problem is not admissible.
Then, characterising the admissibility of Problem 1 corresponds to checking the finiteness of the dual problem, and can be seen as a more elaborate analogue of Strassen's theorem for the classical optimal transport problem.
Theorem 3.6. If Problem 1 is admissible, then where the supremum is taken over all and satisfying (7) and (8) (in the distributional sense), then the infimum of Problem 2 is attained.
Proof. The proof relies on the Fenchel-Rockafellar theorem which plays a key role in the applications of convex analysis. One may note that the objective function (6) is not convex in (ρ, α, β). As we will see below, (6) can be written as the convex conjugate (which is always convex) of another function with respect to (ρ,Ā := αρ,B := βρ) and (Ā,B) are absolutely continuous with respect to ρ. In addition, the constraints (7) and (8) are linear in (ρ,Ā,B). Therefore, throughout the proof, we will work on (ρ,Ā,B) instead. For simplicity, we will write dĀ and dB in short for α(t, x)ρ(t, dx)dt and β(t, x)ρ(t, dx)dt, respectively.
Formulate the constraints (7) and (8) in the following weak form: where C ∞ c (Λ) is the space of smooth functions with compact support on Λ . Thus Problem 2 can be reformulated as the following saddle point problem: The strategy of the proof is to first construct a function Φ whose convex conjugate Φ * is equal to the objective function of Problem 2, and construct another function Ψ whose convex conjugate Ψ * is equal to the rest part inside the infimum of (16) so that V = inf ρ,Ā,B (Φ * + Ψ * )(ρ,Ā,B). Then, the duality is established by applying the Fenchel-Rockafellar theorem.
Adopting the terminology of Huesmann and Trevisan (2019), we say the triple (r, a, b) is represented If we choose (r, a, b) from C b (Λ, X ), by the first equation above, ∂ t φ is a measure because of the presence of the Dirac delta functions. Thus, φ has bounded variation with respect to t on [0, T ] and has possible jump discontinuities at t = τ i . Now, define functionals Φ : (14) and (15), otherwise we can arbitrarily scale (φ, λ) in (16) then V becomes unbounded. Assume that (r, a, b) can be represented by both (φ,λ) and (φ,λ), then we have Integrating (17) with any ρ that satisfies (15) and ρ(0, so the value of Ψ does not depend on the choice of (φ, λ) and hence Ψ is well-defined.
Denote by Φ * and Ψ * the convex conjugates of Φ and Ψ , respectively. For Φ, its convex conjugate

As shown in Lemma
where the supremum is taken over all triples (r, a, b) In terms of (φ, λ), As proved in Lemma A.2, the objective V can be expressed as This shows that , 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 e.g., Brezis, 2011, Chapter 1). We then obtain and the infimum is in fact attained. Consequently, Writing (r, a, b) in terms of (φ, λ) with φ(T, ·) = 0 , we obtain the required result.

Viscosity solutions
Adopting the concept of viscosity solutions, it can be shown that the supremum of the objective with respect to φ is achieved by the viscosity solution of the HJB equation (10). Due to presence of the Dirac delta functions in (10), we shall introduce a suitable definition of the viscosity solution that allows to have jump discontinuities in time.
Definition 3.7. Denote by set(τ ) the set of entries of vector τ and by K the cardinality of set(τ ). Let where t k−1 < t k and t k ∈ set(τ ) for all k = 1, . . . , K .
Remark 3.9 (Comparison principle). The comparison principle still holds for viscosity solutions of (10). Let u and v be a viscosity subsolution and a viscosity supersolution of the equation (10), respectively. At the terminal time T , u(T, ·) ≤ v(T, ·). Since t K = T is in set(τ ) and u, v have the same jump , 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 Remark 3.10 (Existence and uniqueness). As a consequence of the comparison principle, there exists a unique viscosity solution of (10). The uniqueness is a direct consequence of the comparison principle. The existence can be obtained by Perron's method (see Crandall et al., 1992) under which the comparison principle is a key argument.
Now we shall prove Proposition 3.5. The proof relies on a smoothing argument used in Bouchard et al. (2017), which is based on the shaken coefficients technique of Krylov (2000). The proof is similar to Theorem 2.4 in Bouchard et al. (2017), which we sketch here for completeness.
Proof of Proposition 3.5. Denote by ϕ a viscosity solution of the equation (10) with any λ ∈ R m . From Remark 3.10, we know that such ϕ exists and is unique. The first part of the proposition is proved in two steps: Step 1. Assuming that there exists a sequence of supersolutions of (10) in BV ([0, T ], C 2 b (R d )) converging to ϕ pointwise, we can show that ϕ achieves the supremum with respect to φ in the objective of the dual (12). Let φ ∈ BV ([0, T ], C 2 b (R d )) be any solution that satisfies (13), and φ is also a (viscosity) supersolution of (10). By Remark 3.9, we have ϕ(0, The equality can be achieved in (18) by taking the supremum with respect to φ on the left-hand side of (18).
Step 2. Now, we shall construct the sequence of supersolutions required in Step 1. Let us introduce the regularising kernel r ε : where r ′ is some compactly supported function that satisfies R d r ′ (x) dx = 1 . Then we define ϕ ε = ϕ * r ε where the convolution acts only on the variable x. By applying the result of Bouchard et al. (2017) which relies critically on the fact that F * (a, b) is convex in (a, b), it can be shown that ϕ ε are supersolutions of equation (10). If we send ε to 0 , the supersolutions ϕ ε converge to the viscosity solution ϕ pointwise. The desired sequence is then constructed.

LSV Calibration
In this section, we illustrate our method by calibrating a Heston-like LSV model. This method could also be easily extended to other LSV models. We consider the LSV model with following dynamics under the risk-neutral measure: where Z t is the logarithm of the stock price at time t. The interpretations of r and q differ between financial markets. In the equity market, r is the risk-free rate and q is the dividend yield. 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 given and obtained by calibrating a pure Heston model. Note that in the literature, the widely considered LSV model has a volatility function σ(t, Z t , V t ) = L(t, Z t ) √ V t and a constant correlation η , where L(t, Z t ) is known as the leverage function. By contrast, we consider a local-stochastic volatility σ > 0 and a local-stochastic correlation η ∈ [−1, 1] whose values depend on (t, Z t , V t ). Our objective is to 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 of the variable V , the model is equivalent to a local volatility model.
Consider a probability measure P ∈ P 1 and a two-dimensional P-semi-martingale X t . The process X t has dynamics (19), i.e., X t = (Z t , V t ), if P is characterised by (α P , β P ) such that with functions σ t = σ(t, Z t , V t ) and η t = η(t, Z t , V t ). Recall that the parameters (κ, θ, ξ) are assumed to be given. Also, r t and q t are known and V t is a state variable. Hence, the only unknown variables in (20) are σ t and η t . As we will see below, σ t will be the only free variable in the calibration. Given m European options with prices c ∈ R m + , maturities τ = (τ 1 , . . . , τ m ) ∈ (0, T ] m and discounted payoffs G = (G 1 , . . . , G m ) where G i : R 2 → R + (e.g., G i (x) = e − τ i 0 r(s)ds (e x1 − K) + if the i -th option is a European call with strike K and maturity τ i , where x 1 stands for the first element of x). If X t has an initial distribution µ 0 = δ (Z0,V0) and is exactly calibrated to these European options, then P ∈ P(µ 0 , τ, c, G). One way to build a calibrated LSV model is to solve where F is a suitable convex cost function that forces (α P , β P ) to take the form of (20).
One possible way to choose the cost function F is based on the idea of minimising the difference between each element of β P and a reference value while keeping β P in S 2 + . However, it is often impossible to find an explicit formula to approximate F * . Thus numerical optimisation is needed, which makes the method computationally expensive. To overcome this issue, we choose the correlation whereη is a constant correlation obtained (along with κ, θ, ξ ) by calibrating a pure Heston model. In this case, β P t is positive semidefinite if and only if σ 2 t ≥η 2 V t for t ≤ T .
The parameter p is a constant greater than 1, and a, b, c are constants determined to minimise the function at x =x with min H = 0 .
Givenx and s satisfyingx > s, the function H is convex in x and minimised atx . It is finite only when x > s. In the numerical examples (see Section 5.2 and 5.3 below), the parameter p is set to 4 . A plot of H is given in Figure 1. Then, we define the cost function as follows.
Definition 4.3. The cost function F : where the convex set Γ is defined as By applying Proposition 3.5, the dual formulation of (21) with the cost function (23) is as follows: where φ λ is the viscosity solution to the HJB equation with a terminal condition φ λ (T, ·) = 0 .
The optimal λ can be found through a standard optimisation algorithm (see Section 5.1 below). The convergence of the algorithm can be improved by providing the gradient of the objective. Letβ λ denote the optimal β 11 that solves the supremum in (25), which also implicitly depends on λ. In fact, solving the supremum in (25) is equivalent to solving the following equation for σ 2 : for which a closed-form solution is available. We also denote by P λ ∈ P 1 a probability measure charac- . The gradient of J(λ) with respect to λ i can be formulated as: In addition, with the terminal condition φ ′ (τ i , ·) = G i .
Remark 4.6. Note that E P λ G i (X τi ) is the price of the i -th European option calculated by X t under P λ , which we refer to as the model price, and c i is the market price. Instead of solving (30) once for each option, we can perform a Monte Carlo simulation to efficiently calculate the model prices for all options. However, for the sake of accuracy, we still choose to solve (30) in the numerical examples below. Moreover, as the gradient is decreasing to zero while the solution is moving towards the optimal solution, the optimisation process can be interpreted as matching the model X t to market prices.

Numerical method
In this section, we present a numerical method for solving the dual formulation. To shorten notations, we will simply write φ for φ λ from now on. Starting with an initial λ = λ 0 (e.g., setting it to a null vector), we solve the HJB equation (25) to calculate φ(0, Z 0 , V 0 ) and hence calculate J(λ 0 ). Then, J is maximised over λ ∈ R m through an optimisation algorithm. In particular, we employed the L-BFGS algorithm (Liu and Nocedal, 1989) and obtained good convergence. The optimisation process can be accelerated by providing the gradient ∇J(λ) which can be numerically computed by (27). We measure for k = N T − 1, . . . , 0 do 4 if t k+1 is equal to the maturity of any calibrating options. then Approximate σ 2 t k by solving (26) Solve the HJB equation (25)  Calculate the gradient ∇J(λ) by (27) 12 Update λ by the L-BFGS algorithm 13 while ∇J(λ) ∞ > ǫ 1 the optimality by the maximum absolute value on the gradient. In other words, by setting a threshold ǫ 1 , the algorithm terminates when the following stopping criterion is reached: For solving the HJB equation (25), we use an alternating direction implicit (ADI) method together with the central finite difference scheme. In the numerical examples below, we employ the Douglas scheme from In 't Hout and Foulon (2010). Given a λ, we solve the HJB equation backward. Consider a discretisation {t k } of the time interval [0, T ] such that 0 = t 0 < t 1 < · · · < t NT = T, N T ∈ N. Without loss of generality, we assume that set(τ ) ⊂ {t k } . At each time step t k , we approximate σ 2 t k by solving (26) with φ = φ t k+1 for which an analytical solution can be found. At t = t k , with the approximated σ 2 t k , the HJB equation (25) is solved by the ADI finite difference method. Note that this approximation scheme of σ is similar to the one used in Ren et al. (2007) for approximating the leverage function.
Letτ i be an element in set(τ ) such that ∪ K i=1 {τ i } = set(τ ) and 0 =:τ 0 <τ 1 < . . . <τ K = T (see Definition 3.7 for the definitions of set(τ ) and K ). Denote by D the spatial computational domain and by ∂D the boundary of D . When numerically solving the HJB equation (25), we impose the following boundary conditions for the spatial dimensions: In addition, we set a sufficiently large D to reduce the impact of the boundary conditions.
To handle the jump discontinuities caused by the presence of the Dirac delta terms, we can solve the HJB equation interval-wise in the intervals separated by the maturities, and the jump discontinuity can be incorporated into the terminal condition of the HJB equation in each interval. More precisely, if t k+1 is equal to the maturity of any calibrating options, we incorporate the jump discontinuity by adding The numerical method is summarised in Algorithm 1. Due to the non-linearity of the HJB equation, when the time step sizes are too large, it might not be accurate to simply approximate σ 2 t k by solving (26) with φ = φ t k+1 once per time step. Therefore, we slightly modify the algorithm by including an iterative step to improve the accuracy of the approximation of σ 2 t k . In the literature, this iterative step is known as policy iteration, see e.g., Ma and Forsyth (2017). Specifically, at each time step t k , we first approximate σ 2 t k by solving (26) with φ = φ t k+1 . Next, we obtain φ t k by solving the HJB equation (25) with σ 2 t k , and then approximate σ 2 t k again by solving (26) with φ = φ t k . This process is repeated until φ t k converges. For completeness, the modified numerical method with policy iteration is summarised in Algorithm 2. For the sake of accuracy, we use Algorithm 2 in both Section 5.2 and Section 5.3.
In our experiments, we notice that the algorithm can provide satisfactory calibration results even with coarse grids. However, it is crucial to ensure that the grids are fine enough, because we do not want to calibrate the wrong model prices to the calibrating option prices. In fact, we observe that the algorithm converges faster with finder grids, because the numerical approximations of the gradients are more accurate with finer grids.

Numerical example: simulated data
In this section, we provide two numerical examples with simulated data to demonstrate the calibration method. In both examples, the risk-free rate is set to a constant r = 0.05 and the dividend yield is set to q = 0 . Let Z 0 = ln 100 and V 0 = 0.04 for both models. We consider a uniform mesh over the spatial

Example 1
In the first example, we use parameters (κ, θ, ξ,η) = (0.5, 0.04, 0.16, −0.4) for both the OT-LSV model and the Heston generating model. This example represents a trivial case, since if we use the same set of parameters for both models, the optimal solution of the dual formulation is a null vector λ = 0 ∈ R m , and hence V = 0 . In this case, under the optimal measure of Problem 1, σ 2 (t, Z, V ) = V and η(t, Z, V ) =η .
Setting a threshold ǫ 1 = 10 −6 , we obtain the expected results. The plot of σ 2 (t, Z, V ) is provided in Figure 2.

Example 2
In the second example, we give different parameters to the OT-LSV model and the Heston generating model (see Table 1). As noted in Remark 4.1, the OT-LSV model reduces to a LV model if σ 2 (t, Z, V ) is independent of V . Also, it is well known that an LV model can be calibrated to any arbitrage-free option prices. In this example, the Heston generating model has characteristics that are outside of Γ in the cost function F , so the Heston generating model would lead to an infinite cost. However, since the generated option prices are arbitrage free, a finite cost is still achievable by the OT-LSV model and the problem is admissible, i.e., P(µ 0 , τ, c, G) = ∅ and V < +∞.  for a subset of options are given in Table 2. If σ 2 is in the form of σ 2 (t, Z, V ) = L 2 (t, Z)V for some

Numerical example: FX market data
In this example, we calibrate the OT-LSV model to the FX options data provided in Tian et al. (2015).
The options data and the domestic and foreign yields are listed in Table 4 and Table 5. The parameters (κ, θ, ξ,η) are shown in Table 3, which are obtained by (roughly) calibrating a standard Heston model to the market option prices. In this case, 2κθ/ξ 2 = 0.169 ≪ 1 and the Feller condition is strongly violated.
Setting a threshold of ǫ 1 = 6 × 10 −6 , we obtain an exact calibration. The maximum difference between the model implied volatility and the market implied volatility is less than 1 basis point. Figure  6 shows the implied volatility of the short-maturity options (1 month and 3 months) for the market data, the uncalibrated LSV model and the OT-calibrated LSV model. Figure 7 shows the implied volatility of the long-maturity options (2 years and 5 years).

IV for T = 3m
Calibrated OT-LSV Uncalibrated OT-LSV Market options Figure 6: The implied volatility (IV) skews generated by both the uncalibrated and the calibrated OT-LSV model for 1 month and 3 months maturities in the FX market data example.

IV for T = 5Y
Calibrated OT-LSV Uncalibrated OT-LSV Market options Figure 7: The implied volatility (IV) skews generated by both the uncalibrated and the calibrated OT-LSV model for 2 years and 5 years maturities in the FX market data example.

A Appendix
If we restrict the domain of its convex conjugate Φ * : Proof. Let us identify the cases where Φ * < +∞. For any (ρ,Ā,B) ∈ C b (Λ, X ) * , using the definition of convex conjugate, we have To show that one can restrict to ρ ∈ M + (Λ, X ) if Φ * < +∞, we assume that there exists a measurable set E ⊂ Λ such that ρ(E) < 0 . By the fact that C b is dense in L 1 , there exists a sequence of nonnegative functions ζ n ∈ C b (Λ) that converges to 1 E ∈ L 1 (dρ t dt). Let us construct a sequence (r n , a n , b n ) = (−kζ n , O d×1 , O d×d ) ∈ C b (Λ, X ) where k is an arbitrary positive constant and O m×n denotes a null matrix of size m × n. It is clear that the constraint r + F * (a, b) ≤ 0 is satisfied at (r, a, b) = (r n , a n , b n ) as F * (O d×1 , O d×d ) ≤ 0 . Then, by the dominated convergence theorem, we have Φ * (ρ,Ā,B) ≥ lim n→+∞ Λ r n dρ + a n · dĀ + b n : dB If we send k to infinity, the function Φ * becomes unbounded.
To show that it is necessary to have (Ā,B) ≪ ρ if Φ * < +∞, we assume that there exists a measurable set E such that (Ā,B)(E) = 0 but ρ(E) = 0 . Again, by the fact that C b is dense in L 1 , there exists a sequence of functions ζ n ∈ C b (Λ) such that ζ n take values between 0 and 1 and the sequence converges to 1 E ∈ L 1 (dρ t dt). Such sequence can be found by taking convolution of 1 E with a standard regularising kernel. Let us construct a sequence (r n , a n , b n ) = (−F * (k 1 I d×1 , k 2 I d×d )ζ n , k 1 ζ n I d×1 , k 2 ζ n I d×d ) ∈ C b (Λ, X ) where k 1 , k 2 are arbitrary constants and I m×n denotes an all-ones matrix of size m × n. By the convexity of F * and the fact that F * (O d×1 , O d×d ) ≤ 0 , it is clear that the constraint r +F * (a, b) ≤ 0 is satisfied at (r, a, b) = (r n , a n , b n ). Then, by the dominated convergence theorem, we have Φ * (ρ,Ā,B) ≥ lim n→+∞ Λ r n dρ + a n · dĀ + b n : dB The function Φ * goes to infinity if we send k 1 , k 2 to +∞ or −∞, depending on the sign of i (Ā(E)) i and i,j (B(E)) ij . Now, since the integrand of the integral in Φ * is linear in (r, a, b), if Φ * is finite, the supremum must occur at the boundary. Thus, assuming that ρ ∈ M + (Λ, X ) and (Ā,B) ≪ ρ, we have The last equality holds since the convex and lower semi-continuous function F coincides with its biconjugate F * * according to the Fenchel-Moreau theorem (see e.g., Brezis, 2011, Theorem 1.11).
The first equality above is justified by the dominated convergence theorem. The second equality above holds because if Φ * is finite, then the supremum must occur at the boundary.