Decomposition formula for rough Volterra stochastic volatility models

TL;DR

This paper introduces a new decomposition formula for pricing European options under rough fractional stochastic volatility models, providing an efficient approximation method that enhances calibration speed compared to traditional Monte Carlo simulations.

Contribution

It offers a novel decomposition formula and explicit approximation for option pricing in rough volatility models, improving computational efficiency and calibration speed.

Findings

The approximation formula performs well for the rBergomi model.

The hybrid calibration scheme accelerates market calibration.

Numerical tests show good accuracy of the approximation.

Abstract

The research presented in this article provides an alternative option pricing approach for a class of rough fractional stochastic volatility models. These models are increasingly popular between academics and practitioners due to their surprising consistency with financial markets. However, they bring several challenges alongside. Most noticeably, even simple non-linear financial derivatives as vanilla European options are typically priced by means of Monte-Carlo (MC) simulations which are more computationally demanding than similar MC schemes for standard stochastic volatility models. In this paper, we provide a proof of the prediction law for general Gaussian Volterra processes. The prediction law is then utilized to obtain an adapted projection of the future squared volatility -- a cornerstone of the proposed pricing approximation. Firstly, a decomposition formula for European option prices under general Volterra volatility models is introduced. Then we focus on particular models with rough fractional volatility and we derive an explicit semi-closed approximation formula. Numerical properties of the approximation for a popular model -- the rBergomi model -- are studied and we propose a hybrid calibration scheme which combines the approximation formula alongside MC simulations. This scheme can significantly speed up the calibration to financial markets as illustrated on a set of AAPL options.

Tables

Table 1. Table 1 : Model / data settings for sensitivity analyses.

Model params	Values
$\xi$	$\{10\%,50\%,100\%\}$
${\sigma }_0$	8%
$\rho$	-20%
$H$	0.1
Data set specifics	Values
Moneyness	70% – 130% with 5% step
Time to maturity	{1M, 3M, 6M, 1Y}

Table 2. Table 2: Data on AAPL options used in calibration trials

	Data # 1	Data # 2	Data # 3	Data # 4
Date (all EOD)	1-Apr-2015	15-Apr-2015	1-May-2015	15-May-2015
Moneyness range	34%–157%	45%–154%	31%–151%	33%–151%
Time to maturity [Yrs]	0.12-1.81	0.08–1.77	0.04–1.73	0.02–1.69
Total nb. of contracts	113	158	201	194

Equations

$$\mathrm{d}{S}_t=r{S}_t\mathrm{d}t+{\sigma }_t{S}_t\left(\rho \mathrm{d}{W}_t+\sqrt{1-{\rho }^2}\mathrm{d}{\tilde{W}}_t\right),$$

$$\mathrm{d}{X}_t=\left(r-\frac{1}{2}{\sigma }_t^2\right)\mathrm{d}t+{\sigma }_t\left(\rho \mathrm{d}{W}_t+\sqrt{1-{\rho }^2}\mathrm{d}{\tilde{W}}_t\right).$$

$$BS\left(t,x,y\right)=e^x\Phi (d_{+})-K{e}^{-r\tau }\Phi (d_{-}),$$

$$d_{\pm }(y)=\frac{x-\ln K+(r\pm \frac{y^2}{2})\tau }{y\sqrt{\tau }}.$$

$$V_t=e^{-r\tau }{\mathbb{E}}_t[(e^{X_T}-K{)}^{+}]$$

$${\mathcal{L}}_y:={\partial }_t+\frac{1}{2}{y}^2{\partial }_x^2+\left(r-\frac{1}{2}{y}^2\right){\partial }_x-r$$

$$\Gamma BS(t,x,y)$$

$$\Lambda \Gamma BS(t,x,y)$$

$${\Gamma }^{2}BS(t,x,y)$$

$$R_t:=\frac{1}{8}{\mathbb{E}}_t\left[{\int }_t^T\mathrm{d}{⟨M,M⟩}_u\right]$$

$$U_t:=\frac{\rho }{2}{\mathbb{E}}_t\left[{\int }_t^T{\sigma }_u\mathrm{d}{⟨M,W⟩}_u\right],$$

$$M_t:={\int }_0^T{\mathbb{E}}_t\left[{\sigma }_s^2\right]\mathrm{d}s.$$

$$V_t={\mathbb{E}}_t[BS(t,X_t,{\bar{\sigma }}_t)]$$

$${\bar{\sigma }}_t^2:=\frac{1}{T-t}{\int }_t^T{\sigma }_u^2\mathrm{d}u.$$

$$a_t^2:={\int }_t^T{\mathbb{E}}_t[{\sigma }_u^2]\mathrm{d}u$$

$$v_t^2:={\mathbb{E}}_t({\bar{\sigma }}_t^2)=\frac{a_t^2}{T-t}$$

$$\mathrm{d}{v}_t^2=\frac{1}{T-t}\left[\mathrm{d}{M}_t+\left(v_t^2-{\sigma }_t^2\right)\mathrm{d}t\right].$$

$$\left|\mathbb{E}\left({\Lambda }^n\Gamma BS\left(u,X_u,v_u\right)|{\mathcal{G}}_t\right)\right|\le C(a_u^2{)}^{-\frac{1}{2}(n+1)},$$

$${\mathbb{E}}_t\left[e^{-r(T-t)}A(T,X_T,v_T^2)B_T\right]=A(t,X_t,v_t^2)B_t$$

$$V_t$$

$$+\frac{\rho }{2}{\mathbb{E}}_t\left[{\int }_t^T{e}^{-r(u-t)}\Lambda \Gamma BS(u,X_u,v_u){\sigma }_u\mathrm{d}{⟨W_{to2.91666pt\text{scalebox}1.25\cdot },M_{to2.91666pt\text{scalebox}1.25\cdot }⟩}_u\right]$$

$$+\frac{1}{8}{\mathbb{E}}_t\left[{\int }_t^T{e}^{-r(u-t)}{\Gamma }^{2}BS(u,X_u,v_u)\mathrm{d}{⟨M_{to2.91666pt\text{scalebox}1.25\cdot },M_{to2.91666pt\text{scalebox}1.25\cdot }⟩}_u\right]$$

$$=BS(t,X_t,v_t)+(I)+(II)$$

$$A(t,X_t,v_t^2)=\Lambda \Gamma BS(t,X_t,v_t)$$

$$B_t=U_t=\frac{\rho }{2}{\mathbb{E}}_t\left[{\int }_t^T{\sigma }_u\mathrm{d}{⟨W_{to2.91666pt\text{scalebox}1.25\cdot },M_{to2.91666pt\text{scalebox}1.25\cdot }⟩}_u\right]$$

$$A(t,X_t,v_t^2)={\Gamma }^{2}BS(t,X_t,v_t)$$

$$B_t=R_t=\frac{1}{8}{\mathbb{E}}_t\left[{\int }_t^T\mathrm{d}{⟨M_{to2.91666pt\text{scalebox}1.25\cdot },M_{to2.91666pt\text{scalebox}1.25\cdot }⟩}_u\right]$$

$$V_t$$

$$+\Lambda \Gamma BS(t,X_t,v_t)U_t$$

$$+{\Gamma }^{2}BS(t,X_t,v_t)R_t$$

$$+{ϵ}_t.$$

$${\sigma }_t:=g(t,Y_t),t\ge 0,$$

$$Y_t={\int }_0^{t}K(t,s)\mathrm{d}{W}_s,$$

$$\underset{0}{\overset{t}{\int }}{K}^2(t,s)\mathrm{d}s<\infty ,$$

Full text

In this paper, we provide a proof of the prediction law for general Gaussian Volterra processes. The prediction law is then utilized to obtain an adapted projection of the future squared volatility – a cornerstone of the proposed pricing approximation. Firstly, a decomposition formula for European option prices under general Volterra volatility models is introduced. Then we focus on particular models with rough fractional volatility and we derive an explicit semi-closed approximation formula. Numerical properties of the approximation for a popular model – the rBergomi model – are studied and we propose a hybrid calibration scheme which combines the approximation formula alongside MC simulations. This scheme can significantly speed up the calibration to financial markets as illustrated on a set of AAPL options.

Received 1 August 2019

Keywords: Volterra stochastic volatility; rough volatility; rough Bergomi model; option pricing; decomposition formula

MSC classification: 60G22; 91G20; 91G60

JEL classification: G12; C58; C63

1 Introduction

It is well known that the main issue of the Black-Scholes model lies in its assumptions about volatility of the modelled asset. Opposed to the model assumptions, the realized volatility time series tend to cluster, depend on the spot asset level and certainly they do not take a constant value within a reasonable time-frame (see e.g. Cont (2001)).

To deal with the aforementioned inconsistencies, stochastic volatility (SV) models have been proposed originally by Hull and White (1987) and later e.g. by Heston (1993). These models do not only assume that the asset price follows a specific stochastic process, but also that the instantaneous volatility of asset returns is of random nature as well. Especially, the latter approach by Heston became popular in the eyes of both practitioners and academics. Several modifications of this model have been proposed over the last 20 years: models with jump-diffusion dynamics (Bates 1996; Duffie, Pan, and Singleton 2000), with time-dependent parameters (Mikhailov and Nögel 2003; Elices 2008; Benhamou, Gobet, and Miri 2010), with fractionally scaled volatility (Comte and Renault 1998; Alòs, León, and Vives 2007; El Euch and Rosenbaum 2019) and models with several aspects combined (Pospíšil and Sobotka 2016; Baustian, Mrázek, Pospíšil, and Sobotka 2017).

The original pricing approach of Heston (1993) was several times revisited, e.g by Lewis (2000), Attari (2004) and Albrecher, Mayer, Schoutens, and Tistaert (2007) with focus on semi-closed form Fourier transform solutions by Kahl and Jäckel (2006), Alfonsi (2010) with respect to Monte-Carlo simulation techniques and last but not least by Alòs (2012) who introduced an analytical approach to option pricing approximation. This approach improves on the techniques introduced by Hull and White (1987) and shows how an adaptive projection of future volatility can be used to price European options under Heston (1993) model. Many other papers generalized this idea, see e.g. Alòs, de Santiago, and Vives (2015); Merino and Vives (2015), or recently also a paper by Merino, Pospíšil, Sobotka, and Vives (2018). In this article, we revisit the latter approach and we come up with the approximation technique for SV models with volatility process driven by the fractional Brownian motion. This includes exponential rough fractional volatility models introduced by Gatheral, Jaisson, and Rosenbaum (2014, 2018).

Although many SV models have been proposed since the original Hull and White (1987) model, it seems that none of them can be considered as the universal best market practice approach. Several models might perform well for calibration to complex volatility surfaces, but can suffer from over-fitting or they might not be robust in the sense described by Pospíšil, Sobotka, and Ziegler (2018). Also a model with a good fit to implied volatility surface might not be in-line with the observed time-series properties.

Pioneers of the fractional SV models – Comte and Renault (1998), see also Comte, Coutin, and Renault (2012) – assumed the so-called Hurst parameter 111Named after a hydrologist Harold Edwin Hurst, for more information on fractional Brownian motion, see Section 4.3 and e.g. the article by Mandelbrot and Van Ness (1968) ranged within $H\in(1/2,1)$ which implies that the spot variance evolution is represented by a persistent process, i.e. it would have a long-memory property. In Alòs, León, and Vives (2007), a mean reverting fractional stochastic volatility model with $H\in(0,1)$ was presented. Gatheral, Jaisson, and Rosenbaum (2018) and Bayer, Friz, and Gatheral (2016) came up with a more detailed analysis of rough fractional volatility models which should be consistent with market option prices (Bayer, Friz, and Gatheral 2016), with realized volatility time series and also they could provide superior volatility prediction results to several other models (Bennedsen, Lunde, and Pakkanen 2017). An approach considering a two factor fractional volatility model, combining a rough term ($H<\frac{1}{2}$) and a persistent one ($H>\frac{1}{2}$), was presented in Funahashi and Kijima (2017). Recently, also an approximation for target-volatility options under log-normal fractional SABR model was studied by Alòs, Chatterjee, Tudor, and Wang (2019) who use the Malliavin caluculus techniques to derive the decomposition formula. In parallel, short-term at-the-money asymptotics for a class of stochastic volatility models were studied by El Euch, Fukasawa, Gatheral, and Rosenbaum (2019) who use the Edgeworth expansion of the density of an asset price.

A typical problem of rough fractional models lies in their tractability – only cumbersome simulation techniques seem to be available for vanilla European options at the time of writing this article. This serves as a motivation to develop a pricing approximation based on the works of Alòs (2012), whose approach was further generalised by Merino and Vives (2015) and Merino, Pospíšil, Sobotka, and Vives (2018). Without using the Malliavin calculus we derive a decomposition and approximation formula for European option prices under general Volterra volatility models and in particular under a class of models with rough fractional volatility. Obtained results hence give a better understanding of the whole surface of prices than results focusing only on at-the-money volatility skew.

The structure of the paper is the following. Section 2 is devoted to preliminaries. In Section 3 we present a generic decomposition formula of the vanilla European option fair value. In Setion 4 we present Volterra volatility models. In particular in Section 4.1 we prove the prediction law for general Gaussian Volterra processes, in Section 4.2 we obtain the decomposition formula for the exponential Volterra volatility model including the error bound for the approximation formula and in Section 4.3 we provide particular results for exponential fracional volatility models (including the standard Brownian motion case). In Section 5, we provide a numerical comparison of approximated prices against Monte Carlo simulations. All obtained results are concluded in Section 6.

2 Preliminaries and notation

Let $S=(S_{t},t\in[0,T])$ be a strictly positive asset price process under a market chosen risk neutral probability measure $\mathcal{P}$ that follows the model:

[TABLE]

where $S_{0}$ is the current price, $r\geq 0$ is the interest rate, $W_{t}$ and $\tilde{W}_{t}$ are independent standard Wiener processes defined on a probability space $(\Omega,\mathcal{F},\mathcal{P})$ and $\rho\in\left(-1,1\right)$. In the following, we will denote by $\mathcal{F}^{W}$ and $\mathcal{F}^{\tilde{W}}$ the filtrations generated by $W$ and $\tilde{W}$ respectively. Moreover, we define $\mathcal{F}:=\mathcal{F}^{W}\vee\mathcal{F}^{\tilde{W}}$. The volatility process $\sigma_{t}$ is a square-integrable process assumed to be adapted to the filtration generated by $W$ and its trajectories are assumed to be a.s. càdlàg and stricly positive a.e.. Note that $\rho$ is the correlation between the price and the volatility processes.

Without any loss of generality, it will be convenient in the following sections to make the change of variable $X_{t}=\log S_{t},t\in[0,T]$, and write

[TABLE]

Recall that $Z:=\rho W+\sqrt{1-\rho^{2}}\tilde{W}$ is a standard Wiener process.

The following notation will be used throughout the paper:

• •

We will denote by $BS(t,x,y)$ the price of a plain vanilla European call option under the classical Black-Scholes model with constant volatility $y$, current log stock price $x$, time to maturity $\tau=T-t$, strike price $K$ and interest rate $r$. In this case,

[TABLE]

where $\Phi(\cdot)$ denotes the cumulative distribution function of the standard normal law and

[TABLE]

• •

In our setting, the call option price is given by

[TABLE]

where ${\mathbb{E}}_{t}$ is the conditional expectation respect to the $\sigma-$algebra $\mathcal{F}_{t}$.

• •

Recall that from the Feynman-Kac formula for the model (2), the operator

[TABLE]

satisfies ${\mathcal{L}}_{y}BS(t,x,y)=0$.

• •

We define the operators $\Lambda:=\partial_{x}$, $\Lambda^{n}:=\partial^{n}_{x}$, $\Gamma:=\left(\partial^{2}_{x}-\partial_{x}\right)$ and $\Gamma^{2}=\Gamma\circ\Gamma.$ In particular, for the Black-Scholes formula, using straightforward calculations, we get:

[TABLE]

• •

We define

[TABLE]

and

[TABLE]

where $\!\left\langle\cdot,\cdot\right\rangle$ denotes the quadratic covariation process and $M$ is the $\mathcal{F}-$martingale defined by

[TABLE]

3 A generic decomposition formula

In this section, we provide an insight on a generic decomposition formula based on the work of Alòs (2012), Merino and Vives (2015) and Merino, Pospíšil, Sobotka, and Vives (2018). In particular, we recover the results for a generic stochastic volatility model presented in Merino, Pospíšil, Sobotka, and Vives (2018).

It is well known that if the stochastic volatility process is independent of the price process, the pricing formula of a vanilla European call option is given by

[TABLE]

where ${\bar{\sigma}}^{2}_{t}$ is the so called average future variance that is defined by

[TABLE]

Naturally, ${\bar{\sigma}}_{t}$ is called the average future volatility.

We consider the adapted projection of the future variance

[TABLE]

and the average future variance as

[TABLE]

to obtain a decomposition of $V_{t}$ in terms of $v_{t}.$ This idea switches an anticipative problem related to the anticipative process ${\bar{\sigma}}_{t}$ into a non-anticipative one related to the adapted process $v_{t}$.

Taking into account $M$ defined in (4), we can write

[TABLE]

In this paper, we will utilize the following lemma which is proved in Alòs (2012), p. 406.

Lemma 3.1**.**

Let $0\leq t\leq u\leq T$ and $\mathcal{G}_{t}:=\mathcal{F}_{t}\vee\mathcal{F}_{T}^{W}.$ For every $n\geq 0,$ there exists $C=C(n)$ such that

[TABLE]

where $a_{t}^{2}$ is defined by (5).

Remark 3.2**.**

It is easy to see that the previous Lemma holds for put options and for several other non-path dependent options as well (e.g. Gap options).

Now we use a generic decomposition formula proved in Merino, Pospíšil, Sobotka, and Vives (2018).

Theorem 3.3** (Generic decomposition formula).**

Let $B_{t}$ be a continuous semimartingale with respect to the filtration $\mathcal{F}^{W}$, let $A(t,x,y)$ be a $C^{1,2,2}([0,T]\times[0,\infty)\times[0,\infty))$ function and let $v^{2}_{t}$ and $M_{t}$ be defined as above. Then we are able to formulate the expectation of $e^{-rT}A(T,X_{T},v^{2}_{T})B_{T}$ in the following way:

[TABLE]

Proof.

See the proof of Theorem 3.1 in Merino, Pospíšil, Sobotka, and Vives (2018). ∎

Using the previous decomposition formula, we find that

Corollary 3.4** (BS decompostion formula).**

Under assumptions of the Theorem 3.3, we can obtain a decomposition of European option price $V_{t}$ as:

[TABLE]

Proof.

Using Theorem 3.3 with $A(t,X_{t},v^{2}_{t})=BS(t,X_{t},v_{t})$ and $B\equiv 1$, the proof follows in a straightforward way. ∎

The terms $(I)$ and $(II)$ are not easy to evaluate. Therefore, it becomes important to find simpler approximations to $(I)$ and $(II)$ and estimate the error terms. In order to find these approximations, we are going to apply Theorem 3.3 to find a decomposition formula for the terms $(I)$ and $(II)$. Using

[TABLE]

and

[TABLE]

a decomposition of the term $(I)$ can be found, and using

[TABLE]

and

[TABLE]

a decomposition of the term $(II)$ is obtained.

After that process, we can approximate the price of a call option by

[TABLE]

where $\epsilon_{t}$ denotes error terms. Terms of $\epsilon_{t}$ under a general setting for $\sigma_{t}$ are provided in Appendix A. We note that the error term will depend on the assumed volatility dynamics.

4 Volterra volatility models

4.1 General Volterra volatility model

In this section, we apply the generic decomposition formula to model (2) with general Volterra volatility process defined as

[TABLE]

where $g:[0,+\infty)\times\mathbb{R}\mapsto[0,+\infty)$ is a deterministic function such that $\sigma_{t}$ belongs to $L^{1}(\Omega\times[0,+\infty))$ and $Y=(Y_{t},t\geq 0)$ is the Gaussian Volterra process

[TABLE]

where $K(t,s)$ is a kernel such that for all $t>0$

[TABLE]

and

[TABLE]

Let

[TABLE]

denote the autocovariance function of process $Y_{t}$ and

[TABLE]

be the variance function (i.e. the second moment).

Extending the Theorem 3.1 in Sottinen and Viitasaari (2017) enables us to rephrase the adapted projection of the future squared volatility.

Theorem 4.1** (Prediction law for Gaussian Volterra processes).**

Let $(Y_{t},t\geq 0)$ be the Gaussian Volterra process (7) satisfying assumptions (A1) and (A2). Then, the conditional process $(Y_{u}|\mathcal{F}_{t},0\leq t\leq u)$ is Gaussian with $\mathcal{F}_{u}$-measurable mean function

[TABLE]

for $u_{1},u_{2}\geq t$.

Proof.

Let $0\leq t\leq u$. Then

[TABLE]

∎

In the upcoming sections, we will denote $\hat{r}(u|t):=\hat{r}(u,u|t)$.

Under the general volatility process (6), we have

[TABLE]

and the martingale

[TABLE]

In the upcoming lemma, we express the conditional expectation of the future squared volatility in terms of the mean function $\hat{m}_{t}(u)$.

Let us denote:

[TABLE]

Lemma 4.2** (Auxiliary terms in the decomposition formula for the general volatility model).**

Let $0\leq t\leq u$ and $F(t,\hat{m}_{t}(u))=\mathbb{E}_{t}\left[g^{2}(u,Y_{u})\right]$, then

[TABLE]

Proof.

Let $0\leq t\leq u$ and

[TABLE]

Theorem 4.1 implies that

[TABLE]

where

[TABLE]

is a Gaussian density function with stochastic mean $\hat{m}_{t}(u)$ and deterministic variance $\hat{r}(u|t)$. To calculate the quadratic variation, we note that

[TABLE]

where

[TABLE]

Since

[TABLE]

we retrieve the following expression by Itô’s formula,

[TABLE]

and consequently,

[TABLE]

Due to:

[TABLE]

we obtain

[TABLE]

More generally, we have

[TABLE]

and consequently,

[TABLE]

Let $F=F(t,m)$ be a $\mathcal{C}^{1,2}$-function of time $t$ and “spot” $m=\hat{m}_{t}(u)$ of the prediction martingale. Because of the filtrations are the same, we have, in general,

[TABLE]

Now we set

[TABLE]

and applying the Itô formula, we obtain (12). Moreover,

[TABLE]

and

[TABLE]

∎

4.2 Exponential Volterra volatility model

Assume now that $X_{t}$ is the log-price process (2) with $\sigma_{t}$ being the exponential Volterra volatility process

[TABLE]

where $(Y_{t},t\geq 0)$ is the Gaussian Volterra process (7) satisfying assumptions (A1) and (A2), $r(t)$ is its autocovariance function (9), and $\sigma_{0}>0$, $\xi>0$ and \leavevmode\color[rgb]{0,0,0}\alpha\in[0,1] are model parameters.

Lemma 4.3** (Auxiliary terms in the decomposition formula for the exponential Volterra volatility model).**

Let $\sigma_{t}$ be as in (19) and $0\leq t\leq u$. Then

[TABLE]

Proof.

Let $\hat{\varphi}_{t}(u,z)$ be given by (16). Then

[TABLE]

Remaining formulas follow accordingly. ∎

Remark 4.4**.**

Using that $F(t,\hat{m}_{t}(u))=\mathbb{E}_{t}\left[\sigma^{2}_{u}\right]$, it is straightforward to see that

[TABLE]

Lemma 4.5**.**

Let $\sigma_{t}$ be as in (19) and $0\leq t\leq u$. Then, we can re-write $F(t,\hat{m}_{t}(u))$ as

[TABLE]

Moreover, we also have the following equalities

[TABLE]

and

[TABLE]

Proof.

The calculations to obtain these statements are straightforward. ∎

Proposition 4.6** (Terms in the approximation formula for the exponential Volterra volatility model).**

Let $\sigma_{t}$ be as in (19) and $0\leq t\leq u$. Then

[TABLE]

Proof.

We have that

[TABLE]

Similarly, we have that

[TABLE]

∎

For the exponential Volterra volatility model we can determine an upper error bound for the price approximation in the following way.

Theorem 4.7** (Upper error bound for the exponential Volterra volatility model).**

Let $X_{t}$ be a log-price process (2) with $\sigma_{t}$ being the exponential Volterra volatility process (19). Then we can express the call option fair value $V_{t}$ using the processes $R_{t},U_{t}$ from Proposition 4.6. In particular,

[TABLE]

where $\epsilon_{t}$ are error terms of order $O\left(\left(\xi^{2}+\rho\xi\right)^{2}\right)$.

Proof.

Note that using (24) we have that

[TABLE]

Applying the Jensen’s inequality to (27), we can see that

[TABLE]

where the exponent

[TABLE]

is a Gaussian process. Therefore $\frac{1}{a^{2}_{t}}$ has finite moments of all orders.

Using the error terms specified in the Appendix A and Lemma 3.1, we find the following decompositions for each term

[TABLE]

and

[TABLE]

Since all previous conditional expectations are obviously finite, the error order is $O\left(\left(\xi^{2}+\rho\xi\right)^{2}\right)$.

∎

Further, we express differentials with respect to the $n^{\text{th}}$-power of the exponential Volterra volatility process when $Y_{t}$ is a semimartingale.

Lemma 4.8**.**

Let $\sigma_{t}$ be as in (19) and $Y_{t}$ a semimartingale. Let $n\geq 1$, we have that

[TABLE]

Proof.

The formula is an immediate consequence of the Itô formula. ∎

Lemma 4.9**.**

Let $\sigma_{t}$ be as in (19) and $Y_{t}$ is a semimartingale. We can calculate $\,\mathrm{d}U_{t}$ and $\,\mathrm{d}R_{t}$. In order to simplify the notation, we define the two following functions

[TABLE]

Proof.

[TABLE]

∎

Remark 4.10**.**

Note that the Theorem 4.7 requires that the conditional expectations of (4.2) and (4.2) be finite. In the case where $\sigma_{t}$ is as in (19) and $Y_{t}$ is a semimartingale, it easy to see that these conditions are met using Lemma 4.9.

4.3 Exponential fractional volatility model

Let us now focus on a very important example of Gaussian Volterra processes, the fractional Brownian motion (fBm), which is a process with a Hurst parameter $H\in(0,1)$, with covariance function

[TABLE]

and, in particular, with variance

[TABLE]

One of the first applications of fractional Brownian motion is credited to Hurst (1951) who modelled the long term storage capacity of reservoirs along the Nile river. However, the idea of this concept goes back to Kolmogorov (1940), who studied Wiener spirals and some other curves in Hilbert spaces. Later, Lévy (1953) used the Riemann–Liouville fractional integral to define the process as

[TABLE]

where $H$ may be any positive number. This type of integral turned out to be ill-suited to applications of fractional Brownian motion because of its over-emphasis on the origin for many applications. In their highly cited work, Mandelbrot and Van Ness (1968) introduced the Weyl’s representation of the fractional Brownian motion:

[TABLE]

where

[TABLE]

and $W_{t}$ is the standard Wiener process. Nowadays, the most widely used representation of fBm is the one by Molchan and Golosov (1969)

[TABLE]

To understand the connection between Molchan-Golosov and Mandelbrot-Van Ness representations of fBm we refer readers to the paper by Jost (2008).

Despite of the above mentioned arguments, Alòs, Mazet, and Nualart (2000) proposed to consider a process $\hat{B}_{t}=\int_{0}^{t}(t-s)^{H-1/2}\,\mathrm{d}W_{s}$ instead of $B^{H}_{t}$ in fractional stochastic calculus, since $Z_{t}$ has absolutely continuous trajectories. Since $B^{H}_{t}$ is not a semimartingale, the process $\hat{B}_{t}=\Gamma(H+1/2)B^{H}_{t}-Z_{t}$ is also not a semimartingale. Later on, Thao (2006) introduced the so called approximate fractional Brownian motion process as

[TABLE]

and showed that for every $\varepsilon>0$ the process $\hat{B}^{\varepsilon}_{t}$ is a semimartingale and it converges to $\hat{B}_{t}$ in $L^{2}(\Omega)$ when $\varepsilon$ tends to zero. This convergence is uniform with respect to $t\in[0,T]$ (Thao 2006, Theorem 2.1).

Let us now consider the exponential fractional volatility process

[TABLE]

where $(B^{H}_{t},t\geq 0)$ is one of the above mentioned representations of fBm. We are especially interested in the "rough" models, i.e. when $H<1/2$. In this case, we call the model (1) with volatility process (40) the rough fractional stochastic volatility model (\leavevmode\color[rgb]{0,0,0}\alphaRFSV). Note that if \leavevmode\color[rgb]{0,0,0}\alpha=1, we get the rBergomi model (Bayer, Friz, and Gatheral 2016; Gatheral, Jaisson, and Rosenbaum 2018), if \leavevmode\color[rgb]{0,0,0}\alpha=0, we get the original exponential fractional volatility model. Values of \leavevmode\color[rgb]{0,0,0}\alpha between 0 and 1 give us a new degree of freedom that can be viewed as a weight between these two models.

Example 4.11** (Volatility driven by the approximate fractional Brownian motion).**

Let us consider model (2) with volatility process

[TABLE]

Example 4.12** (Volatility driven by the standard Wiener process).**

If in the previous Example 4.11 we take $H=1/2$ and $\varepsilon=0$, we get model (2) with exponential Wiener volatility process

[TABLE]

Remark 4.13**.**

We can do a Taylor expansion of $U_{0}$ and $R_{0}$ to understand better their dependencies, doing that we obtain

[TABLE]

Theorem 4.14** (Decomposition formula for exponential Wiener volatility model).**

Let $X_{t}$ be the log-price process (2) with $\sigma_{t}$ being the exponential Wiener volatility process defined in (41). Assuming without any loss of generality that the options starts at time 0, then we can express the call option fair value $V_{0}$ using the processes $U_{0},R_{0}$ from (44) and (46) respectively. In particular,

[TABLE]

where $\epsilon$ denotes error terms and for \leavevmode\color[rgb]{0,0,0}\alpha\geq 0, $\left|\epsilon\right|$ is at most of the order C\xi(\sqrt{T}+\rho\xi^{2})T^{3/2}\Pi(\leavevmode\color[rgb]{0,0,0}\alpha,T,\xi,\rho). The exact bound is given in Appendix B.

Proof.

The detailed proof is given in Appendix B, where we also examine the order of magnitude by the first Taylor term of the integrals. ∎

Remark 4.15**.**

It is worth to mention that the order of the error bound from Theorem 4.14 is better than the general estimate from Theorem 4.7, where the time dependency is not considered. To get finer estimates also for the exponential fractional model (case $H\neq 1/2$), a proof similar to the one in Appendix B would have to be performed with more complicated but still tractable calculations.

Example 4.16** (Volatility driven by the standard fractional Brownian motion).**

Let us consider a model with volatility process

[TABLE]

where $B^{H}_{t}$ is the standard fractional Brownian motion as defined in (38), i.e. with the Molchan-Golosov kernel (39). Then, the formulas for $U_{0}$ and $R_{0}$ are given in Proposition 4.6 with particular kernel (39), autocovariance function (36) and $\hat{r}(t,s|u)$ as in Theorem 4.1. In this case, we do not give the formulas for $U_{0}$ and $R_{0}$ after substituting the Molchan-Golosov kernel, since these formulas are too long. However, it is worth to mention that the formulas are explicit and numerical evaluation requires only the computation of some multiple Gaussian integrals.

5 Numerical comparison of approximation formula

In this section, we focus on numerical aspects of the introduced approximation formula. We detail on its numerical implementation and a comparison with the Monte Carlo (MC) simulation framework introduced by Bennedsen, Lunde, and Pakkanen (2017) will be provided.

In the second part of this section, we also introduce two interesting outcome analysis for rBergomi model. In particular, we show how the model can be efficiently calibrated using the approximation formula to short maturity smiles. We remark that classical SV models (e.g. Heston model) might fail to fit the short term smiles, unless they exploit high volatility of volatility levels for which they would be typically inconsistent with the long term skew of the volatility surface.

In what follows, we will inspect the approximation quality for rBergomi model and time to maturity / volatility of volatility $\xi$ scenarios. Based on the nature of error terms (see Appendix A) those two factors should play prominent role when it comes to approximation quality.

5.1 On implementation of the approximation formula

We note that for the models studied in this paper, we have obtained either a semi-closed form or analytical formula for standard Wiener case (H=0.5). Moreover, for the class of exponential fractional models – represented by the \leavevmode\color[rgb]{0,0,0}\alphaRFSV model – we only need to numerically evaluate multiple integrals in $R_{0}$ and $U_{0}$.

In our case, this was done using a trapezoid quadrature routine – not necessarily the most efficient approach, but easy to implement. We used a discretisation222Typically we used from 1000 up to 27000 points for 3D integrals. of integrands such that the numerical error doesn’t affect the results in a significant way. I.e. to be lower than standard MC errors when comparing to simulated prices or lower than the expected approximation error.

For benchmarking we use a first-order hybrid MC scheme introduced by Bennedsen, Lunde, and Pakkanen (2017) alongside 50000 MC sample paths. Similarly to the implementation of the approximation formula, we remark that this scheme could be also improved as described in McCrickerd and Pakkanen (2018).

5.2 Sensitivity analysis for rBergomi (\leavevmode\color[rgb]{0,0,0}\alpha=1) approximation w.r.t. increasing $\xi$ and time to maturity $\tau$

In this section, we illustrate the approximation quality for European call options under various model regimes / data set properties as described in Table 1. We use option expiries up to 1Y – we are expecting a loss of approximation quality, based on the nature of the approximation formula. Since we utilized a first order approximation arguments with respect to volatility of volatility, we are also expecting more pronounced differences between MC simulations and the formula for large values of $\xi$.

In Figure 1, we illustrate the approximation quality of the rBergomi approximation for low $\xi$ values. We can observe an expected behaviour: a very good match upto 3M expiry and almost linear deterioration of the quality with increasing $\tau$. Also the approximation formula provides a similar scale of errors across the tested moneyness.

For different moneyness regimes and 1M time to maturity, we obtained the following discrepancies between the MC trials and the introduced formula, measured in the relative option fair value (FV)333Relative FV is the absolute option fair value divided by the initial spot price.:

[TABLE]

In the table above, we can see reasonable approximation error measures which fell below standard 1 MC error for $\xi=10\%$ and $\xi=50\%$ regimes. Due to the theoretical properties of the approximation formula, we observe significant deterioration for high volatility of volatility regimes. This also depends on the time to maturity of the approximated option – the shorter maturity we have, the higher $\xi$ we can allow to obtain reasonable approximation errors (i.e. of the order 1e-04 and lower in terms of FV).

Although the introduced approximation is typically not suitable for calibrations to the whole volatility surface – due to the deterioration of approximation quality when increasing time to maturity – we will illustrate how it can significantly speed-up MC calibration to the provided forward at-the-money (ATMF) backbone.

5.3 Short-tenor calibration and a hybrid calibration to ATMF backbone

Unlike previous analyses, which were based on artificial data / model parameter values, we inspect an application fo the formula on the calibration to real option market data. In particular, we utilize four data sets of AAPL options which were analysed in detail by Pospíšil, Sobotka, and Ziegler (2018). Descriptive statistics of the data sets are provided in Table 2. The following calibration test trials will be considered.

• •

Calibration to short maturity smiles:

This should illustrate how well the model can fit short maturity smiles using the introduced approximation formula without exploiting too high volatility of volatility values ($\xi$). For each data set we selected the shortest maturity slice with more than one traded option. The values were not interpolated by any model, i.e. we calibrated to discrete close mid-prices of traded options. We also confirm, that both MC simulation and the formula reprice the smile with the final calibrated parameters without significant differences.

• •

Hybrid calibration to the ATMF backbone:

In the second trial, we calibrate to each data set ATMF backbone. We note that because we have only a discrete set of traded options we might not have for each maturity an option with strike equal to the corresponding forward. Hence, we take an option with the closest strike to the forward value for each expiry. We use the proposed approximation formula only for $\tau<0.2$, for longer time to maturities we price by MC simulations.

In both cases, the calibration routine was formulated as a standard least-square optimization problem. I.e. to obtain calibrated parameters, we numerically evaluated

[TABLE]

where $N$ is the total number of contracts for the calibration, $\text{Mid}_{i}$ is the mid price of the $i^{th}$ option and $\text{rBergomi}_{i}(\Theta)$ represents the corresponding model price based on parameter set $\Theta$. The model price is either obtained by the approximation formula or by means of MC simulations otherwise. The optimization is performed using Matlab’s local search trust region optimizer which also needs an initial guess to start with.

All following results will be quoted in relative FV: e.g. $\text{rBergomi}_{i}(\Theta)/S_{0}$ and also differences between market and the calibrated model will be denoted using this measure. For the calibration to the whole surface of European options, errors in FV below $0.5\%$ are typically considered to be acceptable, whereas anything exceeding $1\%$ difference is considered as a significant model inconsistency.

Firstly, we display the results for the short-maturity calibration. In Figure 2, we illustrate that even with a not well suited initial guess for Data # 4, we can obtain satisfactory results (i.e. errors significantly below 0.5% mark). In fact, for all tested data sets, obtained values of the calibrated parameters were not very sensitive to the initial guess (only a number of iterations differed). This is a desired feature which typically is not present under classical SV models, see e.g. Mrázek, Pospíšil, and Sobotka (2016). For two other sets we have obtained qualitatively similar fitting errors, however for the data set # 3 we have retrieved 4 errors (out of 22) with absolute FV difference greater than 0.5% and two even greater than 1%, see Figure 3. We conclude that this was partly caused due to high $\xi$ values compared to other three calibration trials and also slightly longer maturity – the data set #3 includes only one option at the shortest maturity, hence the second shortest was used. Still we can conclude that the obtained errors (also verified by using MC simulations) are overall acceptable, although might not be optimal.

We have observed that at least short maturity smiles (< 1M) can be efficiently calibrated using the proposed approximation formula, while for expiries greater than 1M, one would need to stay within a low volatility of volatility regime, otherwise the discrepancies would lead to a non-optimal solution when recomputed using a more precise (but much more costly) MC simulations. To illustrate efficiency of the approximation will show how the approximation can speed-up ATMF-backbone calibration.

We will now inspect the hybrid calibration where we switch between the approximation and a MC pricer based on properties of options being priced. In particular, we focus on data set from 15th May (Data # 4) and for $\tau<0.2$ we will use the approximation formula and MC simulations otherwise. For the calibrated parameters, we will also measure the time spent computing FVs by each pricer. For completeness, we remark that both implementations are not perfect. I.e. for MC simulations under the rBergomi model, the "turbocharging" improvements were introduced by McCrickerd and Pakkanen (2018) recently. On the other hand, numerical integrations in the approximation could be performed by some adaptive quadrature and could be vectorized to improve the computation speed.

In Figure 4, we illustrate calibration fit to the ATMF backbone of the option price surface. We conclude that we have retrieved similar errors for both the prices computed using the proposed approximation and the longer maturity option prices quantified by MC simulations. The final fit of the calibrated model (recomputed by MC simulations) is very good, especially considering that the studied model has only 4 parameters. Moreover, only a fraction of the time spent by MC pricer was needed to computed all FV using the approximation formula. In particular, $98.43\%$ of the pricing time444Excluding any data loading / manipulation routines. we were computing MC simulation estimates of FVs. We also note that $7/9$ of total evaluations were computed by the approximation formula.

6 Conclusion

In previous sections we studied an approximation approach for the pricing of European options under rough stochastic volatility dynamics. Our approach is based on the option price decomposition results obtained by Alòs (2012) for the standard Heston SV model and on its recent generalisation to other SV models by Merino, Pospíšil, Sobotka, and Vives (2018). The main contribution of our research is to derive pricing formulas suitable for various practical applications by applying the general decomposition on a class of Volterra volatility models. Our main focus is laid on the rough volatility models introduced by Gatheral, Jaisson, and Rosenbaum (2018) and Bayer, Friz, and Gatheral (2016).

In particular, a prediction law for Gaussian Volterra processes was proved and an adapted projection of future volatility was obtained in Section 4 for the class of exponential Volterra volatility models, where the volatility process can be expressed as a positive $L^{1}$ function of the time and the Volterra process $Y_{t}$, i.e. $\sigma_{t}=g(t,Y_{t})$. We focused on one particular example of Gaussian Volterra processes, namely on a (rough) fractional Brownian motion $B^{H}_{t}$. "Roughness" of sample paths is determined by the Hurst parameter, $H\in(0,1/2]$, where for $H=0.5$ we would recover a standard Wiener process.

The pricing formula for European options, which is numerically tractable, is then derived under a newly introduced \leavevmode\color[rgb]{0,0,0}\alphaRFSV model, i.e. for

[TABLE]

where \leavevmode\color[rgb]{0,0,0}\alpha\in[0,1] and $r(t)$ is the corresponding autocovariance function defined in Section 4. This newly introduced model seems to have interesting special cases:

• •

for \leavevmode\color[rgb]{0,0,0}\alpha=0 it reverts to a simple exponential RFSV model

• •

and for \leavevmode\color[rgb]{0,0,0}\alpha=1 to rBergomi model respectively.

Both special cases of the \leavevmode\color[rgb]{0,0,0}\alphaRFSV model are studied and praised for their surprising consistency with various financial markets in (Bayer, Friz, and Gatheral 2016) who use the the Bergomi–Guyon expansion to study the ATM skew approximation of the impied volatilities. However, the authors admit that “the Bergomi–Guyon expansion does not converge with values of $\eta$ consistent with the SPX volatility surface, so the Bergomi–Guyon expansion is not useful in practice for calibration of the rBergomi model.” This motived us to derive an approximation formula that will be useful in practice for calibration of the studied models to real market data.

In Section 5, the newly obtained approximation formula was compared to the MC pricing approach introduced by Bennedsen, Lunde, and Pakkanen (2017) and McCrickerd and Pakkanen (2018). This enabled us to numerically verify the obtained solution, to quantify its approximation errors under various settings and, last but not least, to comment on suitability of the rBergomi model for calibration tasks to real market data based on AAPL stock options555These data sets were analysed and described in Pospíšil, Sobotka, and Ziegler (2018). The following conclusions were drawn:

• •

The approximation error is well behaved for short maturities (typically for less than 1M) and the error increases with time to maturity and $\xi$ parameter.

• •

For medium-term expiries we are able to obtain well approximated prices only under low $\xi$ regimes.

• •

Although, the approximation under a rough Volterra process involves several numerical integration procedures, it is much faster than MC simulation approach implemented in the same environment666The numerical trials were implemented in MATLAB environment, see https://www.mathworks.com.. For the standard Wiener case (and in particular for the original Bergomi model), we can have an analytical pricing formula, but also more efficient and well developed MC simulation schemes.

• •

Considering that the modelling approach studied has only few parameters, we were able to fit the sample market data surprisingly well. For calibration to short maturity smiles, we can use just the approximation formula – this was verified by recomputing calibration errors using MC simulations. For calibration to the whole surface, one can utilize a newly introduced hybrid scheme which consist of a combination of approximation and simulation techniques. The idea is quite simple – to use approximation formula for low maturities or for low $\xi$ values and MC simulations for the remaining computations. Suitability of this scheme was judged by a simple calibration to an ATMF-like backbone. We retrieved a well calibrated model, while saving a significant computational time compared to the calibration based on MC simulations only.

However, we remark that implementation of the approximation could be further improved – for simplicity we used a simple trapezoidal quadrature to numerically evaluate integrals appearing in $U_{t}$ and $R_{t}$ expressions.

Another possible improvement could focus on postulating variance dynamics in terms of a forward variance curve – to ensure consistency with variance swaps. This is out-of-scope for the current submission and left for further research. However, some of the main ingredients for forward variance curves under the \leavevmode\color[rgb]{0,0,0}\alphaRFSV model are already introduced in this manuscript.

Also to fit some of the pronounced forward variance curves, one might need a more complex term structure – a simple exponential drift term under the rBergomi model might not be sufficiently flexible. In this case, one can utilize our approach to obtain the approximation formula for a new rough volatility model. In fact, once the new volatility function $g$ is postulated it should be only a matter of algebraic operations to obtain the corresponding approximation formula.

Funding

The work of Jan Pospíšil was partially supported by the Czech Science Foundation (GAČR) grant no. GA18-16680S “Rough models of fractional stochastic volatility”. The work of Josep Vives was partially supported by Spanish grant MEC MTM 2016-76420-P.

Acknowledgements

Computational resources were provided by the CESNET LM2015042 and the CERIT Scientific Cloud LM2015085, provided under the programme “Projects of Large Research, Development, and Innovations Infrastructures”.

Appendix A Decomposition formulas in the general model

In this appendix, we give the error terms for a decomposition of the general model.

The term (I) can be decomposed as

[TABLE]

The term (II) can be decomposed as

[TABLE]

Appendix B Upper-bound for decomposition formula for the exponential Wiener volatility model

In this appendix we obtain the upper-bound for the decomposition formula for the exponential Wiener volatility model in Example 4.12.

B.1 Upper-bound for term (I)

For matter of convenience, we define the function

[TABLE]

If \leavevmode\color[rgb]{0,0,0}\alpha\geq 0, we can find an upper-bound for \chi_{1}(t,T\leavevmode\color[rgb]{0,0,0},\leavevmode\color[rgb]{0,0,0}\alpha) which is

[TABLE]

We can re-write the decomposition formula as

[TABLE]

Applying Lemma 3.1 and using the definition of $a_{u}$, we obtain

[TABLE]

Noting that a_{u}=\sigma_{u}\phi^{1/2}(u,T,\leavevmode\color[rgb]{0,0,0}\alpha), where $\phi$ was defined in (42),

[TABLE]

Being $\sigma_{u}$ the only stochastic component, we can get the expectation inside. Each power of $\sigma_{u}$ has a different forward value, in this case, we can bound all terms by $\exp\left\{\frac{9}{2}\xi^{2}(u-t)\right\}$. We have that

[TABLE]

Substituting \phi(u,T,\leavevmode\color[rgb]{0,0,0}\alpha) and using the upper-bound for \chi_{1}(u,T,\leavevmode\color[rgb]{0,0,0}\alpha) when \leavevmode\color[rgb]{0,0,0}\alpha\geq 0, we have

[TABLE]

The above upper-bound error is difficult to interpret. In order to do this, we derive a Taylor expansion for one of the terms. Then, the following error behaviour is retrieved:

[TABLE]

B.2 Upper-bound for term (II)

For matter of convenience, we define the function

[TABLE]

If \leavevmode\color[rgb]{0,0,0}\alpha\geq 0, we can find an upper-bound for \chi_{2}(t,T,\leavevmode\color[rgb]{0,0,0}\alpha) which is

[TABLE]

We can re-write the decomposition formula as

[TABLE]

Applying Lemma 3.1 and using the definition of $a_{u}$, we obtain

[TABLE]

Noting that a_{u}=\sigma_{u}\phi^{1/2}(u,T,\leavevmode\color[rgb]{0,0,0}\alpha), we have

[TABLE]

Being $\sigma_{u}$ the only stochastic component, we can get the expectation inside. Each power of $\sigma_{u}$ has a different forward value, in this case, we can bound all terms by $\exp\left\{8\xi^{2}(u-t)\right\}$. We have that

[TABLE]

Substituting \phi(u,T,\leavevmode\color[rgb]{0,0,0}\alpha) and using the upper-bound for \chi_{2}(u,T,\leavevmode\color[rgb]{0,0,0}\alpha) when \leavevmode\color[rgb]{0,0,0}\alpha\geq 0, we have

[TABLE]

The above upper-bound error is difficult to interpret. In order to this, we do a Taylor analysis of one term. Then, the following error behavior is retrieved

[TABLE]