Neural Option Pricing for Rough Bergomi Model

Empirical studies of a wild range of assets volatility time-series show that the log-volatility in practice behaves similarly to fractional Brownian motion (FBM) (FBM with a Hurst index $H$ is the unique centered Gaussian process $W^{H}(t)$ with autocorrelation being $\mathbb{E}[W^{H}_{t}W^{H}_{s}]=\tfrac{1}{2}[|t|^{2H}+|s|^{2H}-|t-s|^{2H}]$) with a Hurst index $H\approx 0.1$ at any reasonable time scale [8]. Motivated by this empirical observation, several rough stochastic volatility models have been proposed, all of which are essentially based on fBM and involve the fractional kernel. The rough Bergomi (rBergomi) model [3] is one recent rough volatility model that has a remarkable capability of fitting both historical and implied volatilities.

The rBergomi model can be formulated as follows. Let $S_{t}$ be the price of underlying asset with the time horizon $[0,T]$ defined on a given filtered probability space $(\Omega,\mathcal{F},\{\mathcal{F}_{t}\}_{t\geq 0},\mathbb{Q})$, with $\mathbb{Q}$ being a risk-neutral martingale measure:

$$S_{t}=S_{0}\exp\left(-\frac{1}{2}\int_{0}^{t}V_{s}\mathrm{d}s+\int_{0}^{t}\sqrt{V_{s}}\mathrm{d}Z_{s}\right),\quad t\in[0,T].$$

(1)

Here, $V_{t}$ is the spot variance process satisfying

$$V_{t}=\xi_{0}(t)\exp\left(\eta\sqrt{2H}\int_{0}^{t}(t-s)^{H-\frac{1}{2}}\mathrm{d}W_{s}-\frac{\eta^{2}}{2}t^{2H}\right),\quad t\in[0,T].$$

(2)

$S_{0}>0$ denotes the initial value of the underlying asset, and the parameter $\eta$ is defined by

$\displaystyle\eta:=2\nu\sqrt{\frac{\Gamma(3/2-H)}{\Gamma(H+1/2)\Gamma(2-2H)}},$

with $\nu$ being the ratio between the increment of $\log V_{t}$ and the FBM over $(t,t+\Delta t)$ [3, Equation (2.1)].

The Hurst index $H\in(0,{1}/{2})$ reflects the regularity of the volatility $V_{t}$. $\xi_{0}(\cdot)$ is the so-called initial forward variance curve, defined by $\xi_{0}(t)=\mathbb{E}^{\mathbb{Q}}[V_{t}|\mathcal{F}_{0}]=\mathbb{E}[V_{t}]$ [3]. $Z_{t}$ is a standard Brownian motion:

$$Z_{t}:=\rho W_{t}+\sqrt{1-\rho^{2}}W^{\perp}_{t}.$$

(3)

Here, $\rho\in(-1,0)$ is the correlation parameter and $W^{\perp}_{t}$ is a standard Brownian motion independent of $W_{t}$. The price of a European option with payoff function $h(\cdot)$ and expiry $T$ is given by

$\displaystyle P_{0}=\mathbb{E}[\mathrm{e}^{-rT}h(S_{T})].$

(4)

Despite its significance from a modeling perspective, using a singular kernel in the rBergomi model leads to the loss of Markovian and semimartingale structure. In practice, simulating and pricing options using this model involves several challenges. The primary difficulty arises from the singularity of the fractional kernel

$\displaystyle G(t):=t^{H-\frac{1}{2}},$

(5)

at $t=0$, which impacts the simulation of the Volterra process given by:

$$I_{t}:=\sqrt{2H}\int_{0}^{t}G(t-s)\mathrm{d}W_{s}.$$

(6)

Note that this Volterra process is a Riemann-Liouville FBM or Lévy’s definition of FBM up to a multiplicative constant [18]. The deterministic nature of the kernel $G(\cdot)$ implies that $I_{t}$ is a centered, locally $(H-\epsilon)$-Hölder continuous Gaussian process with $\mathbb{E}[I_{t}^{2}]=t^{2H}$. Furthermore, a straightforward calculation shows

$\displaystyle\mathbb{E}[I_{t_{1}}I_{t_{2}}]=t_{1}^{2H}C\left(\frac{t_{2}}{t_{1}}\right),\quad\text{ for }t_{2}>t_{1},$

where for $x>1$, the univariate function $C(\cdot)$ is given by

$\displaystyle C(x):=2H\int_{0}^{1}(1-s)^{H-1/2}(x-s)^{H-1/2}\mathrm{d}s,$

indicating that $I_{t}$ is non-stationary.

Even though the rBergomi model enjoys remarkable calibration capability with the market data, there are many hidden unknown parameters even functions. For example, the Hurst exponent (in a general SDE) is regarded as an unknown function $H:=H(t)\colon\mathbb{R}_{+}\to(0,1)$ parameterized by a neural network in [21], which is learned using neural SDEs. According to [3], $\xi_{0}(t)$ can be any given initial forward variance swap curve consistent with the market price. In view of this fact and the high expressivity of neural SDEs [22, 23, 16, 14, 17, 10, 9], in this work, we propose to parameterize the initial forward variance curve $\xi_{0}(t)$ using a feedforward neural network:

$\displaystyle\xi_{0}(t)=\xi_{0}(t;\theta),$

(7)

where $\theta$ represents the weight parameter vector of the neural network. The training data are generated via a suitable numerical scheme for a given target initial forward variance curve $\xi_{0}(t)$.

Motivated by the implementation in [7], we adopt the Wasserstein metric as the loss function to train the weights $\theta$. Since the Wasserstein loss is convex and has a unique global minima, any SDE maybe learnt in the infinite data limit [13]. Remarkably, the attained Wasserstein-1 distance during the training is a natural upper bound for the discrepancy between the exact option price $P_{0}$ and the learned option price $P_{0}(\theta^{*})$ given by the neural SDE. This in particular implies that if the approximation of the underlying dynamics of stock price $S_{t}$ is accurate to some level, then the European option price $P_{0}$ with this stock as the underlying asset will be accurate to the same level. In this manner, the loss can realize two optimization goals simultaneously.

To generate the training data to learn $\theta$, we need to jointly simulate the Volterra process $I_{t}$ and the underlying Brownian motion $Z_{t}$ for the stock price. The non-stationarity feature and the joint simulation requirement make Cholesky factorization [3] the only available exact method. However, it has an $\mathcal{O}(n^{3})$ complexity for the Cholesky factorization, $\mathcal{O}(n^{2})$ complexity and $\mathcal{O}(n^{2})$ storage with $n$ being the total number of time steps. Clearly, the method is infeasible for large $n$. The most well known method to reduce the computational complexity and storage cost is the hybrid scheme [4]. The hybrid scheme approximates the kernel $G(\cdot)$ by a power function near zero and by a step function elsewhere, which yields an approximation combining the Wiener integrals of the power function with a Riemann sum. Since the fractional kernel $G(\cdot)$ is a power function, it is exact near zero. Its computational complexity is $\mathcal{O}(n\log n)$ and storage cost is $\mathcal{O}(n)$.

In this work, we propose an efficient modified summation of exponentials (mSOE) based method (16) to simulate (1)-(2) (to facilitate the training of the neural SDE). Similar approximations have already been considered, e.g. by Bayer and Breneis [2], and Abi Jaber and El Euch [1]. We further enhance the numerical performance by keeping the kernel exact near the singularity, which achieves high accuracy with the fewest number of summation terms. Numerical experiments show that the mSOE scheme considerably improve the prevalent hybrid scheme [4] in terms of accuracy, while having less computational and storage costs. Besides the rBergomi model, the proposed approach is applicable to a wide class of stochastic volatility models, especially the rough volatility models with completely monotone kernels.

In sum, the contributions of this work are three-fold. First, we derive an efficient modified sum-of-exponentials (mSOE) based method (16) to solve (1)-(2) even with very small Hurst parameter $H\in(0,1/2)$. Second, we propose to learn the forward variance curve by the neural SDE using the loss function based on the Wasserstein 1-distance, which can learn the underlying dynamics of the stock price as well as the European option price. Third and last, the mSOE scheme is further utilized to obtain the training data to train our proposed neural SDE model and serves as the solver for the neural SDE, which improves the efficiency significantly.

The remaining of the paper is organized as follows. In Section 2, we introduce an approximation of the singular kernel $G(\cdot)$ by the mSOE, and then describe two approaches for obtaining them. Based on the mSOE, we propose a numerical scheme for simulating the rBergomi model (1)-(2). We introduce in Section 3 the neural SDE model and describe the training of the model. We illustrate in Section 4 the numerical performance of the mSOE scheme. Moreover, several numerical experiments on the performance of our neural SDE model are shown for different target initial forward variance curves. Finally, we summarize the main findings in Section 5.

In this section, we develop an mSOE scheme for efficiently simulating the rBergomi model (1)-(2). Throughout, we consider an equidistant temporal grid $0=t_{0}<t_{1}<\cdots<t_{n}=T$ with a time stepping size $\tau:={T}/{n}$ and $t_{i}:=i\tau$.

This section is concerned with learning the forward variance curve $\xi_{0}(t;\theta)$ with $\theta$ being the weight parameters from the feedforward neural network (7) by Neural SDE. Without loss of generality, assuming $S_{0}=1$, we can rewrite (1)-(2) as

$$\left\{\begin{aligned} S_{t}&=1+\int_{0}^{t}S_{s}\exp(X_{s})\mathrm{d}Z_{s},\\ X_{t}&=\frac{1}{2}\log(V_{t}).\end{aligned}\right.$$

(17)

Upon parameterizing the forward variance curve by (7), the dynamics of the price process and variance process are also parameterized accordingly, which is given by the following neural SDE

$$\left\{\begin{aligned} S_{t}(\theta)&=1+\int_{0}^{t}S_{s}(\theta)\exp(X_{s}(\theta))\mathrm{d}Z_{s},\\ X_{t}(\theta)&=\frac{1}{2}\log(V_{t}(\theta)),\\ V_{t}(\theta)&=\xi_{0}(t;\theta)\mathcal{E}(\eta I(t)).\end{aligned}\right.$$

(18)

where $\mathcal{E}(\cdot)$ denotes the (Wick) stochastic exponential. Then we propose to use the Wasserstein-1 distance as the loss function to train this neural SDE. That is, the learned weight parameters $\theta^{*}$ are given by

$\displaystyle\theta^{*}:=\arg\min_{\theta}W_{1}\left(S_{T},S_{T}(\theta)\right).$

(19)

Next, we recall the Wasserstein distance. Let $(M,d)$ be a Radon space. The Wasserstein-$p$ distance for $p\in[1,\infty)$ between two probability measures $\mu$ and $\nu$ on $M$ with finite $p$-moments is defined by

$$W_{p}(\mu,\nu)=\left(\inf\limits_{\gamma\in\Gamma(\mu,\nu)}\mathbb{E}_{(x,y)\sim\gamma}d(x,y)^{p}\right)^{1/p},$$

where $\Gamma(\mu,\nu)$ is the set of all couplings of $\mu$ and $\nu$. A coupling $\gamma$ is a joint probability measure on $M\times M$ whose marginals are $\mu$ and $\nu$ on the first and second factors, respectively. If $\mu$ and $\nu$ are real-valued, then their Wasserstein distance can be simply computed by utilising their cumulative distribution functions [24]:

$$W_{p}(\mu,\nu)=\left(\int_{0}^{1}|F^{-1}(z)-G^{-1}(z)|^{p}\mathrm{d}z\right)^{1/p}.$$

In particular, when $p=1$, according to the Kantorovich-Robenstein duality [24], Wasserstein-1 distance can be represented by

$$W_{1}(\mu,\nu)=\sup\limits_{\mathrm{Lip}(f)\leq 1}\mathbb{E}_{x\sim\mu}[f(x)]-\mathbb{E}_{y\sim\nu}[f(y)],$$

(20)

where $\mathrm{Lip}(f)$ denotes the Lipschitz constant of $f$.

In the experiment, we work with the empirical distributions on $\mathbb{R}$. If $\xi$ and $\eta$ are two empirical measures with $m$ samples $X_{1},\cdots,X_{m}$ and $Y_{1},\cdots,Y_{m}$, then the Wasserstein-$p$ distance is given by

$$W_{p}(\xi,\eta)=\left(\frac{1}{m}\sum_{i=1}^{m}|{X_{(i)}-Y_{(i)}}|^{p}\right)^{1/p},$$

(21)

where $X_{(i)}$ is the $i$th smallest value among the $m$ samples. With $S_{t}(\theta^{*})$ being the price process of the underlying asset after training, the price $P_{0}(\theta^{*})$ of a European option with payoff function $h(x)$ and expiry $T$ is given by

$\displaystyle P_{0}(\theta^{*})=\mathbb{E}[\mathrm{e}^{-rT}h(S_{T}(\theta^{*}))].$

(22)

Since the payoff function $h(\cdot)$ is clearly Lipschitz-1 continuous, according to (20), Wasserstein-1 distance is a natural upper bound for the pricing error of the rBergomi model. Thus, the choice of the training loss is highly desirable.

In this section, we illustrate the performance of the proposed mSOE scheme (16) for simulating (1)-(2), and demonstrate the learning of the forward variance curve (7) using neural SDE (18).