On VIX Futures in the rough Bergomi model

TL;DR

This paper explores the dynamics of VIX and forward variance in the rough Bergomi model, developing efficient pricing methods for VIX derivatives and assessing its ability to jointly model VIX and SPX markets.

Contribution

It introduces new algorithms for pricing VIX futures and options within the rough Bergomi model and evaluates its joint calibration to VIX and SPX data.

Findings

The model accurately reproduces implied volatility smiles for short maturities.

Efficient algorithms for VIX derivatives pricing are developed.

The joint calibration demonstrates the model's capability to fit both VIX and SPX data.

Abstract

The rough Bergomi model introduced by Bayer, Friz and Gatheral has been outperforming conventional Markovian stochastic volatility models by reproducing implied volatility smiles in a very realistic manner, in particular for short maturities. We investigate here the dynamics of the VIX and the forward variance curve generated by this model, and develop efficient pricing algorithms for VIX futures and options. We further analyse the validity of the rough Bergomi model to jointly describe the VIX and the SPX, and present a joint calibration algorithm based on the hybrid scheme by Bennedsen, Lunde and Pakkanen.

Equations

$$\begin{array}[]{rll}\mathrm{d}X_{t}&=\displaystyle-\frac{1}{2}V_{t}\mathrm{d}t+\sqrt{V_{t}}\mathrm{d}W_{t},&X_{0}=0\\ V_{t}&=\displaystyle\xi_{0}(t)\mathcal{E}(2\nu C_{H}\mathcal{V}_{t}),&V_{0}>0,\end{array}$$

$${\mathcal{V}}_t:={\int }_0^t(t-u{)}^{H_{-}}\mathrm{d}{Z}_u,$$

$$\mathbb{E}({\mathcal{V}}_t{\mathcal{V}}_s)=s^{2H}{\int }_0^1{\left(\frac{t}{s}-u\right)}^{H_{-}}(1-u{)}^{H_{-}}\mathrm{d}u,\text{for any }s,t\in [0,1].$$

$${\mathcal{V}}_{t,T}:={\int }_T^t(t-u{)}^{H_{-}}\mathrm{d}{Z}_u\text{and}{\mathcal{V}}_t^T:={\int }_0^T(t-u{)}^{H_{-}}\mathrm{d}{Z}_u.$$

$${\sigma }_T^2(t)=\frac{1}{t-T}{\int }_T^t{\xi }_T(u)\mathrm{d}u,$$

$$\mathbb{V}({\mathcal{V}}_t^T)=\frac{t^{2H}-(t-T{)}^{2H}}{2H},\text{for }t\ge T,$$

$$\mathbb{E}\left({\mathcal{V}}_t^T{\mathcal{V}}_s^T\right)={\int }_0^T{\left[(t-u)(s-u)\right]}^{H_{-}}\mathrm{d}u=\frac{(s-t{)}^{H_{-}}}{H_{+}}\left\{t^{H_{+}}\mathrm{F}\left(\frac{-t}{s-t}\right)-(t-T{)}^{H_{+}}\mathrm{F}\left(\frac{T-t}{s-t}\right)\right\},$$

$$\mathrm{F}(u):={}_2{F}_1\left(-H_{-},H_{+},1+H_{+},u\right),$$

$$[\mathcal{V}{]}_t=\frac{t^{2H}}{2H},\text{for }t\ge 0.$$

$$\left(\begin{array}{c}\mathcal{V}\left(\frac{1}{n}\right)\\ \mathcal{V}\left(\frac{2}{n}\right)\\ ⋮\\ \mathcal{V}\left(\frac{n_T}{n}\right)\end{array}\right)=\left(\begin{array}{ccccc}{\overline{Z}}_{0,1}& 0& \cdots & \cdots & 0\\ {\overline{Z}}_{1,1}& {\overline{Z}}_0& \ddots & \ddots & 0\\ {\overline{Z}}_{2,1}& {\overline{Z}}_1& \ddots & \ddots & ⋮\\ ⋮& \ddots & \ddots & \ddots & ⋮\\ {\overline{Z}}_{n_T-1,1}& {\overline{Z}}_{n_T-2}& \cdots & {\overline{Z}}_1& {\overline{Z}}_0\end{array}\right)\left(\begin{array}{c}1\\ {\left(\frac{1}{n}{b}_1^{\ast }\right)}^{H_{-}}\\ {\left(\frac{1}{n}{b}_2^{\ast }\right)}^{H_{-}}\\ ⋮\\ {\left(\frac{1}{n}{b}_{n_T-1}^{\ast }\right)}^{H_{-}}\end{array}\right),$$

$${\mathrm{VIX}}_T^2:=\mathbb{E}\left(\frac{1}{\Delta }{\int }_T^{T+\Delta }\mathrm{d}⟨X_s,X_s⟩\mathrm{d}s|{\mathcal{F}}_T\right),$$

$${\mathfrak{V}}_T:=\mathbb{E}\left({\mathrm{VIX}}_T|{\mathcal{F}}_0\right)=\mathbb{E}\left(\sqrt{\frac{1}{\Delta }{\int }_T^{T+\Delta }\mathbb{E}\left(\mathrm{d}⟨X_s,X_s⟩|{\mathcal{F}}_T\right)\mathrm{d}s}|{\mathcal{F}}_0\right)=\mathbb{E}\left(\sqrt{\frac{1}{\Delta }{\int }_T^{T+\Delta }{\xi }_T(s)\mathrm{d}s}|{\mathcal{F}}_0\right).$$

$${\eta }_T(t):=\exp \left(2\nu C_H{\mathcal{V}}_t^T\right).$$

$${\mathrm{VIX}}_T={\left\{\frac{1}{\Delta }{\int }_T^{T+\Delta }{\xi }_0(t){\eta }_T(t)\exp \left(\frac{{\nu }^2{C}_H^2}{H}\left[(t-T{)}^{2H}-t^{2H}\right]\right)\mathrm{d}t\right\}}^{1/2}.$$

$${\mathrm{VIX}}_T^2$$

$$=\frac{1}{\Delta }{\int }_T^{T+\Delta }\mathbb{E}\left[{\xi }_0(t){\eta }_T(t)\exp \left(2\nu C_H{\mathcal{V}}_{t,T}-\frac{{\nu }^2{C}_H^2{t}^{2H}}{H}\right)|{\mathcal{F}}_T\right]\mathrm{d}t,$$

$${\mathrm{VIX}}_T^2=\frac{1}{\Delta }{\int }_T^{T+\Delta }{\xi }_0(t){\eta }_T(t)\mathbb{E}\left[\exp \left(2\nu C_H{\mathcal{V}}_{t,T}-\frac{{\nu }^2{C}_H^2{t}^{2H}}{H}\right)|{\mathcal{F}}_T\right]\mathrm{d}t.$$

$${\xi }_T(t)={\xi }_0(t){\eta }_T(t)\exp \left(\frac{{\nu }^2{C}_H^2}{H}\left[(t-T{)}^{2H}-t^{2H}\right]\right),\text{for any }t\ge T.$$

$$\mathbb{E}\left(V_t|{\mathcal{F}}_T\right)={\xi }_0(t){\eta }_T(t)\exp \left(\frac{{\nu }^2{C}_H^2}{H}\left[(t-T{)}^{2H}-t^{2H}\right]\right).$$

$$\frac{1}{\Delta }{\int }_T^{T+\Delta }\sqrt{{\xi }_0(t)}\exp \left(\frac{{\nu }^2{C}_H^2}{4H}\left[(t-T{)}^{2H}-t^{2H}\right]\right)\mathrm{d}t\le {\mathfrak{V}}_T\le {\left\{\frac{1}{\Delta }{\int }_T^{T+\Delta }{\xi }_0(s)\mathrm{d}s\right\}}^{1/2}.$$

$${\mathfrak{V}}_T=\mathbb{E}\left({\mathrm{VIX}}_T|{\mathcal{F}}_0\right)=\mathbb{E}\left(\sqrt{\frac{1}{\Delta }{\int }_T^{T+\Delta }{\xi }_T(s)\mathrm{d}s}|{\mathcal{F}}_0\right)\le \sqrt{\mathbb{E}\left(\frac{1}{\Delta }{\int }_T^{T+\Delta }{\xi }_T(s)\mathrm{d}s|{\mathcal{F}}_0\right)}.$$

$${\mathfrak{V}}_T$$

$$\ge \mathbb{E}\left[\frac{1}{\Delta }{\int }_T^{T+\Delta }\sqrt{{\xi }_0(t){\eta }_T(t)}\exp \left(\frac{{\nu }^2{C}_H^2}{2H}\left[(t-T{)}^{2H}-t^{2H}\right]\right)\mathrm{d}t|{\mathcal{F}}_0\right]$$

$$=\frac{1}{\Delta }{\int }_T^{T+\Delta }\sqrt{{\xi }_0(t)}\mathbb{E}\left(\sqrt{{\eta }_T(t)}\right)\exp \left(\frac{{\nu }^2{C}_H^2}{2H}\left[(t-T{)}^{2H}-t^{2H}\right]\right)\mathrm{d}t$$

$$=\frac{1}{\Delta }{\int }_T^{T+\Delta }\sqrt{{\xi }_0(t)}\exp \left(\frac{{\nu }^2{C}_H^2}{4H}\left[t^{2H}-(t-T{)}^{2H}\right]\right)\exp \left(\frac{{\nu }^2{C}_H^2}{2H}\left[(t-T{)}^{2H}-t^{2H}\right]\right)\mathrm{d}t$$

$$=\frac{1}{\Delta }{\int }_T^{T+\Delta }\sqrt{{\xi }_0(t)}\exp \left(\frac{{\nu }^2{C}_H^2}{4H}\left[(t-T{)}^{2H}-t^{2H}\right]\right)\mathrm{d}t,$$

$$\text{Scenario }1:{\xi }_0(t)=0.234^2;\text{Scenario }2:{\xi }_0(t)=0.234^2(1+t{)}^2;\text{Scenario }3:{\xi }_0(t)=0.234^2\sqrt{1+t}.$$

$$\varphi (T):=\sqrt{\frac{1}{\Delta }{\int }_T^{T+\Delta }f(t)\mathrm{d}t}=\sqrt{f(T)+\frac{\Delta }{2}{f}^{\prime }(T)+\frac{{\Delta }^2}{6}{f}^{\prime \prime }(T)+\mathcal{O}({\Delta }^3)},$$

$$\psi (T):=\frac{1}{\Delta }{\int }_T^{T+\Delta }\sqrt{f(t)}\mathrm{d}t=\sqrt{f(T)}+\frac{\Delta }{4}\frac{f^{\prime }(T)}{\sqrt{f(T)}}+\mathcal{O}({\Delta }^2),$$

$$\varphi (T{)}^2-\psi (T{)}^2={\Delta }^2\left(\frac{f^{\prime \prime }(T)}{6}-\frac{f^{\prime }(T{)}^2}{16f(T)}\right)+\mathcal{O}({\Delta }^3).$$

$$\begin{array}[]{rll}Z_{t_{i}}&=Z_{t_{i-1}}+n^{H_{-}}\left(\mathcal{V}(t_{i})-\mathcal{V}(t_{i-1})\right),&\text{for }i=1,\ldots,\kappa,\\ Z_{t_{i}}&=Z_{t_{i-1}}+\overline{Z}_{i-1},&\text{for }i>\kappa;\end{array}$$

$${\tilde{\mathcal{V}}}_{{\tau }_0}^T:={\mathcal{V}}_T^T\text{and}{\tilde{\mathcal{V}}}_{{\tau }_j}^T:=\sum _{i=1}^{n_T}\frac{Z_{t_i}-Z_{t_{i-1}}}{({\tau }_j-t_{i-1}{)}^{-H_{-}}},\text{for }j=1,\dots ,N;$$

$${\mathrm{VIX}}_T\approx {\left\{\frac{1}{\Delta }\sum _{j=0}^{N-1}\frac{Q_{T,{\tau }_j}^2+Q_{T,{\tau }_{j+1}}^2}{2}({\tau }_j-{\tau }_{j-1})\right\}}^{1/2},$$

Full text

The authors would like to thank Christian Bayer, Jim Gatheral, Mikko Pakkanen and Mathieu Rosenbaum for useful discussions. AJ acknowledges financial support from the EPSRC First Grant EP/M008436/1. The numerical implementations have been carried out on the collaborative platform Zanadu (www.zanadu.io).

1. Introduction

Volatility, though not directly observed nor traded, is a fundamental object on financial markets, and has been the centre of attention of decades of theoretical and practical research, both to estimate it and to use it for trading purposes. The former goal has usually been carried out under the historical measure ($\mathbb{P}$) while the latter, through the introduction of volatility derivatives (VIX and related family), has been evolving under the pricing measure $\mathbb{Q}$. Most models used for pricing purposes (Heston [19], SABR [17], Bergomi [5]) are constructed under $\mathbb{Q}$ and are of Markovian nature (making pricing, and hence calibration, easier). Recently, Gatheral, Jaisson and Rosenbaum [14] broke this routine and introduced a fractional Brownian motion as driving factor of the volatility process. This approach (Rough Fractional Stochastic Volatility, RFSV for short) opens the door to revisiting classical pricing and calibration conundrums. They, together with the subsequent paper by Bayer, Friz and Gatheral, (see also [1, 12]) in particular showed that these models were able to capture the extra steepness of the implied volatility smile in Equity markets for short maturities, which continuous Markovian stochastic volatility models fail to describe. The icing on the cake is the (at last!!) reconciliation between the two measures $\mathbb{P}$ and $\mathbb{Q}$ within a given model, showing remarkable results both for estimation and for prediction.

One of the key issues in Equity markets is, not only to fit the (SPX) implied volatility smile, but to do so jointly with a calibration of the VIX (Futures and ideally options). Gatheral’s [15] double mean reverting process is the leading (Markovian) continuous model in this direction, while models with jumps have been proposed abundantly by Carr and Madan [9] and Kokholm and Stisen [20]. This issue was briefly tackled by Bayer, Friz and Gatheral [3] for a particular rough model (rough Bergomi), and we aim here at providing a deeper analysis of VIX dynamics under this rough model and at implementing pricing schemes for VIX Futures and options. Our main contribution is a precise link between the forward variance curve $(\xi_{T}(\cdot))_{T\geq 0}$ and the initial forward variance curve $\xi_{0}(\cdot)$ in the rough Bergomi model. This in turn, allows us not only to provide simulation methods for the VIX, but also to refine the log-normal approximation of [3] for VIX Futures, matching exactly the first two moments. Finally, we develop an efficient algorithm for VIX Futures calibration, upon which we build a joint calibration method with the SPX. As opposed to the Cholesky approach in [3], we adapt the hybrid-scheme by Bennedsen, Lunde and Pakkanen [4] with better complexity $\mathcal{O}(n\log n)$ . Assuming the universality of the Hurst parameter $H$ across VIX and SPX allows us to compute efficiently prices recursively with complexity $\mathcal{O}(n)$. In passing, we also investigate the joint consistency of VIX and SPX in the market. The organisation of the paper follows accordingly: we first introduce the rough Bergomi model and its main properties (Section 2), before presenting its pricing power for VIX Futures (Section 3), and finally develop the joint calibration algorithm in Section 4.

2. Rough volatility and the rough Bergomi model

Comte and Renault[7] were the first to propose a stochastic volatility model in which the instantaneous volatility is driven by a fractional Brownian motion, with a Hurst index restricted to be greater than $1/2$. Recently Gatheral, Jaisson and Rosenbaum [14] presented a new approach with a Hurst index smaller than $1/2$, producing extremely good fits to observed volatility data under the physical measure $\mathbb{P}$. These models form the so-called Rough Fractional Stochastic Volatility (RFSV) family that is understood as a natural extension of the classical volatility models driven by standard Brownian motion. Our work focuses on the pricing measure $\mathbb{Q}$ and we assume through this paper that the model presented by Gatheral, Jaisson and Rosenbaum [14] under $\mathbb{P}$ is a reasonable model. Finally, and most importantly, we follow the recent paper by Bayer, Friz and Gatheral [3], in order to extend the RFSV model to pricing schemes under the measure $\mathbb{Q}$. More precisely, Bayer, Friz and Gatheral [3] proposed the following model for the log stock price process $X:=\log(S)$:

[TABLE]

with $\nu,\xi_{0}(\cdot)>0$, $\mathcal{E}(\cdot)$ is the Doléans-Dade [10] stochastic exponential and $C_{H}:=\sqrt{\frac{2H\Gamma(2-H_{+})}{\Gamma(H_{+})\Gamma(2-2H)}}$, where, for notational convenience (throughout the paper), we use the symbols $H_{\pm}:=H\pm\frac{1}{2}$. All the processes are defined on a given filtered probability space $(\Omega,\mathcal{F},(\mathcal{F}_{t})_{t\geq 0},\mathbb{Q})$ supporting the two standard Brownian motions $W$ and $Z$ (see below). The initial forward variance curve is observed at inception, and we therefore assume without loss of generality that it is $\mathcal{F}_{0}$-measurable. The process $\mathcal{V}$, defined as

[TABLE]

is a centred Gaussian process with covariance structure

[TABLE]

We shall also introduce, for any $0\leq T\leq t$ the notations

[TABLE]

Note in particular that $\mathcal{V}^{T}_{T}=\mathcal{V}_{T}$. The two standard Brownian motions $W$ and $Z$ are correlated with correlation parameter $\rho\in(-1,1)$. Here, $\xi_{T}(t)$ denotes the forward variance observed at time $T$ for a maturity equal to $t$. More precisely, if $\sigma^{2}_{T}(t)$ denotes the fair strike of a variance swap observed at time $T$ and maturing at $t$, then

[TABLE]

or equivalently $\xi_{T}(t)=\frac{\mathrm{d}}{\mathrm{d}t}\left((t-T)\sigma^{2}_{T}(t)\right)$. For any fixed $t>0$, the process $(\xi_{s}(t))_{s\leq t}$, is a martingale, i.e. $\mathbb{E}[\xi_{s}(t)|\mathcal{F}_{u}]=\xi_{u}(t)$, for all $u\leq s\leq t$. Furthermore, $\mathcal{V}^{T}_{t}$ is a centred Gaussian process with variance

[TABLE]

and covariance structure

[TABLE]

for any $t<s$, where we introduce the function $\mathrm{F}:\mathbb{R}_{-}\to\mathbb{R}$ as

[TABLE]

and $\text{}_{2}F_{1}$ is the hypergeometric function [2, Chapter 15]. Finally, the quadratic variation of $\mathcal{V}$ is given by

[TABLE]

2.1. Hybrid simulation scheme

Bayer, Friz and Gatheral [3] present a Cholesky method to simulate the rough Bergomi model. Although exact, this method is very slow and other approaches need to be considered for calibration purposes. Recently, Bennedsen, Lunde and Pakkanen [4] presented a new simulation scheme for Brownian semistationary ($\mathcal{BSS}$) processes. This method, as opposed to Cholesky, is an approximate method. However, in [4] the authors show that the method yields remarkable results in the case of the rough Bergomi model. In addition, their approach leads to a natural simulation of both the Volterra process $\mathcal{V}$ and the stock price $S$ and yields a computational complexity of order $\mathcal{O}(n\log n)$.

Definition 2.1**.**

Let $W$ be a standard Brownian motion on a given filtered probability space $(\Omega,\mathcal{F},(\mathcal{F}_{t})_{t\geq 0},\mathbb{P})$. A truncated Brownian semistationary ($\mathcal{BSS}$) process is defined as $\mathcal{B}(t)=\int_{0}^{t}g(t-s)\sigma(s)\mathrm{d}W_{s}$, for $t\geq 0$, where $\sigma$ is $(\mathcal{F}_{t})_{t\geq 0}$-predictable with locally bounded trajectories and finite second moments, and $g:(0,\infty)\to[0,\infty)$ is Borel measurable and square integrable. We shall call it a $\mathcal{BSS}(\alpha,W)$ process if furthermore

1. (i)

there exists $\alpha\in\left(-\frac{1}{2},\frac{1}{2}\right)\setminus\{0\}$ such that $g(x)=x^{\alpha}L_{g}(x)$ for all $x\in(0,1]$, where $L_{g}\in\mathcal{C}^{1}((0,1]\to[0,\infty))$, is slowly varying111A measurable function $L:(0,1]\to[0,\infty)$ is slowly varying [6] at [math] if for any $t>0$, $\lim\limits_{x\downarrow 0}L(tx)/L(x)=1.$ at the origin and bounded away from zero. Moreover, there exists a constant $C>0$ such that $|L^{\prime}_{g}(x)|\leq C(1+x^{-1})$ for all $x\in(0,1]$; 2. (ii)

the function $g$ is differentiable on $(0,\infty)$.

Under this assumption, the hybrid scheme, proposed in [4] and recalled in Appendix A, provides an efficient way to simulate $\mathcal{BSS}$ processes. It applies in particular to the rough Bergomi model:

Proposition 2.2**.**

The Volterra process $\mathcal{V}$ in (2.2) is a truncated $\mathcal{BSS}(H_{-},Z)$ process.

Proof.

From (2.2), $g(x)\equiv x^{H_{-}}$ and $\sigma(\cdot)\equiv 1$ as in Definition 2.1, so that $V$ is a $\mathcal{BSS}$ process. Since $H_{-}\in(-\frac{1}{2},0)$, then $L_{g}\equiv 1$, and $\mathcal{V}$ satisfies Definition 2.1(i); Definition 2.1(i) trivially holds, and so does the corollary. ∎

The corollary implies that we can apply the hybrid scheme to $\mathcal{V}$. In particular, for $\kappa=1$ the matrix form representation of the scheme reads (recall that $n_{T}:=\lfloor nT\rfloor$)

[TABLE]

where the coefficients $\{b^{*}_{i}\}$ are defined in (A.1). This matrix multiplication is, by brute force, of order $\mathcal{O}(n^{2})$, however using discrete convolution we may use FFT to reduce it to $\mathcal{O}(n\log n)$ as suggested in [4] .

3. Rough Bergomi and VIX

Bayer, Friz and Gatheral [3] briefly discuss the lack of consistency of the rough Bergomi model with observed VIX options data, leading to an incorrect term structure of the VIX. In this section, we investigate in detail the dynamics of the VIX, and propose a log-normal approximation. Additionally, we investigate the viability of the model in terms of VIX Futures and options, and compare it to the approximation in [3].

3.1. VIX Futures in the rough Bergomi model.

From now on, we fix a given maturity $T\geq 0$, and define the VIX at time $T$ via the continuous-time monitoring formula

[TABLE]

where $\Delta$ is equal to $30$ days. The risk-neutral formula for the VIX future $\mathfrak{V}_{T}$ with maturity $T$ is then given by

[TABLE]

Note that, when $T>0$, $\xi_{T}(s)$ is a market input which is not $\mathcal{F}_{0}$-measurable, and is hence difficult to interpret only knowing $\mathcal{F}_{0}$. We shall make repeated use of the following random variable defined for any $t\geq T$, by

[TABLE]

Proposition 3.1**.**

The VIX dynamics are given by

[TABLE]

Proof.

Using Fubini’s theorem and the instantaneous variance representation in (2.1), we can write

[TABLE]

with $\mathcal{V}_{t,T}$ defined in (2.3). Since $\eta_{T}(t)\in\mathcal{F}_{T}$ and $\xi_{0}(t)\in\mathcal{F}_{0}$, this expression simplifies to

[TABLE]

The proposition follows since $\mathcal{V}_{t,T}$ is centred Gaussian, independent of $\mathcal{F}_{T}$, with variance given in (2.4), and $\mathbb{E}\left(\left.\mathrm{e}^{2\nu C_{H}\mathcal{V}_{t,T}}\right|\mathcal{F}_{T}\right)=\mathbb{E}\left(\mathrm{e}^{2\nu C_{H}\mathcal{V}_{t,T}}\right)=\exp\left(\frac{\nu^{2}C_{H}^{2}}{H}(t-T)^{2H}\right).$ ∎

The main challenge for simulation is $\eta_{T}(t)$. However, since the latter is independent of $\xi_{0}(\cdot)$, robustness of simulation schemes for the VIX will not be affected by the qualitative properties of the initial variance curve $\xi_{0}$.

Proposition 3.2**.**

The forward variance curve $\xi_{T}$ in the rough Bergomi model admits the representation

[TABLE]

Proof.

Since $\mathbb{E}\left(V_{t}|\mathcal{F}_{T}\right)=\xi_{T}(t)$ by (3.1), the proposition follows from Proposition 3.1 and the equality

[TABLE]

∎

Bayer, Friz and Gatheral [3] did not derive such a representation for $\xi_{T}$, and their approach for pricing VIX derivatives relies on an approximation which avoids the computations developed in this section. Proposition 3.2 allows for a better understanding of the process $\xi_{T}$, and for an innovative approach to price VIX derivatives.

3.2. Upper and lower bounds for VIX Futures

Theorem 3.3**.**

The following bounds hold for VIX Futures:

[TABLE]

Proof.

The conditional Jensen’s inequality gives

[TABLE]

Furthermore, since $\xi_{0}$ is $\mathcal{F}_{0}$-adapted, Fubini’s theorem along with the martingale property of $\xi_{T}$ yield the upper bound $\mathfrak{V}_{T}=\mathbb{E}\left(\mathrm{VIX}_{T}|\mathcal{F}_{0}\right)\leq\sqrt{\Delta^{-1}\int_{T}^{T+\Delta}\xi_{0}(s)\mathrm{d}s}$. To obtain a lower bound we use the representation in Proposition 3.2, and Cauchy-Schwarz’s inequality, and Fubini’s theorem, so that

[TABLE]

since $\eta_{T}(t)^{1/2}$ is log-normal (Proposition 3.1), so that $\mathbb{E}(\sqrt{\eta_{T}(t)})=\exp\left(\frac{\nu^{2}C_{H}^{2}}{4H}\left[t^{2H}-(t-T)^{2H}\right]\right)$. ∎

We perform a numerical experiment to check the tightness of the bounds obtained in Proposition 3.3). For this analysis we consider three qualitative scenarios for the initial forward variance curve:

[TABLE]

Figures 1 suggest that the lower bound given in Proposition 3.3) is surprisingly tight for very different shapes of $\xi_{0}$. This can be explained with the following argument: consider a simplified and deterministic version of the VIX futures price in Proposition 3.1), denoted by

[TABLE]

for some strictly positive (deterministic) function $f\in\mathcal{C}^{2}(\mathbb{R})$. We further introduce

[TABLE]

which is the corresponding lower bound by Cauchy-Schwarz’s (or Jensen’s) inequality, so that

[TABLE]

Hence, we observe that for small $\Delta$, as it is the case in VIX futures, the lower bound is very close from the original value which explains (at least for the deterministic case) the behaviour observed in Figure 1.

3.3. Numerical implementation of VIX process

In this section, we investigate different simulation schemes for the VIX in the rough Bergomi model.

3.3.1. Hybrid scheme and forward Euler approach

In order to simulate the process $(\mathcal{V}_{t}^{T})_{t\in[T,T+\Delta]}$ it is important to notice from (2.3) that the kernel has a singularity only for $\mathcal{V}^{T}_{T}$, hence we may overcome this by simulating the process using the hybrid scheme from Section 2.1. Then, we may easily extract $Z$ and simulate $(\mathcal{V}^{T}_{t})_{t\in(T,T+\Delta]}$ using the forward Euler scheme with complexity $\mathcal{O}(n)$. A forward Euler scheme is chosen in this case due to the fact that the kernel in (2.3) no longer has a singularity for $t\in(T,T+\Delta]$ and this method is faster than the hybrid scheme. Once $(\mathcal{V}_{t}^{T})_{t\in[T,T+\Delta]}$ is simulated, numerical integration routines may be used to simulate the VIX process using the expression in Proposition 3.1. It must be pointed out that this approach is computationally expensive and memory consuming since it involves to simulate the Volterra process using the hybrid scheme with complexity of $\mathcal{O}(n\log n)$ and additionally, each $\mathcal{V}^{T}_{t}$ by forward Euler. The simulation algorithm can be summarised as follows:

Algorithm 3.4** (VIX simulation in the rough Bergomi model).**

Fix a grid $\mathcal{T}=\{t_{i}\}_{i=0,\ldots,n_{T}}$ and $\kappa\geq 1$.

1. (1)

Simulate the Volterra process $(\mathcal{V}_{t})_{t\in[0,T]}$ using the hybrid scheme in Appendix A, yielding a sample of the random variable $\mathcal{V}_{T}^{T}=\mathcal{V}_{T}$; 2. (2)

extract the path of the Brownian motion $Z$ driving the Volterra process.

[TABLE] 3. (3)

fix a grid $\mathfrak{T}=\{\tau_{j}\}_{j=0,...,N}$ on $[T,T+\Delta]$ and approximate the continuous-time process $\mathcal{V}^{T}$ by the discrete-time version $\widetilde{\mathcal{V}}^{T}$ defined via the following forward Euler scheme:

[TABLE] 4. (4)

compute the VIX process via numerical integration, for example using a composite trapezoidal rule:

[TABLE]

where $Q^{2}_{T,\tau_{j}}:=\xi_{0}(\tau_{j})\exp\left(2vC_{H}\widetilde{\mathcal{V}}^{T}_{\tau_{j}}\right)\exp\left(\frac{\nu^{2}C_{H}^{2}}{H}\left((\tau_{j}-T)^{2H}-\tau_{j}^{2H}\right)\right)$.

Remark 3.5**.**

Step 4 may obviously be replaced by any available numerical integration routine, but one must then carefully choose the partition in Step 3.

3.3.2. Truncated Cholesky approach

Alternatively, one could use the more expensive, yet exact, Cholesky decomposition to simulate $\mathcal{V}^{T}$ on $[T,T+\Delta]$ since its covariance structure is known from (2.5). However, computational complexity aside, with the same grid $\mathcal{T}$ as in Algorithm 3.4, numerical experiments suggest that the determinant of the covariance matrix is equal to zero when using more than $n_{T}=8$ discretisation points. Hence, although valid in theory, the Cholesky approach is not feasible numerically. This fact implies that there exists strong linear dependence. In fact, for any $\varepsilon>0$, the strict inequality $\text{corr}(\mathcal{V}^{T}_{t_{1}},\mathcal{V}^{T}_{t_{1}+\varepsilon})<\text{corr}(\mathcal{V}^{T}_{t_{2}},\mathcal{V}^{T}_{t_{2}+\varepsilon})$ holds for all $T<t_{1}<t_{2}$, as well as the following:

Proposition 3.6**.**

The limit $\lim\limits_{\varepsilon\downarrow 0}\text{corr}(\mathcal{V}^{T}_{t},\mathcal{V}^{T}_{t+\varepsilon})=1$ holds for any $t\in[T,T+\Delta]$.

Proof.

This follows readily from the continuity in $\L^{2}$ of the map $t\mapsto\mathcal{V}^{T}_{t}$:

[TABLE]

Applying the identity [2, page 564]

[TABLE]

to the case $(a,b,c)=(-H_{+},1+H_{+},2+H_{+})$ and using properties of the Gamma function along with the fact that $\text{}_{2}F_{1}\left(a,0,c,z\right)\equiv 1$, we obtain

[TABLE]

Finally, the series representation of $\text{}_{2}F_{1}$ [2, Chapter 15.1.1] implies that $\text{}_{2}F_{1}\left(-H_{+},2+2H_{+},-2H_{+},\frac{-\varepsilon}{T}\right)$ converges to $1$ as $\varepsilon$ tends to zero, and hence that $E[(\mathcal{V}^{T}_{t+\varepsilon}-\mathcal{V}^{T}_{t})^{2}]$ tends to zero. ∎

In light of Proposition 3.6, we model exactly the dependence structure on the first 8 grid points $t_{1},\ldots,t_{8}$, then truncate the Cholesky decomposition and compute the correlations $\rho_{i}:=\text{corr}(\mathcal{V}^{T}_{t_{8+i}},\mathcal{V}^{T}_{t_{8+i+1}})$, for $i=0,\ldots,n-9$ to approximate the process by adequately rescaling and correlating each pair of subsequent grid points. In contrast to the hybrid + forward Euler scheme, the computational complexity is much lower, since the Cholesky method is truncated with only 8 components. The VIX simulation algorithm therefore reads as follows:

Algorithm 3.7** (VIX simulation (truncated Cholesky)).**

Fix a grid $\mathfrak{T}=\{\tau_{j}\}_{j=0,...,N}$ on $[T,T+\Delta]$,

1. (i)

compute the covariance matrix of $(\mathcal{V}_{\tau_{j}}^{T})_{i=j,\ldots,8}$ using (2.5); 2. (ii)

generate $\{\mathcal{V}^{T}_{\tau_{j}}\}_{j=1,...,N}$ by correlating and rescaling using (2.5):

[TABLE] 3. (iii)

compute the VIX via numerical integration as in Algorithm 3.4(4).

3.3.3. Numerical experiment

We compute the price of VIX Futures using the simulation algorithm introduced in the previous section. We set the same parameters as in [3] and [4]:

[TABLE]

We perform $10^{5}$ simulations for the hybrid scheme + forward Euler (HS+FE) method, while $10^{6}$ simulations are used for the truncated Cholesky. Figure 2 suggests that both methods agree qualitatively and converge to a similar output. In particular, the truncated Cholesky approach seems to suffer from larger oscillations as maturity increases, even with $10^{6}$ simulations. Nevertheless, Figure 2 indicates that the Monte-Carlo variance increases in $T$ for both schemes, which is confirmed in Figure 3, where the error also increases with maturity. On the other hand, Figure 4 shows that the HS+FE methods is slower than the truncated Cholesky method, which is consistent with the computational complexities discussed in the previous section. In particular, the computational time of the Cholesky method is almost constant when using parallel computing. Figure 4 also suggests that large simulations are needed to obtain precise prices. In light of this analysis, both methods seem to approximate the required output in a decent manner. Even if the truncated Cholesky approach gives a considerably fast output for each maturity, it is not considered for calibration, since its computational time grows linearly in the number of maturities, making the algorithm too slow for reasonable calibration. Instead, we will use the truncated Cholesky approach as a benchmark for the upcoming approximations.

3.4. The VIX process and log-normal approximations

We now investigate approximate methods to price VIX futures and options. We define the $\mathcal{F}_{T}$-measurable random variable $\Delta\mathrm{VIX}^{2}_{T}=\int_{T}^{T+\Delta}\xi_{T}(t)\mathrm{d}t$. In [3], Bayer, Friz and Gatheral assumed that $\log(\Delta\mathrm{VIX}^{2}_{T})$ follows a Gaussian distribution, and computed directly its first and second moments, which hence fully characterise the distribution of $\Delta\mathrm{VIX}^{2}_{T}$. However, since the sum of log-normal random variables is not (in general) log-normal, the fact that $\xi_{T}$ is log-normal does not imply that $\Delta\mathrm{VIX}^{2}_{T}$ is. In a different context (geometric Brownian motion and Asian options), Dufresne [11] proved that, under certain conditions, an integral of log-normal variables asymptotically converges to a log-normal. This approximation has been widely used for many applications [11], and Dufresne’s result motivates Bayer-Friz-Gatheral’s assumption. We provide here exact formulae for the mean and variance of this distribution, and compare them numerically to those by Bayer-Friz-Gatheral.

Proposition 3.8**.**

The following holds:

[TABLE]

Furthermore, $\mathbb{E}(\Delta\mathrm{VIX}^{2}_{T})=\int_{T}^{T+\Delta}\xi_{0}(t)\mathrm{d}t$ and

[TABLE]

where $\overline{\Theta}_{u,t}$ is equal to zero if $u=t$ and otherwise equal to $\Theta_{u\vee t,u\wedge t}$, where

[TABLE]

Remark 3.9**.**

Since all the integrals in the proposition are computed over compact intervals, they are finite as long as $\xi_{0}$ is well behaved on $[T,T+\Delta]$. Assuming this is indeed the case is not restrictive in practice as $\xi_{0}$ represents the initial forward variance curve; note in particular that continuity of $\xi_{0}$ is sufficient.

Remark 3.10**.**

The reason why $\overline{\Theta}_{u,t}$ is defined that way is for numerical purposes. Indeed, when $u<t$, $\Theta_{u,t}$ is not well defined since $\mathrm{F}(x)$ only makes sense when $x\leq 1$ (details on the radius of convergence of hypergeometric functions can be found in [2][Page 556]), even though the integral representation (3.6) is still well defined. However, most numerical packages implement hypergeometric functions via series expansions. The trick from $\Theta_{u,t}$ to $\overline{\Theta}_{u,t}$ allows us to bypass this issue.

Proof of Proposition 3.8.

The expectation follows directly from the tower property and Fubini’s theorem. For the second moment, we use the decomposition

[TABLE]

where in the last step we use that $\xi_{T}(s)$ is $\mathcal{F}_{T}$-measurable in order to apply Fubini. Using the representation obtained in Proposition 3.2, we get

[TABLE]

where the random variable

[TABLE]

is Gaussian with zero expectation and variance

[TABLE]

Then,

[TABLE]

and the second moment follows from the immediate computation

[TABLE]

after defining $\Theta_{u,t}:=\exp\left(\frac{1}{2}\mathbb{V}(\vartheta_{u,t})\right)$. ∎

In order to provide closed-form expressions for VIX Futures and options, we enforce the following assumption:

Assumption 3.11**.**

$\Delta\mathrm{VIX}^{2}_{T}$ is log-normal.

In fact, this is almost the same as the assumption by Bayer, Friz and Gatheral [3]; however, they did not compute the variance exactly as in Proposition 3.8, and instead considered the lognormal approximation

Assumption 3.12** (Bayer-Friz-Gatheral).**

$\log(\Delta\mathrm{VIX}^{2}_{T})$ is Gaussian with mean $\tilde{\mu}$ and variance $\tilde{\sigma}^{2}$ given by

[TABLE]

Lemma 3.13**.**

A VIX future is worth

[TABLE]

Proof.

Since, by assumption $\Delta^{1/2}\mathrm{VIX}_{T}\,{\buildrel\Delta\over{=}}\,\exp\left(\frac{\mu+\sigma\mathcal{N}(0,1)}{2}\right)$, the price of a VIX future directly reads

[TABLE]

and the second case follows analogously. ∎

As opposed to the simulation schemes, the log-normal approximation does depend on the qualitative properties of the initial forward variance curve $\xi_{0}(\cdot)$. Hence, different curves should be analysed to check the robustness of the method. We now exploit the approximation in Assumption 3.11 to obtain a closed-form Black-Scholes type formulae for European options on the VIX.

Lemma 3.14**.**

For $0\leq t\leq T$, let $\mathfrak{V}_{T}(t):=\mathbb{E}\left(\mathrm{VIX}_{T}|\mathcal{F}_{t}\right)$ denote the price at time $t$ of a VIX future maturing at $T$. Then, under Assumption 3.11, a European Call option on a VIX future maturing at $T$ is worth

[TABLE]

where $\widetilde{K}:=[\log(K^{2}\Delta)-\log\int_{T}^{T+\Delta}\xi_{0}(t)\mathrm{d}t+\sigma^{2}/2]/\sigma$, $d_{1}:=-\widetilde{K}+\frac{1}{2}\sigma$ and $d_{2}:=-\widetilde{K}$, with $\sigma^{2}$ as in Proposition 3.8

Proof.

Under Assumption 3.11, the lemma follows directly from the following trivial computations:

[TABLE]

∎

3.4.1. Numerical tests of the log-normal approximation

We perform a numerical analysis of the approximation by pricing VIX Futures using the parameters in (3.5). We consider $10^{6}$ simulations and the three qualitative scenarios introduced in (3.4) for the forward variance curve.

Figures 5 suggest that the approximation is accurate being the difference of order $10^{-3}$ or less. The oscillating nature of the difference is also a good sign since it is probably caused by the Monte Carlo error and does not show any monotonicity. Furthermore, the method shows to be robust for different type of curves $\xi_{0}(\cdot)$. Moreover, the log-normal approximation seems to converge to the true mean, avoiding the oscillations of the Truncated Cholesky method. Therefore, we may conclude that the log-normal approximation produces a good output, with desired smoothness properties. On the other hand, we also analyse European VIX Call options repeating the previous parameters and considering the flat forward variance curve from Scenario $1$.

Figure  6 shows that the log-normal approximation is accurate for at-the-money options. Nevertheless, the difference increases with maturity, and hence one must be careful when using this approximation to price options with long maturities. However, in practice, on Equity markets, liquid maturities are only up to two years.

3.4.2. Numerical tests

We benchmark Bayer-Friz-Gatheral’s approximation against ours with the parameters in (3.5). Figures 7 and 8 suggest that both are similar. In particular, we observe that the variance deviates as maturity increases. Nevertheless, for practical purposes, both approximations agree over a four-year horizon, long enough to cover available Futures data. Many functional forms of $\xi_{0}(t)$ were also tested giving similar results as the ones shown in Figures 7 and 8. We conclude that the approximation by Bayer, Friz and Gatheral is accurate enough for practical purposes and yields and reduces significantly the computational costs since the computation of $\tilde{\sigma}^{2}$ involves a single integral while our approximation requires a double integral (for $\sigma^{2}$). In particular, in order to generate Figure 7 our approximation was 40 times slower than Bayer, Friz and Gatheral’s.

3.5. VIX Futures calibration in rBergormi

In light of the promising results obtained in the previous sections, we create a calibration algorithm based on the log-normal approximation by Bayer, Friz and Gatheral [3]. Even if our approximation seems to be more accurate the computational cost is much larger and the difference between both approaches has been shown to be very small. Moreover, the approach by Bayer, Friz and Gatheral allows to compute the gradient of the objective function in a semi-explicit form, which in terms of optimisation and calibration is extremely useful.

3.5.1. Objective function

To calibrate the model to VIX Futures, we first define the objective function

[TABLE]

which we minimise over $(\nu,H)$. Here $(\mathfrak{F}_{i})_{i=1,\ldots,N}$ are the observed Futures prices on the time grid $T_{1}<\ldots<T_{N}$, $\mathfrak{V}_{T_{i}}=\Delta^{-1/2}\sqrt{\int_{T_{i}}^{T_{i}+\Delta}\xi_{0}(t)\mathrm{d}t}\exp\left(-\frac{\tilde{\sigma}_{i}^{2}}{8}\right)$. and

[TABLE]

which is the closed-form expression of the variance in Assumption 3.12. The gradient of the objective function is an important source of information in many optimisation algorithms. To compute it, we differentiate the objective function in (3.8) with respect to $\nu$ and $H$ and apply the chain rule:

[TABLE]

where

[TABLE]

On the other hand, $\displaystyle\frac{\partial\mathcal{L}^{\mathfrak{F}}}{\partial H}(\nu,H)=-\frac{1}{4}\sum_{i=1}^{N}(\mathfrak{V}_{T_{i}}-\mathfrak{F}_{i})\mathfrak{V}_{T_{i}}\frac{\partial\tilde{\sigma}_{i}^{2}}{\partial H}$, with

[TABLE]

where

[TABLE]

and K_{H}=\frac{H_{+}}{H}\Big{\{}1+H\psi\left(2-H_{+}\right)-H\left(\psi(H_{+})+\psi(1-2H_{-})\right)\Big{\}}, where $\psi$ is the digamma function.

3.5.2. Obtaining the initial forward variance curve

The initial forward variance curve plays a crucial role in both the Bergomi [5] and the rough Bergomi models [3], since it is a market input. In particular it depends on the current term structure of variance swaps. Even if variance swaps are not traded in standard exchange markets, they are actively traded over-the-counter (OTC). This in turn means that there is no observable data and we must establish a valuation method for variance swaps. The celebrated static replication formula by Carr and Madan [8], applicable here since the underlying process is a continuous semimartingale, allows us to price any variance swap. Nevertheless, out-of-money Call and Put option prices are needed for all possible strikes. Since this information is not available in practice, one can either adopt a model-free valuation formula or directly propose a parameterisation for the implied volatility surface. The first approach involves a discretisation of the static replication formula, which is how the VIX index is computed by the Chicago Board Options Exchange. It is not easy, however, to extend this computation for large maturities (VIX is a $30$-day ahead index), where liquidity of options may play a major role. On the other hand, the second approach allows to calibrate a model using available data and additionally allows to interpolate / extrapolate available option data to all strikes and maturities. In this work we follow the latter approach, with the eSSVI parameterisation [18] for the implied volatility surface, which is a refinement of the SSVI parametrisation introduced in [16]:

[TABLE]

where $\theta_{t}$ is the observed ATM variance curve, and where the shape function $\varphi(\cdot)$ takes the form $\varphi(\theta)=\eta\theta^{-\lambda}(1+\theta)^{\lambda-1}$. For the correlation parameter $\rho(\cdot)$ we restrict it to the following functional form:

[TABLE]

ensuring that $|\rho(\cdot)|\leq 1$. We shall indistinctly refer to the total implied variance by $\sigma_{\mathrm{BS}}^{2}(\cdot)t$ or $w(\cdot)$, where $\sigma_{\mathrm{BS}}(\cdot)$ represents the implied volatility. Gatheral and Jacquier [16] found (sufficient and almost necessary) conditions on the parameters $\rho$, $\varphi(\cdot)$, $\eta$ and $\lambda$ preventing arbitrage. In the eSSVI formulation, the correlation has a term structure (3.10), and care must be taken in order to preclude arbitrage. Concretely, following  [18], the restriction

[TABLE]

where $\gamma:=\partial_{\theta}(\theta\varphi(\theta))/\varphi(\theta)$, is a necessary condition to preclude calendar spread arbitrage. To prevent butterfly arbitrage, exactly as in the SSVI parameterisation, it is sufficient [16, Theorem 4.2] to check that the inequality $\theta_{t}\varphi^{2}(\theta_{t})(1+|\rho(\theta_{t})|)\leq 4$ holds for all maturity $t$. In this parameterisation, variance swaps can be computed in closed form, as proved by Martini [21], based on earlier works by Gatheral [13] and Fukasawa [12]:

Proposition 3.15**.**

The fair strike (in total variance) of a variance swap in the eSSVI model reads

[TABLE]

where $\chi_{t}:=\frac{1}{4}[1-\rho(\theta_{t})^{2}]\theta_{t}\varphi(\theta_{t})$, and

[TABLE]

Recalling the relation between variance swaps and the forward variance curve, we have

[TABLE]

Remark 3.16**.**

In order to interpolate/extrapolate the eSSVI, it is necessary to also interpolate/extrapolate $\theta_{t}$ for all $t$. However, $\theta_{t}$ is only observed on a discrete set of maturity dates, and consequently, cubic splines are used to interpolate/extrapolate all other maturities.

3.5.3. Calibration algorithm and numerical results

We first calibrate the eSSVI parameterisation (3.9) on the SPX implied volatility surface on the $4$th of December $2015$. Figure 9 shows the fit for the shortest, medium and longest maturities available in the data set. For short maturities the eSSVI is not able to fully capture the volatility smile, however as maturity increases the fit improves remarkably.

For VIX Futures, the calibration algorithm reads as follows:

Algorithm 3.17** (VIX Futures calibration algorithm in the rough Bergomi model).**

1. (i)

Calibrate eSSVI to available SPX option data; 2. (ii)

compute the variance swap term structure $(\sigma_{0}(t)^{2})_{t\geq 0}$ using Proposition 3.15; 3. (iii)

extract the initial forward variance curve, $\xi_{0}(\cdot)$ via $\xi_{0}(t)\approx\sigma^{2}_{0}(t)+\frac{\sigma^{2}_{0}(t+\varepsilon)-\sigma^{2}_{0}(t-\varepsilon)}{2\varepsilon}t$ (with $\varepsilon=1E-8$); 4. (iv)

minimise (over $\nu,H$) the objective function in (3.8).

Figures 10-12 suggest that the model fits very well the observed VIX Futures term structure for different dates. Moreover, we notice that both the model and the observed data are qualitatively equal in terms of convexity/concavity. In the rough Bergomi model this information is obtained from option prices through $\xi_{0}$, which suggests a correspondence in the market between VIX futures and SPX options. However, we also observe that in all three cases the error is greater for short maturities, mimicking the calibration limits of eSSVI for short maturities, as detailed in Section 3.5.2.

Remark 3.18**.**

The reader should recall the importance of the initial forward variance curve $\left(\xi_{0}(t)\right)_{t\geq 0}$ in the VIX Futures process, since $(\mathfrak{V}_{t})_{t\geq 0}$ depends on the whole path of $\xi_{0}$ up to time $t$. Therefore, even if $\xi_{0}$ is only misspecified for short maturities (as is the case of the eSSVI), this affects the whole term structure of the VIX process (for details we refer the reader to Section 3). Therefore, an improved $\xi_{0}$ estimation would not only increase the accuracy of the model for short maturities, but also for the whole term structure.

Remark 3.19**.**

Our calibration involves two different data sets, SPX options and VIX Futures: the Vanilla quotes are extracted from the CBOE delayed option quotes page222CBOE delayed option quotes: http://www.cboe.com/delayedquote/quotetable.aspx and the VIX Futures from the CBOE VIX futures historical data page333CBOE VIX futures historical data: http://cfe.cboe.com/data/historicaldata.aspx. We perform an aggregation of the different Future quotes, since they are quoted on a Future per Future basis, and check the consistency between the two data sets by comparing the left end extrapolation of the Futures curve with the theoretical VIX computed from option prices.

4. From VIX Futures to SPX options

In this final chapter we assess whether the Hurst parameter $H$ obtained through the VIX Futures calibration algorithm is consistent with SPX options. For this purpose, we calibrate the rough Bergomi model to SPX option data by fixing the parameter $H$ and letting the algorithm calibrate $\nu$ and $\rho$. One of the main reasons to fix $H$ is that the hybrid scheme introduced in Section 2.1 remarkably reduces its complexity to $\mathcal{O}(n)$, since the $\mathcal{O}(n\log n)$ complexity of the Volterra is computed only once and reused afterwards. Therefore, by fixing $H$ the pricing scheme is much faster when several valuations are performed, as is the case of a calibration algorithm.

4.1. Pricing in the rough Bergomi model

We present a pricing scheme, where the Volterra process $\mathcal{V}$ is simulated using a hybrid scheme, while a standard Euler scheme generates the paths of the stock process:

Algorithm 4.1** (Simulation of the rough Bergomi model).**

Consider the grid $\mathcal{T}:=\{t_{i}\}_{i=0,\ldots,n_{T}}$, and fix $\kappa\geq 1$.

1. (i)

Simulate the Volterra process $\mathcal{V}$ on the grid $\mathcal{T}$ using the hybrid scheme; 2. (ii)

simulate the variance process as $V_{t}=\xi_{0}(t)\mathcal{E}(2\nu C_{H}\mathcal{V}_{t})$, for $t\in\mathcal{T}$ and where $([\mathcal{V}]_{t})_{t\geq 0}$ is given in (2.7); 3. (iii)

extract the path of the Brownian motion $Z$ driving $\mathcal{V}$:

[TABLE]

compute $\{Z^{\perp}\}_{i=0}^{n_{T}-1}$ where $Z^{\perp}_{i}\,{\buildrel\Delta\over{=}}\,\mathcal{N}(0,1/n_{T})$ is an independent standard Gaussian sample; 4. (iv)

correlate the two Brownian motions via $W_{t_{i}}-W_{t_{i-1}}=\rho\overline{Z}_{i-1}+\sqrt{1-\rho^{2}}\overline{Z}_{i-1}^{\perp};$ 5. (v)

simulate $S_{t_{i}}=\exp(X_{t_{i}})$ using a forward Euler scheme:

[TABLE] 6. (vi)

compute the expectation by averaging the payoff over all terminal values of each path.

4.2. Calibration of SPX options via VIX Futures

We first follow the calibration algorithm 3.17 to obtain $H$ and $\xi_{0}$, and we then aim at minimising, over $(\nu,\rho)$, the objective function

[TABLE]

where $\mathrm{C}_{T_{i,j}}$ is the Call price given by the rough Bergomi model, computed using the scheme introduced in Section 4.1, with maturity $T_{i}$ and strike $K^{(j)}$. On the other hand, $(\mathrm{C}^{\mathrm{obs}}_{i,j})_{i,j}$ is the set observed Call prices in the time grid $T_{1}<\ldots<T_{N}$ and strike grid $K^{(1)}<\ldots<K^{(L)}$. In order to optimise the calibration algorithm, we first compute the Volterra process $\mathcal{V}$ , which will then be used in a forward Euler simulation at each calibration step:

Algorithm 4.2** (Calibration algorithm for SPX options via VIX Futures).**

1. (i)

Calibrate $H$ and $\xi_{0}$ using the VIX Futures; 2. (ii)

compute $M$ paths of the Volterra process, $\{\mathcal{V}^{(u)}\}_{u=1}^{M}$ and extract the Brownian motions $\{Z^{(u)}\}_{u=1}^{M}$ driving each process. Also, compute independent Brownian motions $\{Z^{\perp(u)}\}_{u=1}^{M}$; 3. (iii)

evaluate the Call prices in each calibration step:

[TABLE] 4. (iv)

compute the Call price for each available maturity $\{T_{1},\ldots,T_{N}\}$ and set of strikes $\{K^{(1)},\ldots,K^{(L)}\}$:

[TABLE] 5. (v)

minimise over $(\nu,\rho)$ the objective function $\mathcal{L}^{\mathrm{C}}(\nu,\rho)$ in (4.1).

Remark 4.3**.**

Item (v) in the algorithm above may change the optimal values for $\nu$, which was initially calibrated in (i) from the VIX Futures. Backtesting however shows that the calibration in (i) is not really affected by this.

4.3. Results

We calibrate the model on December $4$, $2015$, fixing $H=0.09237$ obtained previously through VIX, and plot the fit in Figure 13. The model is not fully consistent for short maturities, which may follow from the inability of $\xi_{0}$ to fully capture the smiles for these maturities, but the fit greatly improves with maturity. Interestingly, we observe a $20\%$ difference between the the parameter $\nu$ obtained through VIX calibration and the one obtained through SPX. This suggests that the volatility of volatility in the SPX market is $20\%$ higher when compared to VIX, revealing potential data inconsistencies (arbitrage?). Nevertheless, we emphasise the importance of an accurate $\xi_{0}$ curve to improve the fit to SPX and to provide an efficient joint calibration.

Conclusion

Following the path set by Bayer, Friz and Gatheral [3], we developed here a relatively fast algorithm to calibrate VIX Futures and the VIX smile, consistently with the SPX smile, in the rough Bergomi model. The clear strength of this model is that only a few parameters are needed, making the (re)calibration robust and stable. From a trader’s point of view, we highlight some potential market discrepancy between the VIX and the SPX, and leave a refined analysis thereof for future research.

Appendix A The hybrid scheme

We briefly recall the hybrid scheme developed in [4]. Following the notation in Definition 2.1, we consider a (truncated) Brownian semistationary process $\mathcal{B}(\alpha,W)$, and introduce the truncation parameter $\kappa\in\mathbb{N}$. On an equidistant grid $\mathcal{T}:=\{t_{i}=i/n\}_{i=0,\ldots,n_{T}}$, with $n_{T}:=\lfloor nT\rfloor$, for $n\geq 2$, under Definition 2.1, the hybrid scheme for the BSS process $\mathcal{B}$ is approximated by $\mathcal{B}_{n}(t_{i})=\widetilde{\mathcal{B}}_{n}(i)+\widehat{\mathcal{B}}_{n}(i)$ with

[TABLE]

with $L_{g}$ introduced in Definition 2.1, and where

[TABLE]

For any $i,k$, the random variables $\overline{W}_{i}$ and $\overline{W}_{i,k}$ are centred Gaussian with the following covariance structure:

[TABLE]