A generative adversarial network approach to calibration of local stochastic volatility models

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, b)) 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.

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 $(S_{t})_{t\geq 0}$ of an asset satisfies

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

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 Guyon and Henry-Labordere (2012); Guyon and Henry-Labordère (2013); Cozma et al. (2017) 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-Labordere (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.

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(\cdot,\cdot,\theta)$ and volatility $\sigma(\cdot,\cdot,\theta)$ by neural networks with parameters $\theta$, i.e.,

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

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

where $J$ is the number of considered options and where $\pi^{\text{mod}}_{j}(\theta)$ and $\pi^{\text{mkt}}_{j}$ stand for the respective model and market prices.

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^{\text{mod}}(\theta)$ 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-Labordere (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-Labordere (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-Labordere (2012) (compare Section 6 and the results in Guyon and Henry-Labordere (2012)). Moreover, we also observe good extrapolation and generalization properties of the calibrated leverage function.

2. Generative Adversarial Approaches in Finance

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

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;

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,\text{market 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!).

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 (1.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 5. 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 (1.5).

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 Vidales et al. (2018); Potters et al. (2001).

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

where $(I_{1},\ldots,I_{N})$ are i.i.d with the same distribution as $I$. Then, for any $(h_{t})_{t\in[0,T]}\in L^{2}(Z)$ 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

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

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

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

In many cases, of local stochastic volatility models as of form (1.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(S_{T})$ and let $\pi^{g}_{\text{BS}}(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

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

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

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(Z_{T})$. Then in Markov models it makes sense to specify the hedging strategy via a function

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 NotationWe here use $\delta$ to denote the parameters of the hedging neural networks, as $\theta$ shall be used for the networks of the leverage function. B.4). Following the approach in (Bühler et al., 2019, Remark 3), an optimal hedge for the claim $C$ with given market price $\pi^{\text{mkt}}$ can be computed via

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

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

As always in this article we shall assume that activation functions $\phi$ 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(\delta)$ is given by

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 Section A.

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

where we obtain existence, uniqueness and stability for the second equation by Theorem A.3, 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.

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

Then the derivative $\nabla_{\delta}(h(\cdot,Z_{\cdot-},\delta)\bullet Z)$ in direction $\delta$ at $\delta_{0}$ satisfies

If additionally the derivative in direction $\delta$ at $\delta_{0}$ of $\nabla_{\delta}h^{\varepsilon}(\cdot,Z_{\cdot-},\delta_{0})$ converges ucp to $\nabla_{\delta}h(\cdot,Z_{\cdot-},\delta_{0})$ as $\varepsilon\to 0$, then the directional derivative of the discretized integral, i.e. $\nabla_{\delta}(h^{\varepsilon}(\cdot,Z^{\varepsilon}_{\cdot-},\delta_{0})\bullet Z^{\varepsilon})$ or equivalently $(\nabla_{\delta}h^{\varepsilon}(\cdot,Z^{\varepsilon}_{\cdot-},\delta_{0})\bullet Z^{\varepsilon})$, converges, as the discretization mesh $\varepsilon\to 0$, to

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

by equicontinuity of $\{(t,z)\mapsto\nabla_{\delta}h(t,z,\delta_{0}+\varepsilon\delta)\,|\,\varepsilon\in\}$.

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

which converges by assumption ucp to $f^{0}=\nabla_{\delta}h(\cdot,Z_{\cdot-},\delta_{0})$.

Calibration of LSV Models

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 B.3) 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.,

where $F^{i}\in\mathcal{NN}_{1,1}$ for $i=1,\ldots,n$ (see Notation B.4). For notational simplicity we shall often omit the dependence on $\theta_{i}\in\Theta_{i}$. However, when needed we write for instance $S_{t}(\theta)$, 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(\theta)$, which we prove here rigorously by applying the results of Appendix A.

Let $(t,s,\theta)\mapsto L(t,s,\theta)$ be of form (3.1) where the neural networks $(s,\theta_{i})\mapsto F^{i}(s,\theta_{i})$ are bounded and $C^{1}$, with bounded and Lipschitz continuous derivativesThis just means that the activation function is bounded and $C^{1}$, with bounded and Lipschitz continuous derivatives., for all $i=1,\ldots,n$. Then the directional derivative in direction $\theta$ at $\widehat{\theta}$ satisfies the following equation

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

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

First note that Theorem A.2 implies the existence and uniqueness of

for every $\theta$. Here, the driving process is one-dimensional and given by $Y=\int_{0}^{\cdot}\alpha_{s}dW_{s}$. Indeed, according to Remark A.4, 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_{\cdot}\mapsto S_{\cdot}(\theta)L(\cdot,S_{\cdot}(\theta),\theta)$ is functionally Lipschitz and Theorem A.2 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 A.3 to the following system: consider

In the terminology of Theorem A.3, $Z^{\varepsilon,1}=S(\widehat{\theta})$, $Z^{\varepsilon,2}=S(\widehat{\theta}+\varepsilon\theta)$ and $Z^{\varepsilon,3}=\frac{S_{t}(\widehat{\theta}+\varepsilon\theta)-S_{t}(\widehat{\theta})}{\varepsilon}$. Moreover, $F^{\varepsilon,3}$ is given by

Indeed, for every fixed $t$, the family $\{s\mapsto L(t,s,\widehat{\theta}+\varepsilon\theta),\,|\,\varepsilon\in\}$ 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

whence ucp convergence of the first term in (3.4). 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}+\varepsilon\theta)\,|\,\varepsilon\in\}$ 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 A.2 yields the existence of a unique solution to (3.2) and Theorem A.3 implies convergence. ∎

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

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

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

We solve the minimization problems (3.5) iteratively: we start with maturity $T_{1}$ and fix $\theta_{1}$. This then enters in the computation of $\pi^{\text{mod}}_{2j}(\theta_{2})$ and thus in (3.5) 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$, (3.5) becomes

The calibration task then amounts to finding a minimum of

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

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.

The calibration functionals (3.8) and (3.9), can then simply be defined by replacing $Q_{j}(\theta)(\omega)$ by $X_{j}(\theta)(\omega)$ so that we end up minimizing

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

for some learning rate $\eta_{k}$, i.i.d samples $(\omega^{(k)}_{1},\ldots,\omega^{(k)}_{N})$, where the values

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

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

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

for i.i.d samples $\{\omega_{1},\ldots,\omega_{N}\}\in\Omega$ and some loss function $u$ when $\theta$ is fixed. Here,

This means to iterate the two optimization procedures, i.e., minimizing (3.11) for $\theta$ (with fixed $\delta)$ and (3.14) 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 (3.8), we refer to Appendix C.

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.See http://vollib.org/.

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.

The SABR-LSV model is specified via the SDE,

We shall often work in log-price coordinates for $S$. In particular, we can then consider $L$ as a function of $X:=\log S$ 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

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

For the rest of the paper we shall set $S_{0}=1$.

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.

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(S)$),

where $\sigma_{\text{Dup}}$ denotes Dupire’s local volatility function Dupire (1996). Our ground truth assumption consists of supposing that the function $\sigma_{\text{Dup}}$ (or to be more precise $\sigma_{\text{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 $\widetilde{a}_{\xi}$ indexed by parameters

and $p_{0}=1-(p_{1}+p_{2})$ satisfying the constraints

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

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 $\widetilde{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

We fix the choice of the parameters $\gamma_{i},\beta_{i},\lambda_{i},\kappa$ as given in Table 1. 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^{2}_{\xi}$ is not defined at $t=0$. When doing a Monte Carlo simulation, we simply replace $a^{2}_{\xi}(0,x)$ with $a^{2}_{\xi}(\Delta_{t},x)$, where $\Delta_{t}$ is the time increment of the Monte Carlo simulation.

What is left to be specified are the parameters

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 (4.1) with $\sigma_{\text{Dup}}$ replaced by $a_{\xi}$, i.e.,

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 (4.3). Calibrating the SABR-LSV model, i.e., finding the parameters $\nu,\varrho$, 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,\ldots,J_{i}\}$. To compute the ground truth prices for these European calls we use a Euler-discretization of (4.3) 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,\ldots,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.,

For the smallest and largest strikes per maturity we choose

with the values of $k_{i}$ given in Figure 2b.

We now specify a distribution under which we draw the parameters

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

We can now generate data by the following scheme.

For $m=1,\ldots,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,\ldots,n=4$ and $j=1,\ldots,20$ according to (4.3) using $10^{7}$ Brownian trajectories (for each $m$ we use new trajectories).

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.,

with parameter $0.2$ for the first three hidden layers and $\phi=\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.

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 (3.11) 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 (2.3). The only parameter that remains to be specified, is the number of trajectories used for the Monte Carlo estimator which is done in Algorithm D.1 and Specification D.2 below.

As a first calibration step, we calibrate the SABR model (i.e., (4.1) with $L\equiv 1$) to the synthetic market prices of the first maturity and fix the calibrated SABR parameters $\nu,\varrho$ 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 (3.11) and we compute the full gradient as specified in (3.13). 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,\ldots,4$, of the neural networks we apply Algorithm D.1 in the Appendix D.

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 4.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 4 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 3.

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 minutes. 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.

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

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

with $X_{j,m}$ defined as $X_{j}$ in (3.10) (see also (3.7)) but with synthetic market prices $m=1,\ldots,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.

Plots

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

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-Labordere (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.

Appendix A Variations of Stochastic Differential Equations

We follow here the excellent exposition of Protter (1990) to understand the dependence of solutions of stochastic differential equations on parameters, in particular when we aim to calculate derivatives with respect to parameters of neural networks.

the property $X^{\tau-}=Y^{\tau-}$ implies $F(X)^{\tau-}=F(Y)^{\tau-}$ for any stopping time $\tau$,

there exists an increasing process ${(K_{t})}_{t\geq 0}$ such that for $t\geq 0$

Functional Lipschitz assumptions are sufficient to obtain existence and uniqueness for general stochastic differential equations, see (Protter, 1990, Theorem V 7).

for $t\geq 0$ and $i=1,\ldots,n$. If $J$ is a semimartingale, then $Z$ is a semimartingale as well.

With an additional uniformity assumption on a sequence of stochastic differential equations with converging coefficients and initial data we obtain stability, see (Protter, 1990, Theorem V 15).

for $t\geq 0$ and $i=1,\ldots,n$. If $J^{\varepsilon}\to J^{0}$ in ucp, $F^{\varepsilon}(Z^{0})\to F^{0}(Z^{0})$ in ucp, then $Z^{\varepsilon}\to Z^{0}$ in ucp.

We shall apply these theorems to a local stochastic volatility model of the form

where $\theta\in\Theta$, $(W,\alpha)$ denotes some Brownian motion together with an adapted, càdlàg stochastic process $\alpha$ (all on a given stochastic basis) and $S_{0}>0$ is some real number.

We assume that for each $\theta\in\Theta$

is bounded, càdlàg in $t$ (for fixed $s>0$), and globally Lipschitz in $s$ with a Lipschitz constant independent of $t$ on compact intervals . In this case, the map

is functionally Lipschitz and therefore the above equation has a unique solution for all times $t$ and any $\theta$ by Theorem A.2. If, additionally,

where the $\sup$ is taken over some compact set, then we also have that the solutions $S(\theta)$ converge ucp to $S(\widehat{\theta})$, as $\theta\to\widehat{\theta}$ by Theorem A.3.

Appendix B Preliminaries on Deep Learning

We shall here briefly introduce two core concepts in deep learning, namely artificial neural networks and stochastic gradient descent. The latter is a widely used optimization method for solving maximization or minimization problems involving the first. In standard machine-learning terminology, the optimization procedure is usually referred to as “training”. We shall use both terminologies interchangeably.

We start with the definition of feed-forward neural networks. These are functions obtained by composing layers consisting of an affine map and a componentwise nonlinearity. They serve as universal approximation class which is stated in Theorem B.3. Moreover, derivatives of these functions can be efficiently expressed iteratively (see e.g. Hecht-Nielsen (1992)), which is a desirable feature from an optimization point of view.

is called a feed-forward neural network. Here the activation function $\phi$ is applied componentwise. $M-1$ denotes the number of hidden layers and $N_{1},\ldots,N_{M-1}$ denote the dimensions of the hidden layers and $N_{0}$ and $N_{M}$ the dimension of the input and output layers.

Unless otherwise stated, the activation functions $\phi$ used in this article are always assumed to be smooth, globally bounded with bounded first derivative.

The following version of the so-called universal approximation theorem is due to K. Hornik (Hornik, 1991). An earlier version was proved by G. Cybenko (Cybenko, 1989). To formulate the result, we denote the set of all feed-forward neural networks with activation function $\phi$, input dimension $N_{0}$ and output dimension $N_{M}$ by $\mathcal{NN}^{\phi}_{\infty,N_{0},N_{M}}$.

Suppose $\phi$ is bounded and nonconstant. Then the following statements hold:

We denote by $\mathcal{NN}_{N_{0},N_{M}}$ the set of all neural networks in $\mathcal{NN}^{\phi}_{\infty,N_{0},N_{M}}$ with a fixed architecture, i.e., a fixed number of hidden layers $M-1$, fixed input and output dimensions $N_{m}$ for each hidden layer $m\in\{1,\ldots,M-1\}$ and a fixed activation function $\phi$. This set can be described by

B.2. Stochastic Gradient Descent

In light of Theorem B.3, it is clear that neural networks can serve as function approximators. To implement this, the entries of the matrices $A_{m}$ and the vectors $b_{m}$ for $m\in\{1,\ldots,M\}$ are subject to optimization. If the unknown function can be expressed as the expected value of a stochastic objective function, one widely applied optimization method is stochastic gradient descent, which we shall review below.

Indeed, consider the following minimization problem

The classical method how to solve generic optimization problems for some differentiable objective function $f$ (not necessarily of the expected value form as in (B.1)) is to apply a gradient descent algorithm: starting with an initial guess $\theta^{(0)}$, one iteratively defines

for some learning rate $\eta_{k}$. Under suitable assumptions, $\theta^{(k)}$ converges for $k\to\infty$ to a local minimum of the function $f$.

In the deep learning context, stochastic gradient descent methods, going back to stochastic approximation algorithms proposed by Robbins and Monro (1951), are much more efficient. To apply this, it is crucial that the objective function $f$ is linear in the sampling probabilities. In other words, $f$ needs to be of the expected value form as in (B.1). In the simplest form of stochastic gradient descent, under the assumption that

the true gradient of $f$ is approximated by a gradient at a single sample $Q(\theta)(\omega)$ which reduces the computational cost considerably. In the updating step for the parameters $\theta$ as in (B.2), $f$ is then replaced by $Q(\theta)(\omega_{k})$, hence

The algorithm passes through all samples $\omega_{k}$ of the so-called training data set, possibly several times (specified by the number of epochs), and performs the update until an approximate minimum is reached.

A compromise between computing the true gradient of $f$ and the gradient at a single sample $Q(\theta)(\omega)$ is to compute the gradient of a subsample of size $N_{\text{batch}}$, called (mini)-batch, so that $Q(\theta^{(k)})(\omega_{k})$ used in the update (B.3) is replaced by

where $N$ is the size of the whole training data set. Any other unbiased estimators of $\nabla f(\theta)$ can of course also be applied in (B.3).

Appendix C Alternative Approaches for Minimizing the Calibration Functional

We consider here alternative algorithms for minimizing (3.8).

One alternative is stochastic compositional gradient descent as developed e.g. in Wang et al. (2017). Applied to our problem this algorithm (in its simplest form) works as follows: starting with an initial guess $\theta^{(0)}$, and $y_{j}^{(0)}$, $j=1,\ldots,J$ one iteratively defines

C.2. Estimators Compatible with Stochastic Gradient Descent

for some independent copy $\widetilde{Q}_{j}(\theta)$ of $Q_{j}(\theta)$, which is clearly of the expected value form required in (B.1). A Monte Carlo estimator of $f(\theta)$ is then constructed by

for independent draws $\omega_{1},\ldots,\omega_{N}$ (the same $N$ samples can be used for each strike $K_{j}$). Equivalently we have

for independent draws $\omega_{1},\ldots,\omega_{2N}$. The analog of (B.4) is then given by

for $k\in\{0,1,...,\lfloor N/N_{\text{batch}}\rfloor-1\}$.

Clearly we can now modify and improve the estimator by using again hedge control variates and replace $Q_{j}(\theta)$ by $X_{j}(\theta)$ as defined in (3.10).

Appendix D Algorithms

In this section, we present the calibration algorithm discussed above in form of pseudo code given in Algorithm D.1. Update rules for parameters in Algorithm D.1 are provided in Algorithm D.2. We further provide an implementation in form of a github repository, see https://github.com/wahido/neural_locVol.

In the subsequent pseudo code, the index $i$ stands for the maturities, $N$ for the number of samples used in the variance reduced Monte Carlo estimator as of (3.11) and $k$ for the updating step in the gradient descent:

k1k=k+1 Algorithm D.2. We update the parameters in Algorithm D.1 according to the following rules:

References