A Generative Adversarial Network Approach to Calibration of Local Stochastic Volatility Models

Abstract

We propose a fully data-driven approach to calibrate local stochastic volatility (LSV) models, circumventing in particular the ad hoc interpolation of the volatility surface. To achieve this, we parametrize the leverage function by a family of feed-forward neural networks and learn their parameters directly from the available market option prices. This should be seen in the context of neural SDEs and (causal) generative adversarial networks: we generate volatility surfaces by specific neural SDEs, whose quality is assessed by quantifying, possibly in an adversarial manner, distances to market prices. The minimization of the calibration functional relies strongly on a variance reduction technique based on hedging and deep hedging, which is interesting in its own right: it allows the calculation of model prices and model implied volatilities in an accurate way using only small sets of sample paths. For numerical illustration we implement a SABR-type LSV model and conduct a thorough statistical performance analysis on many samples of implied volatility smiles, showing the accuracy and stability of the method.

1. Introduction

Each day a crucial task is performed in financial institutions all over the world: the calibration of stochastic models to current market or historical data. So far the model choice was not only driven by the capacity of capturing empirically observed market features well, but also by the computational tractability of the calibration process. This is now undergoing a big change since machine-learning technologies offer new perspectives on model calibration.

Calibration is the choice of one model from a pool of models, given current market and historical data. Depending on the nature of data this is considered to be an inverse problem or a problem of statistical inference. We consider here current market data, in particular volatility surfaces, therefore we rather emphasize the inverse problem point of view. We however stress that it is the ultimate goal of calibration to include both data sources simultaneously. In this respect machine learning might help considerably.

We can distinguish three kinds of machine learning-inspired approaches for calibration to current market prices: First, having solved the inverse problem already several times, one can learn from this experience (i.e., training data) the calibration map from market data to model parameters directly. Let us here mention one of the pioneering papers by Hernandez (2017) that applied neural networks to learn this calibration map in the context of interest rate models. This was taken up in Cuchiero et al. (2018) for calibrating more complex mixture models. Second, one can learn the map from model parameters to model prices (compare e.g., Liu et al. (2019a, 2019b)) and then invert this map possibly with machine learning technology. In the context of rough volatility modeling, see Gatheral et al. (2018), such approaches turned out to be very successful: we refer here to Bayer et al. (2019) and the references therein. Third, the calibration problem is considered to be the search for a model which generates given market prices and where additionally technology from generative adversarial networks, first introduced by Goodfellow et al. (2014), can be used. This means parameterizing the model pool in a way which is accessible for machine learning techniques and interpreting the inverse problem as a training task of a generative network, whose quality is assessed by an adversary. We pursue this approach in the present article and use as generative models so-called neural stochastic differential equations (SDE), which just means to parameterize the drift and volatility of an Itô-SDE by neural networks.

1.1. Local Stochastic Volatility Models as Neural SDEs

We focus here on calibration of local stochastic volatility (LSV) models, which are in view of existence and uniqueness still an intricate model class. LSV models, going back to Jex (1999); Lipton (2002); Ren et al. (2007), combine classical stochastic volatility with local volatility to achieve both a good fit to time series data and in principle a perfect calibration to the implied volatility smiles and skews. In these models, the discounted price process ${\left(S_t\right)}_{t\ge 0}$ of an asset satisfies

$$\begin{array}{c}\hfill d{S}_t=S_{t}L(t,S_t){\alpha }_{t}d{W}_t,\end{array}$$

(1)

where ${\left({\alpha }_t\right)}_{t\ge 0}$ is some stochastic process taking values in $\mathbb{R}$, and (a sufficiently regular function) $L(t,s)$ the so-called leverage function depending on time and the current value of the asset and W a one-dimensional Brownian motion. Note that the stochastic volatility process $\alpha$ can be very general and could for instance be chosen as rough volatility model. By slight abuse of terminology we call $\alpha$ stochastic volatility even though it is strictly speaking not the volatility of the log price of S.

For notational simplicity we consider here the one-dimensional case, but the setup easily translates to a multivariate situation with several assets and a matrix valued analog of $\alpha$ as well as a matrix valued leverage function.

The leverage function L is the crucial part in this model. It allows in principle to perfectly calibrate the implied volatility surface seen on the market. To achieve this goal L must satisfy

$$L^2(t,s)=\frac{{\sigma }_{\mathrm{Dup}}^2(t,s)}{\mathbb{E}\left[{\alpha }_t^2\right|S_t=s]},$$

(2)

where ${\sigma }_{\mathrm{Dup}}$ denotes Dupire’s local volatility function (see Dupire (1994, 1996)). For the derivation of (2), we refer to Guyon and Henry-Labordѐre (2013). Please note that (2) is an implicit equation for L as it is needed for the computation of $\mathbb{E}\left[{\alpha }_t^2\right|S_t=s]$. This in turn means that the SDE for the price process ${\left(S_t\right)}_{t\ge 0}$ is actually a McKean–Vlasov SDE, since the law of $(S_t,{\alpha }_t)$ enters in the characteristics of the equation. Existence and uniqueness results for this equation are not at all obvious, since the coefficients do not satisfy any kind of standard conditions like for instance Lipschitz continuity in the Wasserstein space. Existence of a short-time solution of the associated nonlinear Fokker-Planck equation for the density of ${\left(S_t\right)}_{t\ge 0}$ was shown in Abergel and Tachet (2010) under certain regularity assumptions on the initial distribution. As stated in Guyon and Henry-Labordѐre (2012) a very challenging and still open problem is to derive the set of stochastic volatility parameters for which LSV models exist uniquely for a given market implied volatility surface. We refer to Jourdain and Zhou (2016) and Lacker et. al. (2019), where recent progress in solving this problem has been made.

Despite these intriguing existence issues, LSV models have attracted—due to their appealing feature of a potentially perfect smile calibration and their econometric properties—a lot of attention from the calibration and implementation point of view. We refer to Cozma et al. (2017); Guyon and Henry-Labordѐre (2012); Guyon and Henry-Labordѐre (2013); for Monte Carlo (MC) methods (see also Guyon (2014, 2016) for the multivariate case), to Ren et al. (2007); Tian et al. (2015) for PDE methods based on nonlinear Fokker-Planck equations and to Saporito et al. (2017) for inverse problem techniques. Within these approaches the particle approximation method for the McKean–Vlasov SDE proposed in Guyon and Henry-Labordѐre (2012); Guyon and Henry-Labordѐre (2013) works impressively well, as very few paths must be used to achieve very accurate calibration results.

In the current paper we propose an alternative, fully data-driven approach circumventing in particular the interpolation of the volatility surface, being necessary in several other approaches in order to compute Dupire’s local volatility. This means that we only take the available discrete data into account and do not generate a continuous surface interpolating between the given market option prices. Indeed, we just learn or train the leverage function L to generate the available market option prices accurately. Although in principle the method allows for calibration to any traded options, we work here with vanilla derivatives.

Setting $T_0=0$ and denoting by $T_1<T_2\cdots <T_n$ the maturities of the available options, we parametrize the leverage function $L(t,s)$ via a family of neural networks $F^i:\mathbb{R}\to \mathbb{R}$ with weights ${\theta }_i\in {\Theta }_i$, i.e.

$$L(t,s,\theta )=1+F^i(s,{\theta }_i),t\in [T_{i-1},T_i),i\in \{1,\dots ,n\}.$$

We here consider for simplicity only the univariate case. The multivariate situation just means that 1 is replaced by the identity matrix and the neural networks $F^i(·,{\theta }_i)$ are maps from ${\mathbb{R}}^d\to {\mathbb{R}}^{d\times d}$.

This then leads to the generative model class of neural SDEs (see Gierjatowicz et al. (2020) for related work), which in the case of time-inhomogeneous Itô-SDEs, just means to parametrize the drift $\mu (·,·,\theta )$ and volatility $\sigma (·,·,\theta )$ by neural networks with parameters $\theta$, i.e.,

$$\begin{array}{c}\hfill d{X}_t\left(\theta \right)=\mu (X_t\left(\theta \right),t,\theta )dt+\sigma (X_t\left(\theta \right),t,\theta )d{W}_t,X_0\left(\theta \right)=x.\end{array}$$

(3)

In our case, there is no drift and the volatility (for the price) reads as

$$\begin{array}{c}\hfill \sigma (S_t\left(\theta \right),t,\theta )=S_t\left(\theta \right)\left(1+\sum _{i=1}^n{F}^i(S_t\left(\theta \right),{\theta }_i)1_{[T_{i-1},T_i)}\left(t\right)\right){\alpha }_t.\end{array}$$

(4)

Progressively for each maturity, the parameters of the neural networks are learned by optimizing the following calibration criterion

$$\begin{array}{c}\hfill \underset{\theta }{inf}\underset{\gamma }{sup}\sum _{j=1}^J{w}_j^{\gamma }{\ell }^{\gamma }({\pi }_j^{\mathrm{mod}}\left(\theta \right)-{\pi }_j^{\mathrm{mkt}}),\end{array}$$

(5)

where J is the number of considered options and where ${\pi }_j^{\mathrm{mod}}\left(\theta \right)$ and ${\pi }_j^{\mathrm{mkt}}$ stand for the respective model and market prices.

Moreover, for every fixed $\gamma$, ${\ell }^{\gamma }$ is a nonlinear non-negative convex function with ${\ell }^{\gamma }\left(0\right)=0$ and ${\ell }^{\gamma }\left(x\right)>0$ for $x\ne 0$, measuring the distance between model and market prices. The terms $w_j^{\gamma }$, for fixed $\gamma$, denote some weights, in our case of vega type (compare Cont and Ben Hamida (2004)). Using such vega type weights allows to match implied volatility data, our actual goal, very well. The parameters $\gamma$ take here the role of the adversarial part. Indeed, by considering a family of weights and loss functions parameterized by $\gamma$ enables us to take into account the uncertainty of the loss function. In this sense this constitutes a discriminative model as it modifies the distribution of the target, i.e., the different market prices, given $\theta$ and thus the model prices. This can for instance mean that the adversary chooses the weights $w^{\gamma }$ in such a way to put most mass on those options where the fit is worst or that it modifies ${\ell }^{\gamma }$ by choosing ${\pi }^{\mathrm{mkt}}$ within the bid-ask spread with the largest possible distance to the current model prices. In a concrete implementation in Section 4.3, we build a family of loss functions like this, using different market implied volatilities lying in the bid-ask spread. We can therefore also solve such kind of robust calibration problems.

The precise algorithms are outlined in Section 3 and Section 4, where we also conduct a thorough statistical performance analysis. Notice that as is somehow quite typical for financial applications, we need to guarantee a very high accuracy, whence a variance reduction technique to compute the model prices via Monte Carlo is crucial for this learning task. This relies on hedging and deep hedging, which allows the computation of accurate model prices ${\pi }^{\mathrm{mod}}\left(\theta \right)$ for training purposes with only up to $5\times {10}^4$ trajectories. Let us remark that we do not aim to compete with existing algorithms, as e.g., the particle method by Guyon and Henry-Labordѐre (2012); Guyon and Henry-Labordѐre (2013), in terms of speed but rather provide a generic data-driven algorithm that is universally applicable for all kind of options, also in multivariate situations, without resorting to Dupire type volatilities. This general applicability comes at the expense of a higher computation time compared to Guyon and Henry-Labordѐre (2012); Guyon and Henry-Labordѐre (2013). In terms of accuracy, we achieve an average calibration error of about 5 to 10 basis points, whence our method is comparable or in some situations even better than Guyon and Henry-Labordѐre (2012) (compare Section 6 and the results in Guyon and Henry-Labordѐre (2012)). Moreover, we also observe good extrapolation and generalization properties of the calibrated leverage function.

1.2. Generative Adversarial Approaches in Finance

The above introduced generative adversarial approach might seem at first sight unexpected as generative adversarial models or networks are rather applied in areas such as photorealistic image generation. From an abstract point of view, however, a generative network is nothing else than a neural network (mostly of recurrent type) $G^{\theta }$ depending on parameters $\theta$, which transports a standard input law ${\mathbb{P}}_I$ to a target output law ${\mathbb{P}}_O$. In our case, ${\mathbb{P}}_I$ corresponds to the law of the Brownian motion W and the stochastic volatility process $\alpha$. The target law ${\mathbb{P}}_O$ is given in terms of certain functionals, namely the set of market prices, and is thus not fully specified.

Denoting the push-forward of ${\mathbb{P}}_I$ under the transport map $G^{\theta }$ by $G_{*}^{\theta }{\mathbb{P}}_I$, the goal is to find parameters $\theta$ such that $G_{*}^{\theta }{\mathbb{P}}_I\approx {\mathbb{P}}_O$. For this purpose, appropriate distance functions must be used. Standard examples include entropies, integral distances, Wasserstein or Radon distances, etc. The adversarial character appears when the chosen distance is represented as supremum over certain classes of functions, which can themselves be parameterized via neural networks of a certain type. This leads to a game, often of zero-sum type, between the generator and the adversary. As mentioned above, one example, well known from industry, is calibrating a model to generate prices within a bid-ask spread: in this case there is more than one loss function, each of them representing a distance to a possible price structure, and the supremum over these loss functions is the actual distance between generated price structure and the target. In other words: the distance to the worst possible price structure within the bid-ask spread should be as small as possible (see Section 4.3).

In our case the solution measure of the neural SDE as specified in (3) and (4) corresponds to the transport $G_{*}^{\theta }{\mathbb{P}}_I$ and we measure the distance by (5), which can be rewritten as

$$\begin{array}{cc}\hfill & \underset{\theta }{inf}\underset{\gamma }{sup}\sum _{j=1}^J{w}_j^{\gamma }{\ell }^{\gamma }(\underset{\mathrm{model}\mathrm{price}}{\underbrace{{\mathbb{E}}_{G_{*}^{\theta }{\mathbb{P}}_I}\left[C_j\right]}}-\underset{\mathrm{market}\mathrm{price}}{\underbrace{{\mathbb{E}}_{{\mathbb{P}}_O}\left[C_j\right]}}),\hfill \end{array}$$

(6)

where $C_j$ are the corresponding option payoffs.

In general, one could consider distance functions $d^{\gamma }$ such that the game between generator and adversary appears as

$$\begin{array}{c}\hfill \underset{\theta }{inf}\underset{\gamma }{sup}{d}^{\gamma }(G_{*}^{\theta }{\mathbb{P}}_I,{\mathbb{P}}_O).\end{array}$$

The advantage of this point of view is two-fold:

• we have access to the unreasonable effectiveness of modeling by neural networks, due to their good generalization and regularization properties;
• the game theoretic view disentangles realistic price generation from discriminating with different loss functions, parameterized by $\gamma$. This reflects the fact that it is not necessarily clear which loss function one should use. Notice that (6) is not the usual form of generative adversarial network (GAN) problems, since the adversary distance ${\ell }^{\gamma }$ is nonlinear in ${\mathbb{P}}_I$ and ${\mathbb{P}}_O$, but we believe that it is worth taking this abstract point of view.

There is no reason these generative models, if sufficient computing power is available, should not take market price data as inputs, too. This would correspond, from the point of view of generative adversarial networks, to actually learn a map $G^{\theta ,\mathrm{market}\mathrm{prices}}$, such that for any price configuration of market prices one has instantaneously a generative model given, which produces those prices. This requires just a rich data source of typical market prices (and computing power!).

Even though it is usually not considered like that, one can also view the generative model as an engine producing a likelihood function on probability measures on path space: given historic data, ${\mathbb{P}}_O$ is then just the empirical measure of the one observed trajectory that is inserted in the likelihood function. This would allow, with precisely the same technology, a maximum likelihood approach, where one searches for those parameters of the generative network that maximize the likelihood of the historic trajectory. This then falls in the realm of generative approaches that appear in the literature under the name “market generators”. Here the goal is to precisely mimic the behavior and features of historical market trajectories. This line of research has been recently pursued in e.g., Kondratyev and Schwarz (2019); Wiese et al. (2019); Acciaio and Xu (2020); Bühler et al. (2020); Henry-Labordѐre (2019).

From a bird’s eye perspective this machine-learning approach to calibration might just look like a standard inverse problem with another parameterized family of functions. We, however, insist on one important difference, namely implicit regularizations (see e.g., Heiss et al. (2019)), which always appear in machine-learning applications and which are cumbersome to mimic in classical inverse problems.

Finally, let us comment more generally on machine-learning approaches in mathematical finance, which become more and more prolific. Concrete applications include hedging Bühler et al. (2019), portfolio selection Gao et al. (2019), stochastic portfolio theory Samo and Vervuurt (2016); Cuchiero et al. (2020), optimal stopping Becker et al. (2019), optimal transport and robust finance Eckstein and Kupper (2019), stochastic games and control problems Huré et al. (2018) as well as high-dimensional nonlinear partial differential equations (PDEs) Han et al. (2017); Huré et al. (2019). Machine learning also allows for new insights into structural properties of financial markets as investigated in Sirignano and Cont (2019). For an exhaustive overview of machine-learning applications in mathematical finance, in particular for option pricing and hedging we refer to Ruf and Wang (forthcoming).

The remainder of the article is organized as follows. Section 2 introduces the variance reduction technique based on hedge control variates, which is crucial in our optimization tasks. In Section 3 we explain our calibration method, in particular how to optimize (5). The details of the numerical implementation and the results of the statistical performance analysis are then given in Section 4 as well as Section 6. In Appendix A we state stability theorems for stochastic differential equations depending on parameters. This is applied to neural SDEs when calculating derivatives with respect to the parameters of the neural networks. In Appendix B we recall preliminaries on deep learning by giving a brief overview of universal approximation properties of artificial neural networks and briefly explaining stochastic gradient descent. Finally, Appendix C contains alternative optimization approaches to (5).

2. Variance Reduction for Pricing and Calibration Via Hedging and Deep Hedging

This section is dedicated to introducing a generic variance reduction technique for Monte Carlo pricing and calibration by using hedging portfolios as control variates. This method will be crucial in our LSV calibration presented in Section 3. For similar considerations we refer to Potters et al. (2001); Vidales et al. (2018).

Consider on a finite time horizon $T>0$, a financial market in discounted terms with r traded instruments ${\left(Z_t\right)}_{t\in [0,T]}$ following an ${\mathbb{R}}^r$-valued stochastic process on some filtered probability space $(\Omega ,{\left({\mathcal{F}}_t\right)}_{t\in [0,T]},\mathcal{F},\mathbb{Q})$. Here, $\mathbb{Q}$ is a risk neutral measure and ${\left({\mathcal{F}}_t\right)}_{t\in [0,T]}$ is supposed to be right continuous. In particular, we suppose that ${\left(Z_t\right)}_{t\in [0,T]}$ is an r-dimensional square integrable martingale with càdlàg paths.

Let C be an ${\mathcal{F}}_T$-measurable random variable describing the payoff of some European option at maturity $T>0$. Then the usual Monte Carlo estimator for the price of this option is given by

$$\pi =\frac{1}{N}\sum _{n=1}^N{C}_n,$$

(7)

where $(C_1,\dots ,C_N)$ are i.i.d with the same distribution as C and $N\in \mathbb{N}$. This estimator can easily be modified by adding a stochastic integral with respect to Z. Indeed, consider a strategy ${\left(h_t\right)}_{t\in [0,T]}\in L^2\left(Z\right)$ and some constant c. Denote the stochastic integral with respect to Z by $I={(h•Z)}_T$ and consider the following estimator

$$\begin{array}{c}\hfill \widehat{\pi }=\frac{1}{N}\sum _{n=1}^N(C_n-c{I}_n),\end{array}$$

(8)

where $(I_1,\dots ,I_N)$ are i.i.d with the same distribution as I. Then, for any ${\left(h_t\right)}_{t\in [0,T]}\in L^2\left(Z\right)$ and c, this estimator is still an unbiased estimator for the price of the option with payoff C since the expected value of the stochastic integral vanishes. If we denote by

$$H=\frac{1}{N}\sum _{n=1}^N{I}_n,$$

then the variance of $\widehat{\pi }$ is given by

$$\mathrm{Var}\left(\widehat{\pi }\right)=\mathrm{Var}\left(\pi \right)+c^2\mathrm{Var}\left(H\right)-2c\mathrm{Cov}(\pi ,H).$$

This becomes minimal by choosing

$$c=\frac{\mathrm{Cov}(\pi ,H)}{Var\left(H\right)}.$$

With this choice, we have

$$\mathrm{Var}\left(\widehat{\pi }\right)=(1-{\mathrm{Corr}}^2(\pi ,H))\mathrm{Var}\left(\pi \right).$$

In particular, in the case of a perfect pathwise hedge, where $\pi =H$ a.s., we have $\mathrm{Corr}(\pi ,H)=1$ and $\mathrm{Var}\left(\widehat{\pi }\right)=0$, since in this case

$$\mathrm{Var}\left(\pi \right)=\mathrm{Var}\left(H\right)=\mathrm{Cov}(\pi ,H).$$

Therefore, it is crucial to find a good approximate hedging portfolio such that ${\mathrm{Corr}}^2(\pi ,H)$ becomes large. This is subject of Section 2.1 and Section 2.2 below.

2.1. Black–Scholes Delta Hedge

In many cases, of local stochastic volatility models as of form (1) and options depending only on the terminal value of the price process, a Delta hedge of the Black–Scholes model works well. Indeed, let $C=g\left(S_T\right)$ and let ${\pi }_{\mathrm{BS}}^g(t,T,s,\sigma )$ be the price at time t of this claim in the Black–Scholes model. Here, s stands for the price variable and $\sigma$ for the volatility parameter in the Black–Scholes model. Moreover, we indicate the dependency on the maturity T as well. Then choosing as hedging instrument only the price S itself and as approximate hedging strategy

$$h_t={\partial }_s{\pi }_{\mathrm{BS}}^g(t,T,S_t,L(t,S_t){\alpha }_t)$$

(9)

usually already yields a considerable variance reduction. In fact, it is even sufficient to consider ${\alpha }_t$ alone to achieve satisfying results, i.e., one has

$$h_t={\partial }_s{\pi }_{\mathrm{BS}}^g(t,T,S_t,{\alpha }_t),$$

(10)

This reduces the computational costs for the evaluation of the hedging strategies even further.

2.2. Hedging Strategies as Neural Networks—Deep Hedging

Alternatively, in particular when the number of hedging instruments becomes higher, one can learn the hedging strategy by parameterizing it via neural networks. For a brief overview of neural networks and relevant notation used below, we refer to Appendix B.

Let the payoff be again a function of the terminal values of the hedging instruments, i.e., $C=g\left(Z_T\right)$. Then in Markov models it makes sense to specify the hedging strategy via a function $h:{\mathbb{R}}_{+}\times {\mathbb{R}}^r\to {\mathbb{R}}^r$

$$h_t=h(t,z),$$

which in turn will correspond to an artificial neural network $(t,z)\mapsto h(t,z,\delta )\in {\mathcal{NN}}_{r+1,r}$ with weights denoted by $\delta$ in some parameter space $\Delta$ (see Notation1 A1). Following the approach in (Bühler et al. 2019, Remark 3), an optimal hedge for the claim C with given market price ${\pi }^{\mathrm{mkt}}$ can be computed via

$$\underset{\delta \in \Delta }{inf}\mathbb{E}\left[u\left(-C+{\pi }^{\mathrm{mkt}}+{(h(·,Z_{·-},\delta )•Z_{·})}_T\right)\right]$$

for some convex loss function $u:\mathbb{R}\to {\mathbb{R}}_{+}$. Recall that ${(h•Z)}_T$ denotes the stochastic integral with respect to Z. If $u\left(x\right)=x^2$, which is often used in practice, this then corresponds to a quadratic hedging criterion.

To tackle this optimization problem, we can apply stochastic gradient descent, because we fall in the realm of problem (A3). Indeed, the stochastic objective function $Q\left(\delta \right)\left(\omega \right)$ is given by

$$Q\left(\delta \right)\left(\omega \right)=u(-C\left(\omega \right)+{\pi }^{\mathrm{mkt}}+{(h(·,Z_{·-},\delta )\left(\omega \right)•Z_{·}\left(\omega \right))}_T).$$

The optimal hedging strategy $h(·,·,{\delta }^{*})$ for an optimizer ${\delta }^{*}$ can then be used to define

$${(h(·,Z_{·-},{\delta }^{*})•Z_{·})}_T$$

which is in turn used in (8).

As always in this article we shall assume that activation functions $\varphi$ of the neural network as well as the convex loss function u are smooth, hence we can calculate derivatives with respect to $\delta$ in a straight forward way. This is important to apply stochastic gradient descent, see Appendix B.2. We shall show that the gradient of $Q\left(\delta \right)$ is given by

$${\nabla }_{\delta }Q\left(\delta \right)\left(\omega \right)=u^{\prime }(-C\left(\omega \right)+{\pi }^{\mathrm{mkt}}+{(h(·,Z_{·-},\delta )\left(\omega \right)•Z_{·}\left(\omega \right))}_T){({\nabla }_{\delta }h(·,Z_{·-},\delta )\left(\omega \right)•Z_{·}\left(\omega \right))}_T,$$

i.e., we are allowed to move the gradient inside the stochastic integral, and that approximations with simple processes, as we shall do in practice, converge to the correct quantities. To ensure this property, we shall apply the following theorem, which follows from results in Appendix A.

Theorem 1.

For $\epsilon \ge 0$, let $Z^{\epsilon }$ be a solution of a stochastic differential equation as described in Theorem A2 with drivers $Y=(Y^1,\dots ,Y^d)$, functionally Lipschitz operators $F_j^{\epsilon ,i}$, $i=1,\dots ,r$, $j=1,\dots ,d$ and a process $(J^{\epsilon ,1},\dots J^{\epsilon ,r})$, which is here for all $\epsilon \ge 0$ simply $J{1}_{\{t=0\}}\left(t\right)$ for some constant vector $J\in {\mathbb{R}}^r=J$, i.e.

$$Z_t^{\epsilon ,i}=J^i+\sum _{j=1}^d{\int }_0^t{F}_j^{\epsilon ,i}{\left(Z^{\epsilon }\right)}_{s-}d{Y}_s^j,t\ge 0.$$

Let $(\epsilon ,t,z)\mapsto f^{\epsilon }(t,z)$ be a map, such that the bounded càglàd process $f^{\epsilon }:=f^{\epsilon }(.-,Z_{.-}^0)$ converges ucp to $f^0:=f^0(.-,Z_{.-}^0)$, then

$$\underset{\epsilon \to 0}{lim}(f^{\epsilon }•Z^{\epsilon })=\left(f^0•Z^0\right)$$

holds true.

Proof. 

Consider the extended system

$$d(f^{\epsilon }•Z^{\epsilon })=\sum _{j=1}^d{f}^{\epsilon }(t-,Z_{t-}^{\epsilon })F_j^{\epsilon ,i}{\left(Z^{\epsilon }\right)}_{t-}d{Y}_t^j$$

and

$$d{Z}_t^{\epsilon ,i}=\sum _{j=1}^d{F}_j^{\epsilon ,i}{\left(Z^{\epsilon }\right)}_{t-}d{Y}_t^j,$$

where we obtain existence, uniqueness and stability for the second equation by Theorem A2, and from where we obtain ucp convergence of the integrand of the first equation: since stochastic integration is continuous with respect to the ucp topology we obtain the result. □

The following corollary implies the announced properties, namely that we can move the gradient inside the stochastic integral and that the derivatives of a discretized integral with a discretized version of Z and approximations of the hedging strategies are actually close to the derivatives of the limit object.

Corollary 1.

Let, for $\epsilon >0$, $Z^{\epsilon }$ denote a discretization of the process of hedging instruments $Z\equiv Z^0$ such that the conditions of Theorem 1 are satisfied. Denote, for $\epsilon \ge 0$, the corresponding hedging strategies by $(t,z,\delta )\mapsto h^{\epsilon }(t,z,\delta )$ given by neural networks ${\mathcal{NN}}_{r+1,r}$, whose activation functions are bounded and $C^1$, with bounded derivatives.

1. 

Then the derivative ${\nabla }_{\delta }(h(·,Z_{·-},\delta )•Z)$ in direction δ at ${\delta }_0$ satisfies

$${\nabla }_{\delta }(h(·,Z_{·-},{\delta }_0)•Z)=({\nabla }_{\delta }h(·,Z_{·-},{\delta }_0)•Z).$$

2. 

If additionally the derivative in direction δ at ${\delta }_0$ of ${\nabla }_{\delta }{h}^{\epsilon }(·,Z_{·-},{\delta }_0)$ converges ucp to ${\nabla }_{\delta }h(·,Z_{·-},{\delta }_0)$ as $\epsilon \to 0$, then the directional derivative of the discretized integral, i.e.

${\nabla }_{\delta }(h^{\epsilon }(·,Z_{·-}^{\epsilon },{\delta }_0)•Z^{\epsilon })$ or equivalently $({\nabla }_{\delta }{h}^{\epsilon }(·,Z_{·-}^{\epsilon },{\delta }_0)•Z^{\epsilon })$, converges, as the discretization mesh $\epsilon \to 0$, to

$$\underset{\epsilon \to 0}{lim}({\nabla }_{\delta }{h}^{\epsilon }(·,Z_{·-}^{\epsilon },{\delta }_0)•Z^{\epsilon })=({\nabla }_{\delta }h(·,Z_{·-},{\delta }_0)•Z).$$

Proof. 

To prove (1), we apply Theorem 1 with

$$f^{\epsilon }(·,Z_{·-})=\frac{h(·,Z_{·-},{\delta }_0+\epsilon \delta )-h(·,Z_{·-},{\delta }_0)}{\epsilon },$$

which converges ucp to $f^0={\nabla }_{\delta }h(·,Z_{·-},{\delta }_0)$. Indeed, by the neural network assumptions, we have (with the sup over some compact set)

$$\underset{\epsilon \to 0}{lim}\underset{(t,z)}{sup}∥\frac{h(t,z,{\delta }_0+\epsilon \delta )-h(t,z,{\delta }_0)}{\epsilon }-{\nabla }_{\delta }h(t,z,{\delta }_0)∥=0,$$

by equicontinuity of $\{(t,z)\mapsto {\nabla }_{\delta }h(t,z,{\delta }_0+\epsilon \delta )|\epsilon \in [0,1]\}$.

Concerning (2) we apply again Theorem 1, this time with

$$f^{\epsilon }(·,Z_{·-})={\nabla }_{\delta }{h}^{\epsilon }(·,Z_{·-},{\delta }_0),$$

which converges by assumption ucp to $f^0={\nabla }_{\delta }h(·,Z_{·-},{\delta }_0)$. □

3. Calibration of LSV Models

Consider an LSV model as of (1) defined on some filtered probability space $(\Omega ,{\left({\mathcal{F}}_t\right)}_{t\in [0,T]},\mathcal{F},\mathbb{Q})$, where $\mathbb{Q}$ is a risk neutral measure. We assume the stochastic process $\alpha$ to be fixed. This can for instance be achieved by first approximately calibrating the pure stochastic volatility model with $L\equiv 1$, so to capture only the right order of magnitude of the parameters and then fixing them.

Our main goal is to determine the leverage function L in perfect accordance with market data. We here consider only European call options, but our approach allows in principle to take all kind of other options into account.

Due to the universal approximation properties outlined in Appendix B (Theorem A3) and in spirit of neural SDEs, we choose to parameterize L via neural networks. More precisely, set $T_0=0$ and let $0<T_1\cdots <T_n=T$ denote the maturities of the available European call options to which we aim to calibrate the LSV model. We then specify the leverage function $L(t,s)$ via a family of neural networks, i.e.,

$$\begin{array}{c}\hfill L(t,s,\theta )=\left(1+\sum _{i=1}^n{F}^i(s,{\theta }_i)1_{[T_{i-1},T_i)}\left(t\right)\right),\end{array}$$

(11)

where $F^i\in {\mathcal{NN}}_{1,1}$ for $i=1,\dots ,n$ (see Notation A1). For notational simplicity we shall often omit the dependence on ${\theta }_i\in {\Theta }_i$. However, when needed we write for instance $S_t\left(\theta \right)$, where $\theta$ then stands for all parameters ${\theta }_i$ used up to time t.

For purposes of training, similarly as in Section 2.2, we shall need to calculate derivatives of the LSV process with respect to $\theta$. The following result can be understood as the chain rule applied to ${\nabla }_{\theta }S\left(\theta \right)$, which we prove here rigorously by applying the results of Appendix A.

Theorem 2.

Let $(t,s,\theta )\mapsto L(t,s,\theta )$ be of form (11) where the neural networks $(s,{\theta }_i)\mapsto F^i(s,{\theta }_i)$ are bounded and $C^1$, with bounded and Lipschitz continuous derivatives2, for all $i=1,\dots ,n$. Then the directional derivative in direction θ at $\widehat{\theta }$ satisfies the following equation

$$\begin{array}{cc}\hfill d\left({\nabla }_{\theta }{S}_t\left(\widehat{\theta }\right)\right)& =({\nabla }_{\theta }{S}_t\left(\widehat{\theta }\right)L(t,S_t\left(\widehat{\theta }\right),\widehat{\theta })+S_t\left(\widehat{\theta }\right){\partial }_{s}L(t,S_t\left(\widehat{\theta }\right),\widehat{\theta }){\nabla }_{\theta }{S}_t\left(\widehat{\theta }\right)\hfill \\ \hfill & +S_t\left(\widehat{\theta }\right){\nabla }_{\theta }L(t,S_t\left(\widehat{\theta }\right),\widehat{\theta })){\alpha }_{t}d{W}_t,\hfill \end{array}$$

(12)

with initial value 0. This can be solved by variation of constants, i.e.

$$\begin{array}{c}\hfill {\nabla }_{\theta }{S}_t\left(\widehat{\theta }\right)={\int }_0^t{P}_{t-s}{S}_s\left(\widehat{\theta }\right){\nabla }_{\theta }L(s,S_s\left(\widehat{\theta }\right),\widehat{\theta }){\alpha }_{s}d{W}_s,\end{array}$$

(13)

where

$$P_t=\mathcal{E}\left({\int }_0^t\left(L(s,S_s\left(\widehat{\theta }\right),\widehat{\theta })+S_s\left(\widehat{\theta }\right){\nabla }_{s}L(s,S_s\left(\widehat{\theta }\right),\widehat{\theta })\right){\alpha }_{s}d{W}_s\right)$$

with $\mathcal{E}$ denoting the stochastic exponential.

Proof. 

First note that Theorem A1 implies the existence and uniqueness of

$$d{S}_t\left(\theta \right)=S_t\left(\theta \right)L(t,S_t\left(\theta \right),\theta ){\alpha }_{t}d{W}_t,$$

for every $\theta$. Here, the driving process is one-dimensional and given by $Y={\int }_0^{·}{\alpha }_{s}d{W}_s$. Indeed, according to Remark A1, if $(t,s)\mapsto L(t,s,\theta )$ is bounded, càdlàg in t and Lipschitz in s with a Lipschitz constant independent of t, $S_{·}\mapsto S_{·}\left(\theta \right)L(·,S_{·}\left(\theta \right),\theta )$ is functionally Lipschitz and Theorem A1 implies the assertion. These conditions are implied by the form of $L(t,s,\theta )$ and the conditions on the neural networks $F^i$.

To prove the form of the derivative process we apply Theorem A2 to the following system: consider

$$d{S}_t\left(\widehat{\theta }\right)=S_t\left(\widehat{\theta }\right)L(t,S_t\left(\widehat{\theta }\right),\widehat{\theta }){\alpha }_{t}d{W}_t,$$

together with

$$d{S}_t(\widehat{\theta }+\epsilon \theta )=S_t(\widehat{\theta }+\epsilon \theta )L(t,S_t(\widehat{\theta }+\epsilon \theta ),\widehat{\theta }+\epsilon \theta ){\alpha }_{t}d{W}_t,$$

as well as

$$\begin{array}{cc}\hfill d\frac{S_t(\widehat{\theta }+\epsilon \theta )-S_t\left(\widehat{\theta }\right)}{\epsilon }& =\frac{S_t(\widehat{\theta }+\epsilon \theta )L(t,S_t(\widehat{\theta }+\epsilon \theta ),\widehat{\theta }+\epsilon \theta )-S_t\left(\widehat{\theta }\right)L(t,S_t\left(\widehat{\theta }\right),\widehat{\theta })}{\epsilon }{\alpha }_{t}d{W}_t\hfill \\ \hfill & =(\frac{S_t(\widehat{\theta }+\epsilon \theta )-S_t\left(\widehat{\theta }\right)}{\epsilon }L(t,S_t(\widehat{\theta }+\epsilon \theta ),\widehat{\theta }+\epsilon \theta )\hfill \\ \hfill & +S_t\left(\widehat{\theta }\right)\frac{L(t,S_t(\widehat{\theta }+\epsilon \theta ),\widehat{\theta }+\epsilon \theta )-L(t,S_t\left(\widehat{\theta }\right),\widehat{\theta })}{\epsilon }){\alpha }_{t}d{W}_t.\hfill \end{array}$$

In the terminology of Theorem A2, $Z^{\epsilon ,1}=S\left(\widehat{\theta }\right)$, $Z^{\epsilon ,2}=S(\widehat{\theta }+\epsilon \theta )$ and $Z^{\epsilon ,3}=\frac{S_t(\widehat{\theta }+\epsilon \theta )-S_t\left(\widehat{\theta }\right)}{\epsilon }$. Moreover, $F^{\epsilon ,3}$ is given by

$$\begin{array}{cc}\hfill F^{\epsilon ,3}\left(Z_t^0\right)& =Z_t^{0,3}L(t,Z_t^{0,2},\widehat{\theta }+\epsilon \theta )+Z_t^{0,1}{\partial }_{s}L(t,Z_t^{0,1},\widehat{\theta })Z_t^{0,3}+\mathcal{O}\left(\epsilon \right)\hfill \\ \hfill & +Z_t^{0,1}\frac{L(t,Z_t^{0,1},\widehat{\theta }+\epsilon \theta )-L(t,Z_t^{0,1},\widehat{\theta })}{\epsilon },\hfill \end{array}$$

(14)

which converges ucp to

$$\begin{array}{c}\hfill F^{0,3}\left(Z_t^0\right)=Z_t^{0,3}L(t,Z_t^{0,2},\widehat{\theta })+Z_t^{0,1}{\partial }_{s}L(t,Z_t^{0,1},\widehat{\theta })Z_t^{0,3}+Z_t^{0,1}{\nabla }_{\theta }L(t,Z_t^{0,1},\widehat{\theta }).\end{array}$$

Indeed, for every fixed t, the family $\{s\mapsto L(t,s,\widehat{\theta }+\epsilon \theta ),|\epsilon \in [0,1]\}$ is due to the form of the neural networks equicontinuous. Hence pointwise convergence implies uniform convergence in s. This together with $L(t,s,\theta )$ being piecewise constant in t yields

$$\underset{\epsilon \to 0}{lim}\underset{(t,s)}{sup}|L(t,s,\widehat{\theta }+\epsilon \theta )-L(t,s,\widehat{\theta })|=0,$$

whence ucp convergence of the first term in (14). The convergence of term two is clear. The one of term three follows again from the fact that the family $\{s\mapsto {\nabla }_{\theta }L(t,s,\widehat{\theta }+\epsilon \theta )|\epsilon \in [0,1]\}$ is equicontinuous, which is again a consequence of the form of the neural networks.

By the assumptions on the derivatives, $F^{0,3}$ is functionally Lipschitz. Hence Theorem A1 yields the existence of a unique solution to (12) and Theorem A2 implies convergence. □

Remark 1.

1. 

For the pure existence and uniqueness of

$$d{S}_t\left(\theta \right)=S_t\left(\theta \right)L(t,S_t\left(\theta \right),\theta ){\alpha }_{t}d{W}_t,$$

with $L(t,s,\theta )$ of form (11), it suffices that the neural networks $s\mapsto F^i(s,{\theta }_i)$ are bounded and Lipschitz, for all $i=1,\dots ,n$ (see also Remark A1).

2. 

Formula (13) can be used for well-known backward propagation schemes.

Theorem 2 guarantees the existence and uniqueness of the derivative process. This thus allows the setting up of gradient-based search algorithms for training.

In view of this let us now come to the precise optimization task as already outlined in Section 1.1. To ease the notation, we shall here omit the dependence of the weights w and the loss function ℓ on the parameter $\gamma$. For each maturity $T_i$, we assume to have $J_i$ options with strikes $K_{ij}$, $j\in \{1,\dots ,J_i\}$. The calibration functional for the i-th maturity is then of the form

$$\begin{array}{c}\hfill \underset{{\theta }_i\in {\Theta }_i}{argmin}\sum _{j=1}^{J_i}{w}_{ij}\ell ({\pi }_{ij}^{\mathrm{mod}}\left({\theta }_i\right)-{\pi }_{ij}^{\mathrm{mkt}}),i\in \{1,\dots ,n\}.\end{array}$$

(15)

Recall from the introduction that ${\pi }_{ij}^{\mathrm{mod}}\left({\theta }_i\right)$ (${\pi }_{ij}^{\mathrm{mkt}}$ respectively) denotes the model (market resp.) price of an option with maturity $T_i$ and strike $K_{ij}$. Moreover, $\ell :\mathbb{R}\to {\mathbb{R}}_{+}$ is some non-negative, nonlinear, convex loss function (e.g., square or absolute value) with $\ell \left(0\right)=0$ and $\ell \left(x\right)>0$ for $x\ne 0$, measuring the distance between market and model prices. Finally, $w_{ij}$ denote some weights, e.g., of vega type (compare Cont and Ben Hamida (2004)), which we use to match implied volatility data rather than pure prices. Notice that we here omit for notational convenience the dependence of $w_{ij}$ and ℓ on parameters $\gamma$ which describe the adversarial part.

We solve the minimization problems (15) iteratively: we start with maturity $T_1$ and fix ${\theta }_1$. This then enters in the computation of ${\pi }_{2j}^{\mathrm{mod}}\left({\theta }_2\right)$ and thus in (15) for maturity $T_2$, etc. To simplify the notation in the sequel, we shall therefore leave the index i away so that for a generic maturity $T>0$, (15) becomes

$$\begin{array}{c}\hfill \underset{\theta \in \Theta }{argmin}\sum _{j=1}^J{w}_j\ell ({\pi }_j^{\mathrm{mod}}\left(\theta \right)-{\pi }_j^{\mathrm{mkt}}).\end{array}$$

Since the model prices are given by

$$\begin{array}{c}\hfill {\pi }_j^{\mathrm{mod}}\left(\theta \right)=\mathbb{E}\left[{(S_T\left(\theta \right)-K_j)}^{+}\right],\end{array}$$

(16)

we have ${\pi }_j^{\mathrm{mod}}\left(\theta \right)-{\pi }_j^{\mathrm{mkt}}=\mathbb{E}\left[Q_j\left(\theta \right)\right]$ where

$$\begin{array}{c}\hfill Q_j\left(\theta \right)\left(\omega \right):={(S_T\left(\theta \right)\left(\omega \right)-K_j)}^{+}-{\pi }_j^{\mathrm{mkt}}.\end{array}$$

(17)

The calibration task then amounts to finding a minimum of

$$\begin{array}{c}\hfill f\left(\theta \right):=\sum _{j=1}^J{w}_j\ell \left(\mathbb{E}\left[Q_j\left(\theta \right)\right]\right).\end{array}$$

(18)

As ℓ is a nonlinear function, this is not of the expected value form of problem (A3). Hence standard stochastic gradient descent, as outlined in Appendix B.2, cannot be applied in a straightforward manner.

We shall tackle this problem via hedge control variates as introduced in Section 2. In the following we explain this in more detail.

Minimizing the Calibration Functional

Consider the standard Monte Carlo estimator for $\mathbb{E}\left[Q_j\left(\theta \right)\right]$ so that (18) is estimated by

$$\begin{array}{c}\hfill f^{\mathrm{MC}}\left(\theta \right):=\sum _{j=1}^J{w}_j\ell \left(\frac{1}{N}\sum _{n=1}^N{Q}_j\left(\theta \right)\left({\omega }_n\right)\right),\end{array}$$

(19)

for i.i.d samples $\{{\omega }_1,\dots ,{\omega }_N\}\in \Omega$. Since the Monte Carlo error decreases as $\frac{1}{\sqrt{N}}$, the number of simulations N must be chosen large ($\approx {10}^8$) to approximate well the true model prices in (16). Note that implied volatility to which we actually aim to calibrate is even more sensitive. As stochastic gradient descent is not directly applicable due to the nonlinearity of ℓ, it seems necessary at first sight to compute the gradient of the whole function $\widehat{f}\left(\theta \right)$ to minimize (19). As $N\approx {10}^8$, this is however computationally very expensive and leads to numerical instabilities as we must compute the gradient of a sum that contains ${10}^8$ terms. Hence with this method an (approximative) minimum in the high-dimensional parameter space $\Theta$ cannot be found in a reasonable amount of time.

One very expedient remedy is to apply hedge control variates as introduced in Section 2 as variance reduction technique. This allows the reduction of the number of samples N in the Monte Carlo estimator considerably to only up to $5\times {10}^4$ sample paths.

Assume that we have r hedging instruments (including the price process S) denoted by ${\left(Z_t\right)}_{t\in [0,T]}$ which are square integrable martingales under $\mathbb{Q}$ and take values in ${\mathbb{R}}^r$. Consider, for $j=1,\dots ,J$, strategies $h_j:[0,T]\times {\mathbb{R}}^r\to {\mathbb{R}}^r$ such that $h(·,Z_{·})\in L^2\left(Z\right)$ and some constant c. Define

$$\begin{array}{c}\hfill X_j\left(\theta \right)\left(\omega \right):=Q_j\left(\theta \right)\left(\omega \right)-c{(h_j(·,Z_{·-}\left(\theta \right)\left(\omega \right))•Z_{·}\left(\theta \right)\left(\omega \right))}_T\end{array}$$

(20)

The calibration functionals (18) and (19), can then simply be defined by replacing $Q_j\left(\theta \right)\left(\omega \right)$ by $X_j\left(\theta \right)\left(\omega \right)$ so that we end up minimizing

$$\begin{array}{c}\hfill \widehat{f}\left(\theta \right)({\omega }_1,\dots ,{\omega }_N)=\sum _{j=1}^J{w}_j\ell \left(\frac{1}{N}\sum _{n=1}^N{X}_j\left(\theta \right)\left({\omega }_n\right)\right).\end{array}$$

(21)

To tackle this task, we apply the following variant of gradient descent: starting with an initial guess ${\theta }^{\left(0\right)}$, we iteratively compute

$$\begin{array}{c}\hfill {\theta }^{(k+1)}={\theta }^{\left(k\right)}-{\eta }_{k}G\left({\theta }^{\left(k\right)}\right)({\omega }_1^{\left(k\right)},\dots ,{\omega }_N^{\left(k\right)}),\end{array}$$

(22)

for some learning rate ${\eta }_k$, i.i.d samples $({\omega }_1^{\left(k\right)},\dots ,{\omega }_N^{\left(k\right)})$, where the values

$$G\left({\theta }^{\left(k\right)}\right)({\omega }_1^{\left(k\right)},\dots ,{\omega }_N^{\left(k\right)})$$

are gradient-based quantities that remain to be specified. These samples can either be chosen to be the same in each iteration or to be newly sampled in each update step. The difference between these two approaches is negligible, since N is chosen so as to yield a small Monte Carlo error, whence the gradient is nearly deterministic. In our numerical experiments we newly sample in each update step.

In the simplest form, one could simply set

$$G\left({\theta }^{\left(k\right)}\right)({\omega }_1^{\left(k\right)},\dots ,{\omega }_N^{\left(k\right)})=\nabla \widehat{f}\left(\theta \right)({\omega }_1^{\left(k\right)},\dots ,{\omega }_N^{\left(k\right)}).$$

(23)

Note however that the derivative of the stochastic integral term in (20) is in general quite expensive. We thus implement the following modification.

We set

$$\begin{array}{cc}\hfill {\omega }^N& =({\omega }_1,\dots ,{\omega }_N),\hfill \\ \hfill Q_j^N\left(\theta \right)\left({\omega }^N\right)& =\frac{1}{N}\sum _{n=1}^N{Q}_j\left(\theta \right)\left({\omega }_n\right),\hfill \\ \hfill Q^N\left(\theta \right)\left({\omega }^N\right)& =(Q_1^N\left(\theta \right)\left({\omega }^N\right),\dots ,Q_J^N\left(\theta \right)\left({\omega }^N\right)),\hfill \end{array}$$

and define $\tilde{f}:{\mathbb{R}}^J\to \mathbb{R}$ via

$$\tilde{f}\left(x\right)=\sum _{j=1}^J{w}_j\ell \left(x_j\right).$$

We then set

$$G\left(\theta \right)\left({\omega }^N\right)=D_x\left(\tilde{f}\right)\left(X^N\left(\theta \right)\left({\omega }^N\right)\right)D_{\theta }\left(Q^N\right)\left(\theta \right)\left({\omega }^N\right).$$

Please note that this quantity is actually easy to compute in terms of backpropagation. Moreover, leaving the stochastic integral away in the inner derivative is justified by its vanishing expectation. During the forward pass, the stochastic integral terms are included in the computation; however the contribution to the gradient (during the backward pass) is partly neglected, which can e.g.,  be implemented via the tensorflow stop_gradient function.

Concerning the choice of the hedging strategies, we can parameterize them as in Section 2.2 via neural networks and find the optimal weights $\delta$ by computing

$$\begin{array}{c}\hfill \underset{\delta \in \Delta }{argmin}\frac{1}{N}\sum _{n=1}^{N}u(-X_j(\theta ,\delta )\left({\omega }_n\right)).\end{array}$$

(24)

for i.i.d samples $\{{\omega }_1,\dots ,{\omega }_N\}\in \Omega$ and some loss function u when $\theta$ is fixed. Here,

$$X_j(\theta ,\delta )\left(\omega \right)={(S_T\left(\theta \right)\left(\omega \right)-K_j)}^{+}-{(h_j(·,Z_{·-}\left(\theta \right)\left(\omega \right),\delta )•Z_{·}\left(\theta \right)\left(\omega \right))}_T-{\pi }_j^{\mathrm{mkt}}.$$

This means to iterate the two optimization procedures, i.e., minimizing (21) for $\theta$ (with fixed $\delta )$ and (24) for $\delta$ (with fixed $\theta$). Clearly the Black–Scholes hedge ansatz as of Section 2.1 works as well, in this case without additional optimization with respect to the hedging strategies.

For alternative approaches how to minimize (18), we refer to Appendix C.

4. Numerical Implementation

In this section, we discuss the numerical implementation of the proposed calibration method. We implement our approach via tensorflow, taking advantage of GPU-accelerated computing. All computations are performed on a single-gpu Nvidia GeForce^® GTX 1080 Ti machine. For the implied volatility computations, we rely on the python py_vollib library.3

Recall that a LSV model is given on some filtered probability space $(\Omega ,{\left({\mathcal{F}}_t\right)}_{t\in [0,T]},\mathcal{F},\mathbb{Q})$ by

$$d{S}_t=S_t{\alpha }_{t}L(t,S_t)d{W}_t,S_0>0,$$

for some stochastic process $\alpha$. When calibrating to data, it is, therefore, necessary to make further specifications. We calibrate the following SABR-type LSV model.

Definition 1.

The SABR-LSV model is specified via the SDE,

$$\begin{array}{cc}\hfill d{S}_t& =S_{t}L(t,S_t){\alpha }_{t}d{W}_t,\hfill \\ \hfill d{\alpha }_t& =\nu {\alpha }_{t}d{B}_t,\hfill \\ \hfill d{\langle W,B\rangle }_t& =\rho dt,\hfill \end{array}$$

with parameters $\nu \in \mathbb{R}$, $\rho \in [-1,1]$ and initial values ${\alpha }_0>0,S_0>0$. Here, B and W are two correlated Brownian motions.

Remark 2.

We shall often work in log-price coordinates for S. In particular, we can then consider L as a function of $X:=logS$ rather then S. By denoting this parametrization again with L, we therefore have $L(t,X)$ instead of $L(t,S)$ and the model dynamics read as

$$\begin{array}{cc}\hfill d{X}_t& ={\alpha }_{t}L(t,X_t)d{W}_t-\frac{1}{2}{\alpha }_t^2{L}^2(t,X_t)dt,\hfill \\ \hfill d{\alpha }_t& =\nu {\alpha }_{t}d{B}_t,\hfill \\ \hfill d{\langle W,B\rangle }_t& =\rho dt.\hfill \end{array}$$

Please note that α is a geometric Brownian motion, in particular, the closed form solution for α is available and given by

$${\alpha }_t={\alpha }_{0}exp\left(-\frac{{\nu }^2}{2}t+\nu B_t\right).$$

For the rest of the paper we shall set $S_0=1$.

4.1. Implementation of the Calibration Method

We now present a proper numerical test and demonstrate the effectiveness of our approach on a family of typical market smiles (instead of just one calibration example). We consider as ground truth a situation where market smiles are produced by a parametric family. By randomly sampling smiles from this family we then show that they can be calibrated up to small errors, which we analyze statistically.

4.1.1. Ground Truth Assumption

We start by specifying the ground truth assumption. It is known that a discrete set of prices can be exactly calibrated by a local volatility model using Dupire’s volatility function, if an appropriate interpolation method is chosen. Hence, any market observed smile data can be reproduced by the following model (we assume zero riskless rate and define $X=log\left(S\right)$),

$$d{S}_t={\sigma }_{\mathrm{Dup}}(t,X_t)S_{t}d{W}_t,$$

or equivalently

$$d{X}_t=-\frac{1}{2}{\sigma }_{\mathrm{Dup}}^2(t,X_t)dt+{\sigma }_{\mathrm{Dup}}(t,X_t)d{W}_t,$$

(25)

where ${\sigma }_{\mathrm{Dup}}$ denotes Dupire’s local volatility function Dupire (1996). Our ground truth assumption consists of supposing that the function ${\sigma }_{\mathrm{Dup}}$ (or to be more precise ${\sigma }_{\mathrm{Dup}}^2$) can be chosen from a parametric family. Such parametric families for local volatility models have been discussed in the literature, consider e.g., Carmona and Nadtochiy (2009) or Carmona et al. (2007). In the latter, the authors introduce a family of local volatility functions ${\tilde{a}}_{\xi }$ indexed by parameters

$$\xi =(p_1,p_2,{\sigma }_0,{\sigma }_1,{\sigma }_2)$$

and $p_0=1-(p_1+p_2)$ satisfying the constraints

$${\sigma }_0,{\sigma }_1,{\sigma }_2,p_1,p_2>0\mathrm{and}{p}_1+p_2\le 1.$$

Setting $k(t,x,\sigma )=exp\left(-x^2/\left(2t{\sigma }^2\right)-t{\sigma }^2/8\right)$, ${\tilde{a}}_{\xi }$ is then defined as

$${\tilde{a}}_{\xi }^2(t,x)=\frac{{\sum }_{i=0}^2{p}_i{\sigma }_{i}k(t,x,{\sigma }_i)}{{\sum }_{i=0}^2(p_i/{\sigma }_i)k(t,x,{\sigma }_i)}.$$

In Figure 1a we show plots of implied volatilities for different slices (maturities) for a realistic choice of parameters. As one can see, the produced smiles seem to be unrealistically flat. Hence we modify the local volatility function ${\tilde{a}}_{\xi }$ to produce more pronounced and more realistic smiles. To be precise, we define a new family of local volatility functions $a_{\xi }$ indexed by the set of parameters $\xi$ as

$$a_{\xi }^2(t,x)=\frac{1}{4}\times min(2,\left|\frac{\left({\sum }_{i=0}^2{p}_i{\sigma }_{i}k(t,x,{\sigma }_i)+\Lambda (t,x)\right)\left(1-0.6\times {\mathbb{1}}_{\left(t>0.1\right)}\right)}{{\sum }_{i=0}^2(p_i/{\sigma }_i)k(t,x,{\sigma }_i)+0.01}\right|),$$

(26)

with

$$\Lambda (t,x):={\left(\frac{{\mathbb{1}}_{\left(t\le 0.1\right)}}{1+0.1t}\right)}^{{\lambda }_2}min\left\{{\left({\gamma }_1{\left(x-{\beta }_1\right)}_{+}+{\gamma }_2{\left(-x-{\beta }_2\right)}_{+}\right)}^{\kappa },{\lambda }_1\right\}.$$

We fix the choice of the parameters ${\gamma }_i,{\beta }_i,{\lambda }_i,\kappa$ as given in

Table 1. Fixed Parameters for the ground truth assumption $a_{\xi }^2$.

${\mathit{\gamma }}_1$	${\mathit{\gamma }}_2$	${\mathit{\lambda }}_1$	${\mathit{\lambda }}_2$	${\mathit{\beta }}_1$	${\mathit{\beta }}_2$	$\mathit{\kappa }$
$1.1$	20	10	10	$0.005$	$0.001$	$0.5$

. By taking absolute values above, we can drop the requirement $p_0>0$ which is what we do in the sequel. Please note that $a_{\xi }^2$ is not defined at $t=0$. When doing a Monte Carlo simulation, we simply replace $a_{\xi }^2(0,x)$ with $a_{\xi }^2({\Delta }_t,x)$, where ${\Delta }_t$ is the time increment of the Monte Carlo simulation.

What is left to be specified are the parameters

$$\xi =(p_1,p_2,{\sigma }_0,{\sigma }_1,{\sigma }_2)$$

with $p_0=1-p_1-p_2$. This motivates our statistical test for the performance evaluation of our method. To be precise, our ground truth assumption is that all observable market prices are explained by a variation of the parameters $\xi$. For illustration, we plot implied volatilities for this modified local volatility function in Figure 1b for a specific parameter set $\xi$.

Our ground truth model is now specified as in (25) with ${\sigma }_{\mathrm{Dup}}$ replaced by $a_{\xi }$, i.e.,

$$d{X}_t=-\frac{1}{2}{a}_{\xi }^2(t,X_t)dt+a_{\xi }(t,X_t)d{W}_t.$$

(27)

4.1.2. Performance Test

We now come to the evaluation of our proposed method. We want to calibrate the SABR-LSV model to synthetic market prices generated by the previously formulated ground truth assumption. This corresponds to randomly sampling the parameter $\xi$ of the local volatility function $a_{\xi }$ and to compute prices according to (27). Calibrating the SABR-LSV model, i.e., finding the parameters $\nu ,\rho$, the initial volatility ${\alpha }_0$ and the unknown leverage function L, to these prices and repeating this multiple times then allows for a statistical analysis of the errors.

As explained in Section 3, we consider European call options with maturities $T_1<\cdots <T_n$ and denote the strikes for a given maturity $T_i$ by $K_{ij}$, $j\in \{1,\dots ,J_i\}$. To compute the ground truth prices for these European calls we use a Euler-discretization of (27) with time step ${\Delta }_t=1/100$. Prices are then obtained by a variance reduced Monte Carlo estimator using ${10}^7$ Brownian paths and a Black–Scholes delta hedge variance reduction as described previously. For a given parameter set $\xi$, we use the same Brownian paths for all strikes and maturities.

Overall, in this test, we consider $n=4$ maturities with $J_i=20$ strike prices for all $i=1,\dots ,4$. The values for $T_i$ are given in Figure 2a. For the choice of the strikes $K_i$, we choose evenly spaced points, i.e.,

$$K_{i,j+1}-K_{i,j}=\frac{K_{i,20}-K_{i,1}}{19}.$$

For the smallest and largest strikes per maturity we choose

$$K_{i,1}=exp\left(-k_i\right),K_{i,20}=exp\left(k_i\right),$$

with the values of $k_i$ given in Figure 2b.

We now specify a distribution under which we draw the parameters

$$\xi =(p_1,p_2,{\sigma }_0,{\sigma }_1,{\sigma }_2,)$$

for our test. The components are all drawn independently from each other under the uniform distribution on the respective intervals given below.

- $I_{p_1}=[0.4,0.5]$ - $I_{p_2}=[0.4,0.7]$

- $I_{{\sigma }_0}=[0.5,1.7]$ - $I_{{\sigma }_1}=[0.2,0.4]$

- $I_{{\sigma }_2}=[0.5,1.7]$

We can now generate data by the following scheme.

• For $m=1,\dots ,200$ simulate parameters ${\xi }_m$ under the law described above.
• For each m, compute prices of European calls for maturities $T_i$ and strikes $K_{ij}$ for $i=1,\dots ,n=4$ and $j=1,\dots ,20$ according to (27) using ${10}^7$ Brownian trajectories (for each m we use new trajectories).
• Store these prices.

Remark 3.

In very few cases, the simulated parameters were such that the implied volatility computation for model prices failed at least for one maturity due to the remaining Monte Carlo error. In those cases, we simply skip that sample and continue with the next, meaning that we will perform the statistical test only on the samples for which these implied volatility computations were successful.

The second part consists of calibrating each of these surfaces and storing pertinent values for which we conduct a statistical analysis. In the following we describe the procedure in detail:

Recall that we specify the leverage function $L(t,x)$ via a family of neural networks, i.e.,

$$\begin{array}{c}\hfill L(t,x)=1+F^i\left(x\right)t\in [T_{i-1},T_i),i\in \{1,\dots ,n=4\},\end{array}$$

where $F^i\in {\mathcal{NN}}_{1,1}$ (see Notation A1). Each $F^i$ is specified as a 4-hidden layer feed-forward network where the dimension of each of the hidden layers is 64. As activation function we choose leaky-ReLU4 with parameter $0.2$ for the first three hidden layers and $\varphi =tanh$ for the last hidden layer. This choice means of course a considerable overparameterization, where we deal with much more parameters than data points. As is well known from the theory of machine learning, this however allows a profit to be made from implicit regularizations for the leverage function, meaning that the variations of higher derivatives are small.

Remark 4.

In our experiments, we tested different network architectures. Initially, we used networks with three to five hidden layers with layer dimensions between 50 and 100 and activation function tanh in all layers. Although the training was successful, we observed that training was significantly slower with significant lower calibration accuracy compared to the final architecture. We also tried classical ReLU, but observed that the training sometimes got stuck due to flat gradients. In case of pure leaky-ReLU activation functions, we observed numerical instabilities. By adding a final tanh activation, this computation was regularized leading to the results we present here.

Since closed form pricing formulas are not available for such an LSV model, let us briefly specify our pricing method. For the variance reduced Monte Carlo estimator as of (21) we always use a standard Euler-SDE discretization with step size ${\Delta }_t=1/100$. As variance reduction method, we implement the running Black–Scholes Delta hedge with instantaneous running volatility of the price process, i.e., $L(t,X_t){\alpha }_t$ is plugged in the formula for the Black–Scholes Delta as in (9). The only parameter that remains to be specified, is the number of trajectories used for the Monte Carlo estimator which is done in Algorithms A1 and A2 below.

As a first calibration step, we calibrate the SABR model (i.e., (1) with $L\equiv 1$) to the synthetic market prices of the first maturity and fix the calibrated SABR parameters $\nu ,\rho$ and ${\alpha }_0$. This calibration is not done by the SABR formula, but rather in the same way the LSV model calibration is implemented: we use a Monte Carlo simulation based engine where gradients are computed via backpropagation. The calibration objective function is analog to (21) and we compute the full gradient as specified in (23). We only use a maximum of 2000 trajectories and the running Black–Scholes hedge for variance reduction per gradient computation, as we are only interested in an approximate fit. In fact, when compared to a better initial SABR fit achieved by the SABR formula, we observed that the calibration fails more often due to local minima becoming an issue.

For training the parameters ${\theta }_i$, $i=1,\dots ,4$, of the neural networks we apply Algorithm A1 in the Appendix D.

4.2. Numerical Results for the Calibration Test

We now discuss the results of our test. We start by pointing out that from the 200 synthetic market smiles generated, four smiles caused difficulties, in the sense that our implied volatility computation failed due to the remaining Monte Carlo error in the model price computation, compare Remark 3. By increasing the training parameters slightly (in particular the number of trajectories used in the training), this issue can be mitigated but the resulting calibrated implied volatility errors stay large out of the money where the smiles are extreme, and the training will take more time. Hence, we opt to remove those four samples from the following statistical analysis as they represented unrealistic market smiles.

In Figure 3 we show calibration results for a typical example of randomly generated synthetic market data. From this it is already visible that the worst-case calibration error (which occurs out of the money) ranges typically between 5 and 15 basis points. The corresponding calibration result for the square of the leverage function $L^2$ is given in Figure 4.

Let us note that our method achieves a very high calibration accuracy for the considered range of strikes across all considered maturities. This can be seen in the results of a worst-case analysis of calibration errors in Figure 5. There we show the mean as well as different quantiles of the data. Please note that the mean always lies below 10 basis point across all strikes and maturities.

Regarding calibration times, we can report that from the 196 samples, 191 finished within 26 to 27 min. In all these cases, the abort criterion was active on the first time it was checked, i.e., after 5000 iterations. The other five samples are examples of smiles comparable to the four where implied volatility computation itself failed. In those cases, more iteration steps where needed resulting in times between 46 and 72 min. These samples also correspond to the less successful calibration results.

To perform an out of sample analysis, we check for extra- and interpolation properties of the learned leverage function. This means that we compute implied volatilities on an extended range and compare to the implied volatility of the ground truth assumption. The strikes of these ranges are again computed by taking 20 equally spaced points as before, but with parameters $k_i$ as of table Figure 2b multiplied with 1.5. This has also the effect that the strikes inside the original range do not correspond to the strikes considered during training, which allows for an additional analysis of the interpolation properties. These results are illustrated in Figure 6, from which we see that extrapolation is very close to the local volatility model.

4.3. Robust Calibration—An Instance of the Adversarial Approach

Let us now describe a robust version of our calibration methodology realized in an adversarial manner. We start by assuming that there are multiple “true” option prices which correspond to the bid-ask spreads observed on the market. The way we realize this in our experiment is to use several local volatility functions that generate equally plausible market implied volatilities. Recall that the local volatility functions in our statistical test above are functions of the parameters $(p_0,p_1,{\sigma }_0,{\sigma }_1,{\sigma }_2)$. We fix these parameters and generate 4 smiles from local volatility functions with slightly perturbed parameters

$$(p_0+u_{i1},p_1+u_{i2},{\sigma }_0+u_{i3},{\sigma }_1+u_{i4},{\sigma }_2+u_{i5})\mathrm{for}i=1,\dots ,4,$$

where $u_{ij}$ are i.i.d. uniformly distributed random variables, i.e., $u_{ij}\sim {\mathcal{U}}_{[-u,u]}$ with $u=0.01$. The loss function for maturity $T_i$ in the training part now changes to

$$\underset{\theta }{inf}\sum _{j=1}^{J_i}{w}_j\underset{m=1,\dots ,4}{sup}\ell \left(\frac{1}{N}\sum _{n=1}^N{X}_{j,m}\left(\theta \right)\left({\omega }_n\right)\right),$$

(28)

with $X_{j,m}$ defined as $X_j$ in (20) (see also (17)) but with synthetic market prices $m=1,\dots ,4$ generated by the m-th local volatility function. We are thus in an adversarial situation as described in the introduction: we have several possibilities for the loss function corresponding to the different market prices and we take the supremum over these (individually for each strike). In our toy example we can simply compute the gradient of this supremum function with respect to $\theta$. In a more realistic situation, where we do not only have 4 smiles but a continuum we would iterate the inf and sup computation, meaning that we would also perform a gradient step with respect to m. This corresponds exactly to the adversary part. For a given parameter set $\theta$, the adversary tries to find the worst loss function.

In Figure 7, we illustrate the result of this robust calibration, where find that the calibrated model lies between the four different smiles over which we take the supremum.

5. Conclusions

We have demonstrated how the parametrization by means of neural networks can be used to calibrate local stochastic volatility models to implied volatility data. We make the following remarks:

• The method we presented does not require any form of interpolation for the implied volatility surface since we do not calibrate via Dupire’s formula. As the interpolation is usually done ad hoc, this might be a desirable feature of our method.
• Similar to Guyon and Henry-Labordѐre (2012); Guyon and Henry-Labordѐre (2013), it is possible to “plug in” any stochastic variance process such as rough volatility processes as long as an efficient simulation of trajectories is possible.
• The multivariate extension is straight forward.
• The level of accuracy of the calibration result is of a very high degree. The average error in our statistical test is of around 5 to 10 basis points, which is an interesting feature in its own right. We also observe good extrapolation and generalization properties of the calibrated leverage function.
• The method can be significantly accelerated by applying distributed computation methods in the context of multi-GPU computational concepts.
• The presented algorithm is further able to deal with path-dependent options since all computations are done by means of Monte Carlo simulations.
• We can also consider the instantaneous variance process of the price process as short end of a forward variance process, which is assumed to follow (under appropriate assumptions) a neural SDE. This setting, as an infinite-dimensional version of the aforementioned “multivariate” setting, then qualifies for joint calibration to S&P and VIX options. This is investigated in a companion paper.
• We stress again the advantages of the generative adversarial network point of view. We believe that this is a crucial feature in the joint calibration of S&P and VIX options.

6. Plots

This section contains the relevant plots for the numerical test outlined in Section 4.