# Local Volatility Calibration by Optimal Transport

###### Abstract

The calibration of volatility models from observable option prices is a fundamental problem in quantitative finance. The most common approach among industry practitioners is based on the celebrated Dupire’s formula Dupire94pricingwith (), which requires the knowledge of vanilla option prices for a continuum of strikes and maturities that can only be obtained via some form of price interpolation. In this paper, we propose a new local volatility calibration technique using the theory of optimal transport. We formulate a time continuous martingale optimal transport problem, which seeks a martingale diffusion process that matches the known densities of an asset price at two different dates, while minimizing a chosen cost function. Inspired by the seminal work of Benamou and Brenier benamou2000computational (), we formulate the problem as a convex optimization problem, derive its dual formulation, and solve it numerically via an augmented Lagrangian method and the alternative direction method of multipliers (ADMM) algorithm. The solution effectively reconstructs the dynamic of the asset price between the two dates by recovering the optimal local volatility function, without requiring any time interpolation of the option prices.

## 1 Introduction

A fundamental assumption of the classical Black-Scholes option pricing framework is that the underlying risky asset has a constant volatility. However, this assumption can be easily dispelled by the option prices observed in the market, where the implied volatility surfaces are known to exhibit “skews” or “smiles”. Over the years, many sophisticated volatility models have been introduced to explain this phenomenon. One popular class of model is the local volatility models. In a local volatility model, the volatility function is a function of time as well as the asset price . The calibration of the local volatility function involves determining from available option prices.

One of the most prominent approaches for calibrating local volatility is introduced by the path-breaking work of Dupire Dupire94pricingwith (), which provides a method to recover the local volatility function if the prices of European call options are known for a continuum of maturities and strikes . In particular, the famous Dupire’s formula is given by

(1) |

However, in practice, option prices are only available at discrete strikes and maturities, hence interpolation is required in both variables to utilize this formula, leading to many inaccuracies. Furthermore, the numerical evaluation of the second derivative in the denominator can potentially cause instabilities in the volatility surface as well as singularities. Despite these drawbacks, Dupire’s formula and its variants are still used prevalently in the financial industry today.

In this paper, we introduce a new technique for the calibration of local volatility functions that adopts a variational approach inspired by optimal transport. The optimal transport problem was first proposed by Monge monge1781memoire () in 1781 in the context of civil engineering. The basic problem is to transfer material from one site to another while minimizing transportation cost. In the 1940’s, Kantorovich kantorovich2006problem () provided a modern treatment of the problem based on linear programming techniques, leading to the so-called Monge-Kantorovich problem. Since then, the theory of optimal transport has attracted considerable attention with applications in many areas such as fluid dynamics, meteorology and econometrics (see, e.g., evans1997partial () and villani2008optimal ()). Recently, there have been a few studies extending optimal transport to stochastic settings with applications in financial mathematics. For instance, Tan and Touzi tan2013optimal () studied an extension of the Monge-Kantorovich problem for semimartingales, while Dolinsky and Soner dolinsky2014martingale () applied martingale optimal transport to the problem of robust hedging.

In our approach, we begin by recovering the density of the underlying asset at times and from the prices of European options expiring at and . Then, instead of interpolating between different maturities, we seek a martingale diffusion process which transports the density from to , while minimizing a particular cost function. This is similar to the classical optimal transport problem, with the additional constraint that the diffusion process must be a martingale driven by a local volatility function. In the case where the cost function is convex, we find that the problem can be reformulated as a convex optimization problem under linear constraints. Theoretically, the stochastic control problem can be solved by a non-linear PDE closely connected with the ones studied in Bouchard et al. bouchard2015hedging (); bouchard2016almost () and Loeper loeper2013option () in the context of option pricing with market impact. For this paper, we approach the problem via the augmented Lagrangian method and the alternative direction method of multipliers (ADMM) algorithm, which was also used in Benamou and Brenier benamou2000computational () for classical optimal transport problems.

The paper is organized as follows. In Section 2, we introduce the classical optimal transport problem as formulated by Benamou and Brenier benamou2000computational (). In Section 3, we introduce the martingale optimal transport problem and its augmented Lagrangian. The numerical method is detailed in Section 4 and numerical results are given in Section 5.

## 2 Optimal Transport

In this section, we briefly outline the optimal transport problem as formulated by Benamou and Brenier benamou2000computational (). Given density functions with equal total mass . We say that a map is an admissible transport plan if it satisfies

(2) |

for all bounded subset . Let denote the collection of all admissible maps. Given a cost function , which represents the transportation cost of moving one unit of mass from to , the optimal transport problem is to find an optimal map that minimizes the total cost

(3) |

In particular, when where denotes the Euclidean norm, this problem is known as the Monge-Kantorovich problem (MKP).

The MKP is reformulated in benamou2000computational () in a fluid mechanic framework. In the time interval , consider all possible smooth, time-dependent, densities and velocity fields , that satisfy the continuity equation

(4) |

and the initial and final conditions

(5) |

In benamou2000computational (), it is proven that the MKP is equivalent to finding an optimal pair that minimizes

(6) |

subject to the constraints (4) and (5). This problem is then solved numerically in benamou2000computational () via an augmented Lagrangian approach. The specific numerical algorithm used is known as the alternative direction method of multipliers (ADMM), which has applications in statistical learning and distributed optimization.

## 3 Definition of the martingale problem

Let be a probability space, where is the risk-neutral measure. Suppose the dynamic of an asset price on is given by the local volatility model

(7) |

where is a local volatility function and is a one-dimensional Brownian motion. For the sake of simplicity, suppose the interest and dividend rates are zero. Denote by the density function of . It is well known that follows the Fokker-Planck equation

(8) |

Suppose that the initial and the final densities are given by and , which are recovered from European option prices via the Breeden-Litzenberger breeden1978prices () formula. Let be a convex and differentiable cost function. We are interested in minimizing the quantity

where is the support of . Unlike the classical optimal transport problem, the existence of a solution here requires an additional condition: for any convex function , there exists a martingale transport plan if and only if and satisfy:

(9) |

This is known as Strassen’s Theorem strassen1965existence ().

###### Remark 1

The formulation here is actually quite general and it can be easily adapted to a large family of models. For example, the case of a geometric Brownian motion with local volatility can be recovered by substituting everywhere, including in the Fokker-Planck equation. The cost function would then also be dependent on . The later arguments involving convex conjugates still hold since remains a convex function of .

Since is not convex in (which is crucial for our method), the substitution is applied. So we obtain the following martingale optimal transport problem:

(10) |

subject to the constraints:

(11) | |||

(12) |

Using the convexity of , the term can be easily verified to be convex in . Also note that we have the natural restrictions of and . Note that in the current version of the algorithm, for the sake of simplicity, only the condition will be encoded in the convex conjugate formulation. On the other hand, will be numerically enforced during the iterations.

Next, introduce a time-space dependent Lagrange multiplier for the constraints (11) and (12) . Hence the associated Lagrangian is

(13) |

Integrating (13) by parts and letting vanish on the boundaries of , we can rewrite the Lagrangian as follows:

(14) |

Thus the martingale optimal transport problem can be reformulated as the following saddle point problem:

(15) |

The optimality conditions can be obtained by setting the functional derivatives of with respect to , and to zero, resulting in the following equations:

(16) | ||||

(17) | ||||

(18) |

Note that (17) and (18) imply the following condition,

where is the convex conjugate of (see (19) and Proposition 1).

### Augmented Lagrangian Approach

Similar to benamou2000computational (), we solve the martingale optimal transport problem using the augmented Lagrangian approach. Let us begin by briefly recalling the well-known definition and properties of the convex conjugate. For more details, the readers are referred to Section 12 of Rockafellar rockafellar2015convex ().

Fix , let be a proper convex and lower semi-continuous function. Then the convex conjugate of is the function defined by

(19) |

The convex conjugate is also often known as the Legendre-Fenchel transform.

###### Proposition 1

We have the following properties:

(i) is a proper convex and lower semi-continuous function with ;

(ii) if is differentiable, then .

Returning to the problem at hand, recall that is convex in . By adopting the convention of whenever , it can be expressed in terms of the convex conjugate, as shown in the following proposition.

###### Proposition 2

Denote by the convex conjugate of .

(i) Let , the convex conjugate of is given by:

(20) |

(ii) For , We have the following equality,

(21) |

###### Proof

(i) By definition, the convex conjugate of is given by

(22) |

Differentiating the right hand side of (22) with respect to and setting it to zero, we obtain . Applying Proposition 1 (ii), we have

(23) |

Substituting (23) into (22) gives us

(24) |

If , the supremum is achieved by , otherwise, becomes unbounded as increases. This establishes part (i).

(ii) The required equality follows immediately from part (i) and the fact that

∎

Now we are in a position to present the augmented Lagrangian. First, let us introduce the following notations:

(25) | |||

(26) | |||

(27) | |||

(28) | |||

(29) |

By using the above notations, we can express the equality from Proposition 2 (ii) in the following way,

(30) |

Since the restriction is checked point-wise for every , we can exchange the supremum with the integrals in the following equality

(31) |

Therefore, the saddle point problem specified by (3) and (15) can be rewritten as

(32) |

Note that in the new saddle point problem (32), is the Lagrange multiplier of the new constraint . In order to turn this into a convex problem, we define the augmented Lagrangian as follows:

(33) |

where is a penalization parameter. The saddle point problem then becomes

(34) |

which has the same solution as (15).

## 4 Numerical Method

In this section, we describe in detail the alternative direction method of multipliers (ADMM) algorithm for solving the saddle point problem given by (33) and (34). In each iteration, using as a starting point, the ADMM algorithm performs the following three steps:

Step A: | (35) | ||||

Step B: | (36) | ||||

Step C: | (37) |

Step A:

To find the function that minimizes , we set the functional derivative of with respect to to zero:

(38) |

By integrating by parts, we arrive at the following variational equation

(39) |

with Neumann boundary conditions in time :

(40) | ||||

(41) |

For the boundary conditions in space, let . By setting the functional derivative of (3) with respect to to zero, we have the following equality,

(42) |

In benamou2000computational (), periodic boundary conditions were used in the spatial dimension and a perturbed equation was used to yield a unique solution. Since periodic boundary conditions are inappropriate for martingale diffusions and we are dealing with a bi-Laplacian term in space, we impose the following additional boundary conditions in order to enforce a unique solution:

(43) |

Now, the 4th order linear PDE (39) can be numerically solved by the finite difference method or the finite element method.

Step B:

Since is not differentiable, we cannot differentiate with respect to . Nevertheless, we can simply obtain by solving the minimization problem

(44) |

This is equivalent to solving

(45) |

Now, let us define

(46) |

then we can find by solving

(47) |

point-wise in space and time. This is a simple one-dimensional projection problem. If satisfies the constraint , then it is also the minimum. Otherwise, the minimum must occur on the boundary . In this case we substitute the condition into (47) to obtain

(48) |

which must be solved point-wise.
The minimum of (48) can be found using standard root finding methods such as Newton’s method. In some simple cases it is even possible to compute the solution analytically.

Step C:

Begin by computing the gradient by differentiating the augmented Lagrangian respect to . Then, simply update by moving it point-wise along the gradient as follows,

(49) |

It is important here to ensure that stays non-negative. A simple way to achieve this is to check the sign of after has been updated, and replace it by .

Stopping criteria:

Since (16) is always satisfied, combining (17) and (18) gives the following equation,

(50) |

Note that (50) holds even if the optimum occurs at . We use (50) to check the optimality. Define the residual:

(51) |

The residual is weighted by the density to alleviate any potential issues when dividing by small values of . The algorithm iterates until a particular level of tolerance is reached, in other words, .

## 5 Numerical Results

The algorithm was implemented and tested on the following simple example. Consider the computational domain and the time interval . We set the initial and final distributions to be and respectively, where denotes the normal distribution. The following cost function was chosen:

(52) |

Then we discretized the space-time domain as a lattice. The penalization parameter is set to and the tolerance is set to . In general, a higher penalization parameter may lead to faster convergence in earlier iterations, but it may produce large errors and hinder the algorithm in later iterations. In our test, satisfactory results were reached after 3000 iterations. These results are shown in Figures 1 and 2. The unstable tails in Figure 2 are numerical errors caused by small values of when computing . Otherwise, remains constant at approximately 0.0075. This matches the analytical solution of the problem, where the optimal value of is the constant .

## 6 Summary

This paper focuses on a new approach for the calibration of local volatility models. Given the distributions of the asset price at two fixed dates, the technique of optimal transport is applied to interpolate the distributions and recover the local volatility function, while maintaining the martingale property of the underlying process. Inspired by benamou2000computational (), the problem is first converted into a saddle point problem, and then solved numerically by an augmented Lagrangian approach and the alternative direction method of multipliers (ADMM) algorithm. The algorithm performs well on a simple case in which the numerical solution matches its analytical counterpart. The main drawback of this method is due to the slow convergence rate of the ADMM algorithm. Further research is required to improve the efficiency of the algorithm and apply it to more complex cases.

## Bibliography

- (1) Benamou, J.D., Brenier, Y.: A computational fluid mechanics solution to the monge-kantorovich mass transfer problem. Numerische Mathematik 84(3), 375–393 (2000)
- (2) Bouchard, B., Loeper, G., Zou, Y.: Hedging of covered options with linear market impact and gamma constraint. arXiv preprint arXiv:1512.07087 (2015)
- (3) Bouchard, B., Loeper, G., Zou, Y.: Almost-sure hedging with permanent price impact. Finance and Stochastics 20(3), 741–771 (2016)
- (4) Breeden, D.T., Litzenberger, R.H.: Prices of state-contingent claims implicit in option prices. Journal of business pp. 621–651 (1978)
- (5) Dolinsky, Y., Soner, H.M.: Martingale optimal transport and robust hedging in continuous time. Probability Theory and Related Fields 160(1-2), 391–427 (2014)
- (6) Dupire, B., (see Black, T.B.M., Options, G.: Pricing with a smile. Risk Magazine pp. 18–20 (1994)
- (7) Evans, L.C.: Partial differential equations and monge-kantorovich mass transfer. Current developments in mathematics 1997(1), 65–126 (1997)
- (8) Kantorovich, L.V.: On a problem of monge. Journal of Mathematical Sciences 133(4), 1383–1383 (2006)
- (9) LOEPER, G.: Option pricing with market impact and non-linear black and scholes pdesâs. arXiv preprint arXiv:1301.6252 (2013)
- (10) Monge, G.: Mémoire sur la théorie des déblais et des remblais. Histoire de l’Académie Royale des Sciences de Paris (1781)
- (11) Rockafellar, R.T.: Convex analysis. Princeton university press (2015)
- (12) Strassen, V.: The existence of probability measures with given marginals. The Annals of Mathematical Statistics pp. 423–439 (1965)
- (13) Tan, X., Touzi, N., et al.: Optimal transportation under controlled stochastic dynamics. The annals of probability 41(5), 3201–3240 (2013)
- (14) Villani, C.: Optimal transport: old and new, vol. 338. Springer Science & Business Media (2008)