Heston Stochastic Vol-of-Vol Model for Joint Calibration of VIX and S&P 500 Options

TL;DR

This paper introduces a stochastic volatility-of-volatility extension to the Heston model, enabling efficient joint calibration of S&P 500 options and VIX, with practical approximation formulas derived via perturbation techniques.

Contribution

It proposes a new stochastic volatility-of-volatility model based on the Heston framework, with a first-order approximation for option pricing that simplifies calibration.

Findings

Model accurately fits S&P 500 and VIX data

Approximation formulas are computationally efficient

Calibration process is streamlined and effective

Abstract

A parsimonious generalization of the Heston model is proposed where the volatility-of-volatility is assumed to be stochastic. We follow the perturbation technique of Fouque et al (2011, CUP) to derive a first order approximation of the price of options on a stock and its volatility index. This approximation is given by Heston's quasi-closed formula and some of its Greeks. It can be very efficiently calculated since it requires to compute only Fourier integrals and the solution of simple ODE systems. We exemplify the calibration of the model with S&P 500 and VIX data.

Equations

$$\text{mbox}{VIX}_t^2\approx \frac{2}{{\tau }_0}{\int }_0^{+\infty }Q(t,t+{\tau }_0,K)\frac{dK}{K^2},$$

$$\displaystyle\left\{\begin{array}[]{l}dS_{t}=(r-q)S_{t}dt+\sqrt{V_{t}}S_{t}dW_{t}^{S},\\ \\ dV_{t}=\kappa(m-V_{t})dt+\eta_{t}\sqrt{V_{t}}dW_{t}^{V},\\ \\ dW_{t}^{S}dW_{t}^{V}=\raisebox{1.5pt}{$\rho$}_{SV}dt\end{array}\right.$$

$$\text{mbox}{VIX}_t^2=\mathbb{E}\left[\frac{1}{{\tau }_0}{\int }_t^{t+{\tau }_0}{V}_{u}du|{\mathcal{F}}_t\right],$$

$$\displaystyle\left\{\begin{array}[]{l}\eta_{t}=\eta(Y^{\varepsilon}_{t},Z^{\delta}_{t})\\ \\ \displaystyle dY^{\varepsilon}_{t}=\frac{V_{t}}{\varepsilon}\alpha(Y_{t}^{\varepsilon})dt+\sqrt{\frac{V_{t}}{\varepsilon}}\beta(Y_{t}^{\varepsilon})dW_{t}^{Y},\\ \\ dZ^{\delta}_{t}=V_{t}\delta c(Z^{\delta}_{t})dt+\sqrt{\delta V_{t}}g(Z^{\delta}_{t})dW_{t}^{Z},\end{array}\right.$$

$$d{W}_t^{i}d{W}_t^j=\rho dt,i,j=S,V,Y,Z.$$

$$2\kappa m\ge \eta (y,z)\ge c>0,$$

$$P_S^{\epsilon ,\delta }(t,s,v,y,z)$$

$$P_V^{\epsilon ,\delta }(t,v,y,z)$$

$$|g^{\epsilon ,\delta }-f^{\epsilon ,\delta }|\le C(\epsilon +\delta ),$$

$$g^{\epsilon ,\delta }-f^{\epsilon ,\delta }=O(\epsilon +\delta ).$$

$$P_S^{\epsilon ,\delta }(t,s,v,y,z)$$

$$P_V^{\epsilon ,\delta }(t,v,y,z)$$

$$P_{S_0}+P_{S_{1,0}}^{\epsilon }+P_{S_{0,1}}^{\delta }=\frac{e^{-r\tau }}{\pi }{\int }_0^{+\infty }Re\left(e^{-i\xi x(t,s)}(1+h_{S_0}^{\epsilon ,\delta }+v{h}_{S_1}^{\epsilon ,\delta }+v^2{h}_{S_2}^{\epsilon ,\delta })G_S(\tau ,\xi ,v,z)\hat{\phi }(\xi )\right)d{\xi }_r,$$

$$\tau (t)=T-t,$$

$$x(t,s)=(r-q)(T-t)+\log s,$$

$$\hat{\phi }(\xi )={\int }_{\mathbb{R}}\phi (e^{x-(r-q)\tau })e^{-i\xi x}dx,\xi ={\xi }_r+i{\xi }_i,$$

$$G_S(\tau ,\xi ,v,z)=e^{C(\tau ,\xi ,z)+vD(\tau ,\xi ,z)},$$

$$C(\tau ,\xi ,z)=\frac{\kappa m}{{\overline{\eta }}^2(z)}\left((\kappa +i\overline{\rho }(z)\overline{\eta }(z)\xi -d(\xi ,z))\tau -2\log \left(\frac{e^{-d(\xi ,z)\tau }/g(\xi ,z)-1}{1/g(\xi ,z)-1}\right)\right),$$

$$D(\tau ,\xi ,z)=\frac{\kappa +i\overline{\rho }(z)\overline{\eta }(z)\xi +d(\xi ,z)}{{\overline{\eta }}^2(z)}\left(\frac{1-e^{d(\xi ,z)\tau }}{1-g(\xi ,z)e^{d(\xi ,z)\tau }}\right),$$

$$g(\xi ,z)=\frac{\kappa +i\overline{\rho }(z)\overline{\eta }(z)\xi +d(\xi ,z)}{\kappa +i\overline{\rho }(z)\overline{\eta }(z)\xi -d(\xi ,z)},$$

$$d(\xi ,z)=\sqrt{{\overline{\eta }}^2(z)({\xi }^2-i\xi )+(\kappa +i\overline{\rho }(z)\overline{\eta }(z)\xi {)}^2},$$

$$h_{S_0}^{\epsilon ,\delta }=f_0^{\epsilon }+g_0^{\delta },h_{S_1}^{\epsilon ,\delta }=f_1^{\epsilon }+g_1^{\delta },h_{S_2}^{\epsilon ,\delta }=g_2^{\delta },$$

$$\displaystyle\left\{\begin{array}[]{l}\displaystyle\frac{\partial f^{\varepsilon}_{1}}{\partial\tau}(\tau,\xi,z)=(\overline{\eta}^{2}(z)D(\tau,\xi,z)-(\kappa+\overline{\rho}(z)\overline{\eta}(z)i\xi))f^{\varepsilon}_{1}(\tau,\xi,z)\\ \\ \displaystyle\hskip 28.45274pt-i\xi V_{1,2}^{\varepsilon}(z)D^{2}(\tau,\xi,z)-\xi^{2}V_{2,1}^{\varepsilon}(z)D(\tau,\xi,z)+V_{0,3}^{\varepsilon}(z)D^{3}(\tau,\xi,z),\\ \\ \displaystyle\frac{\partial f^{\varepsilon}_{0}}{\partial\tau}(\tau,\xi,z)=\kappa mf^{\varepsilon}_{1}(\tau,\xi,z),\\ \\ f^{\varepsilon}_{0}(0,\xi,z)=f^{\varepsilon}_{1}(0,\xi,z)=0,\end{array}\right.$$

$$\displaystyle\left\{\begin{array}[]{l}\displaystyle\frac{\partial g^{\delta}_{2}}{\partial\tau}(\tau,\xi,z)=-2(\kappa+\overline{\rho}(z)\overline{\eta}(z)i\xi-\overline{\eta}^{2}(z)D(\tau,\xi,z))g^{\delta}_{2}(\tau,\xi,z)\\ \\ \hskip 28.45274pt\displaystyle(V_{0,1,\eta}^{\delta}(z)-i\xi V_{1,0,\eta}^{\delta}(z))\frac{\partial D}{\partial\eta}+(V_{0,1,\rho}^{\delta}(z)\ -i\xi V_{1,0,\rho}^{\delta}(z))\frac{\partial D}{\partial\rho},\\ \\ \displaystyle\frac{\partial g^{\delta}_{1}}{\partial\tau}(\tau,\xi,z)=-(\kappa+\overline{\rho}(z)\overline{\eta}(z)i\xi-\overline{\eta}^{2}(z)D(\tau,\xi,z))g^{\delta}_{1}(\tau,\xi,z)\\ \\ \displaystyle\hskip 28.45274pt+\overline{\eta}^{2}(z)g^{\delta}_{2}(\tau,\xi,z)+(V_{0,1,\eta}^{\delta}(z)-i\xi V_{1,0,\eta}^{\delta}(z))\frac{\partial C}{\partial\eta}+(V_{0,1,\rho}^{\delta}(z)-i\xi V_{1,0,\rho}^{\delta}(z))\frac{\partial\overline{C}}{\partial\rho},\\ \\ \displaystyle\hskip 28.45274pt+V_{0,1,\eta}^{\delta}(z)\frac{\partial D}{\partial\eta}+V_{0,1,\rho}^{\delta}(z)\frac{\partial D}{\partial\rho}\\ \\ \displaystyle\frac{\partial g^{\delta}_{0}}{\partial\tau}(\tau,\xi,z)=\kappa mg^{\delta}_{1}(\tau,\xi,z),\\ \\ g^{\delta}_{0}(0,\xi,z)=g^{\delta}_{1}(0,\xi,z)=g^{\delta}_{2}(0,\xi,z)=0.\end{array}\right.$$

$$\overline{\eta }(z)=\sqrt{⟨{\eta }^2(\cdot ,z)⟩},$$

$$\overline{\rho }(z)=\rho \frac{⟨\eta (\cdot ,z)⟩}{\overline{\eta }(z)},$$

$$V_{1,2}^{\epsilon }(z)=-\sqrt{\epsilon }\frac{\rho }{2}⟨\beta \frac{\partial \varphi }{\partial y}(\cdot ,z)⟩-\sqrt{\epsilon }\rho \rho ⟨\eta (\cdot ,z)\beta \frac{\partial \psi }{\partial y}(\cdot ,z)⟩,$$

$$V_{2,1}^{\epsilon }(z)=-\sqrt{\epsilon }\rho \rho ⟨\beta \frac{\partial \psi }{\partial y}(\cdot ,z)⟩,$$

$$V_{0,3}^{\epsilon }(z)=-\sqrt{\epsilon }\frac{\rho }{2}⟨\eta (\cdot ,z)\beta \frac{\partial \varphi }{\partial y}(\cdot ,z)⟩,$$

$$V_{0,1,\eta }^{\delta }(z)=\sqrt{\delta }\rho g(z)⟨\eta (\cdot ,z)⟩{\overline{\eta }}^{\prime }(z),$$

$$V_{0,1,\rho }^{\delta }(z)=\sqrt{\delta }\rho g(z)⟨\eta (\cdot ,z)⟩{\overline{\rho }}^{\prime }(z),$$

$$V_{1,0,\eta }^{\delta }(z)=\sqrt{\delta }\rho g(z){\overline{\eta }}^{\prime }(z),$$

$$V_{1,0,\rho }^{\delta }(z)=\sqrt{\delta }\rho g(z){\overline{\rho }}^{\prime }(z).$$

$$P_{V_0}+P_{V_{1,0}}^{\epsilon }+P_{V_{0,1}}^{\delta }=\frac{e^{-r\tau }}{\pi }{\int }_0^{+\infty }Re\left((1+h_{V_0}^{\epsilon ,\delta }+v{h}_{V_1}^{\epsilon ,\delta }+v^2{h}_{V_2}^{\epsilon ,\delta })G_V(t,v,\nu )\hat{\phi }(\nu )\right)d{\nu }_i,$$

$$\hat{\phi }(\nu )={\int }_0^{+\infty }{e}^{i\nu v}\phi (\gamma (v))dv,\nu ={\nu }_r+i{\nu }_i,$$

$$\gamma (v)=\sqrt{m(1-\theta )+\theta v},$$

$$\theta =\frac{1-e^{-\kappa {\tau }_0}}{\kappa {\tau }_0},{\tau }_0=30/360,$$

$$G_V(t,v,\nu ,z)=e^{A(\tau ,\nu ,z)+vB(\tau ,\nu ,z)},$$

$$A(\tau ,\nu ,z)=-\frac{2\kappa m}{{\overline{\eta }}^2(z)}\log \left(\nu \frac{{\overline{\eta }}^2(z)}{2\kappa }(1-e^{-\kappa \tau })+1\right),$$

$$B(\tau ,\nu ,z)=\frac{\nu e^{-\kappa \tau }}{\nu \frac{{\overline{\eta }}^2(z)}{2\kappa }(1-e^{-\kappa \tau })+1}.$$

$$\displaystyle\left\{\begin{array}[]{l}\displaystyle\frac{\partial h_{V_{2}}^{\varepsilon,\delta}}{\partial\tau}(\tau,\nu,z)=2(-\kappa+B(\tau,\nu)\overline{\eta}^{2}(z))h_{V_{2}}^{\varepsilon,\delta}(\tau,\nu,z)+V_{1}^{\delta}(z)B(\tau,\nu)\frac{\partial B}{\partial\eta}(\tau,\nu),\\ \\ \displaystyle\frac{\partial h_{V_{1}}^{\varepsilon,\delta}}{\partial\tau}(\tau,\nu,z)=(-\kappa+B(\tau,\nu)\overline{\eta}^{2}(z))h_{V_{1}}^{\varepsilon,\delta}(\tau,\nu,z)+(2\kappa m+\overline{\eta}^{2}(z))h_{V_{2}}^{\varepsilon,\delta}(\tau,\nu,z)\\ \\ \displaystyle\hskip 28.45274pt+V_{3}^{\varepsilon}(z)B^{3}(\tau,\nu)+V_{1}^{\delta}(z)\left(\frac{\partial B}{\partial\eta}(\tau,\nu)+B(\tau,\nu)\frac{\partial A}{\partial\eta}(\tau,\nu)\right),\\ \\ \displaystyle\frac{\partial h_{V_{0}}^{\varepsilon,\delta}}{\partial\tau}(\tau,\nu,z)=\kappa mh_{V_{1}}^{\varepsilon,\delta}(\tau,\nu,z),\\ \\ h_{V_{0}}^{\varepsilon,\delta}(0,\nu,z)=h_{V_{1}}^{\varepsilon,\delta}(0,\nu,z)=h_{V_{2}}^{\varepsilon,\delta}(0,\nu,z)=0.\end{array}\right.$$

$$V_3^{\epsilon }(z)$$

Full text

1 Introduction

The volatility index of the S&P 500, acronymed VIX and also known as the fear index, has drawn the attention of researchers and practitioners alike since its first introduction in the US market in 1993, see CBOE (2003) for its current definition. In 2004, future contracts on VIX began to trade at CBOE Futures Exchange and later on, in 2006, options on VIX were firstly negotiated.

From its definition, the VIX index is computed using the price of liquid options on the S&P 500. In fact,

[TABLE]

where $\tau_{0}$ is 30 days and $Q(t,T,K)$ denotes the price of the out-the-money option at time $t$ with maturity $T$ and strike $K$. The approximation sign $\approx$ appears in the equation above because, obviously, the index is computed by discretizing the integral on the left-hand side. Moreover, the square of the VIX is linearly interpolated between the two closest maturities in order to have its value 30 days from $t$.

Hence, the implied volatility surfaces of the S&P 500 and VIX are highly connected, and as a consequence this dependence is very complex. There are few models proposed to solve this calibration issue, see, for instance, Baldeaux and Badran (2014); Carr and Madan (2014); Pacati et al. (2015); Papanicolaou and Sircar (2013) and Cont and Kokholm (2013). The common aspect of the models described in these references is the presence of jumps in the stock price and/or its spot volatility. Differently from these models, we consider here a continuous diffusion model. More precisely, we propose a simple generalization of Heston (1993) where volatility-of-volatility is stochastic.

Furthermore, to the best of our knowdelege, the only continuous models proposed to joint calibrate stock and volatility options appeared in Gatheral (2008). While the calibration of these models rely on Monte Carlo or PDE methods, ours grants us quasi-closed formulas for the first-order approximation of option prices on both markets. Additionally, in the direction of model-free results, there is the work De Marco and Henry-Labordère (2015).

Our approach is based on the multiscale stochastic volatility perturbation technique proposed by Fouque, Papanicolaou, Sircar and Sølna, see Fouque et al. (2011). This approach allows us to approximate option prices under the full model by their prices and Greeks under a simpler model. This method is very flexible and can be adapted to a large number of models and options, as one will be able to see in this paper.

A different perturbation technique was applied in the context of joint calibration in Papanicolaou and Sircar (2013). In this paper, the authors proposed a regime-switching generalization of the Heston model. The perturbation is done in the jump-process that brings the regime-switching feature into the model. This is fundamentally different from the solution proposed here.

Additionally, the Heston model was also generalized in the lines of the multiscale stochastic volatility modelling in Fouque and Lorig (2011). This generalization is also very different from the one pursued here. The main issue is that the simple formula found in our model, shown in Equation (6.58), is not verified in the aforesaid model.

The paper is organized as follows: we describe our model in Section 2 and the main results are stated in Section 3. We discuss the calibration of the proposed model and exemplify it in Section 4. Some generalizations of our model are outlined in Section 5. Finally, the rationale and computations that justify our first-order approximation are shown in Section 6 and in the Appendices A and B.

2 The Model

We will assume that the stock price $S$, under a risk-neutral probability, follows a Heston dynamics with a stochastic volatility-of-volatility (vol-vol):

[TABLE]

where $\eta_{t}$ will be specified later in this section. We will denominate this model by Heston SVV model.

The volatility index of $S$ at $t$, which it will be denoted, because of obvious reasons, by $\mbox{VIX}_{t}$, is defined as

[TABLE]

where $\tau_{0}=30/360$, i.e. 30 calendar days. The main example to have in mind is the S&P 500 and the VIX. The expected value above, as all the other expected values in this work, is under the chosen risk-neutral measure. This risk-neutral measure is taken to match the vanilla option prices for both the stock and its volatility index markets.

Remark 2.1**.**

As it was discussed in the introduction, the VIX is computed by discretizing Equation (1.1). However, under any continuous model with spot variance $V$, the left-hand side of Equation (1.1) can be written as Equation (2.7). Hence, we are actually incurring in a small discretization error when using the non-discretized version of the volatility index.

Let us now specify the particular formula for the process $\eta_{t}$. In order to be able to find a computationally efficient approximation for the price of options on $S$ and VIX, we choose $\eta_{t}$ to be governed by fast and slow time scales. More precisely,

[TABLE]

where $(W_{t}^{S},W_{t}^{V},W_{t}^{Y},W_{t}^{Z})$ is a correlated Brownian motion with

[TABLE]

Assumption 2.2**.**

The assumptions of this model are:

• •

there exists a unique strong solution of the stochastic differential equations (2.6)-(2.13) for fixed $(\varepsilon,\delta)$;

• •

the covariance matrix of $(W_{t}^{S},W_{t}^{V},W_{t}^{Y},W_{t}^{Z})$ is positive-definite;

• •

$\alpha$ and $\beta$ are such that the process $Y^{1}$ has a unique invariant distribution and is mean-reverting as in (Fouque et al., 2011, Section 3.2);

• •

$\eta(y,z)$ is a positive function, smooth in $z$ and such that $\eta^{2}(\cdot,z)$ is integrable with respect to the invariant distribution of $Y^{1}$.

Remark 2.3** (Feller Condition).**

In order to guarantee that $V_{t}>0$ a.s. for all $t\in[0,T]$, one needs to assume

[TABLE]

for every $(y,z)$ and a constant $c>0$. See, for example, Benhamou et al. (2010).

3 Main Results

In this section we present the first-order approximation for the price of derivative contracts on the stock price, $S$, and on its volatility index, VIX. The derivation of these results and formulas will be fully developed in the sections to follow.

We start by fixing two European derivatives with maturity $T$ and payoff functions \raisebox{1.5pt}{\varphi}_{S} and \raisebox{1.5pt}{\varphi}_{V} that depend only on the terminal values $S_{T}$ and $\mbox{VIX}_{T}$, respectively. The no-arbitrage prices, under the chosen risk-neutral measure, of these derivative contracts are given by the conditional expectations:

[TABLE]

In the case of VIX futures (i.e. \raisebox{1.5pt}{\varphi}_{V}(v)=v), there is no discounting term.

We are interested in jointly calibrating our model to options on $S$ and on VIX. Below, we present the first-order approximation for these option prices.

Remark 3.1**.**

More precisely, we say that a function $g^{\varepsilon,\delta}$ is a first-order approximation to the function $f^{\varepsilon,\delta}$ if

[TABLE]

for some constant $C>0$ and for sufficiently small $\varepsilon,\delta>0$. We use the notation

[TABLE]

We have then the following theorem, see Section 6.3.

Theorem** (Accuracy Theorem).**

Under Assumption 2.2 and if the payoff functions \raisebox{1.5pt}{\varphi}_{S} and \raisebox{1.5pt}{\varphi}_{V} are continuous and piecewise smooth, then

[TABLE]

More importantly, each of the functions on the right-hand side of the equations above can be efficiently computed. Indeed, for options on $S$, we find that $P_{S_{0}}$ is the no-arbitrage price of \raisebox{1.5pt}{\varphi}_{S} under the Heston model with constant vol-vol equals $\overline{\eta}(z)$ and effective correlation equals $\overline{\rho}(z)$. Specifically, we have

[TABLE]

where

[TABLE]

with $f^{\varepsilon}_{0}$ and $f^{\varepsilon}_{1}$ satisfying the ODE system:

[TABLE]

and $g^{\delta}_{0}$, $g^{\delta}_{1}$ and $g^{\delta}_{2}$ satisfying the ODE system:

[TABLE]

The market group parameters $(\overline{\eta}(z),\overline{\rho}(z),V_{1,2}^{\varepsilon}(z),V_{2,1}^{\varepsilon}(z),V_{0,3}^{\varepsilon}(z),V_{1,0,\eta}^{\delta}(z),V_{0,1,\eta}^{\delta}(z)$, $V_{1,0,\rho}^{\delta}(z),V_{0,1,\rho}^{\delta}(z))$ are related to the functions describing the model (2.6)-(2.13) through the equations:

[TABLE]

Now, for options on VIX, we find

[TABLE]

where

[TABLE]

with $h_{V_{0}}^{\varepsilon,\delta}$, $h_{V_{1}}^{\varepsilon,\delta}$ and $h_{V_{2}}^{\varepsilon,\delta}$ satisfying the ODE system:

[TABLE]

Moreover, the group market parameters $(V_{1}^{\delta}(z),V_{3}^{\varepsilon}(z))$ are given by

[TABLE]

4 Calibration

Firstly, we would like to point out that all the parameters related to the first-order approximation of $P_{V}^{\varepsilon,\delta}$, i.e. the market group parameters $(\kappa,m,\overline{\eta}(z),V_{1}^{\delta}(z),V_{3}^{\varepsilon}(z))$, appear in the first-order approximation of $P_{S}^{\varepsilon,\delta}$. Indeed, notice that $V_{0,3}^{\varepsilon}(z)=V_{3}^{\varepsilon}(z)$ and $V_{0,1,\eta}^{\delta}(z)=V_{1}^{\delta}(z)$, see Equations (3.8), (3.13), (3.9) and (3.14), respectively. The market group parameters $(\overline{\rho}(z),V_{1,2}^{\varepsilon}(z),V_{2,1}^{\varepsilon}(z),V_{1,0,\eta}^{\delta}(z),$ $V_{1,0,\rho}^{\delta}(z),V_{0,1,\rho}^{\delta}(z))$ appear exclusively in the first-order approximation of $P_{S}^{\varepsilon,\delta}$.

Assume there are available $M_{S}$ and $M_{V}$ options on $S$ and on VIX, respectively. If we denote all parameters by simply $\Theta$ and the implied volatility of options on $S$ by $\widehat{\sigma}_{S}^{i}$ and of options on VIX by $\widehat{\sigma}_{V}^{i}$, we will consider the following calibration problem:

[TABLE]

The choice of the initial guess for the optimization problem above is important in order to avoid its many local minima. In this paper, we first consider the standard Heston model and calibrate it to the implied volatility seen in the market, by solving the optimization problem (4.1). Then, we use these values and set the $V^{\varepsilon,\delta}$ parameters to zero as the initial guess.

We illustrate in Figures 1 and 2 the effect of the $V^{\varepsilon,\delta}$’s on implied volatilities of the SPX and VIX. We used the parameter values $\kappa=15$, $m=0.04$, $\overline{\eta}(z)=2.0$, $\overline{\rho}(z)=-0.5$, $S_{0}=2000$ and $\mbox{VIX}_{0}=0.15$. Interest and dividend rates were set to zero. Moreover, we considered unusually high values of the $V^{\varepsilon,\delta}$’s to accentuate the impact in the implied volatility.

4.1 Real Data Example

In this section, we follow the calibration procedure outlined above on real data. We consider implied volatility surfaces of the S&P 500 and the VIX on August, 21 of 2015. On this day, the S&P 500 closed at 1970.89, and the VIX at 28.03. In order to compute the corresponding future price of S&P 500 and VIX, we use the Put-Call Parity with ATM options, see Bardgett et al. (2014). Implied volatilities are computed using these future prices. See Figures 4 and 4.

We cleaned the data based on the adjustments described in Bardgett et al. (2014). Namely, we removed implied volatilities:

of in-the-money options; 2. 2.

with moneyness below 75% and above 125%; 3. 3.

with zero traded volume; 4. 4.

with zero open interest. 5. 5.

with maturity larger than 1 year.

In order to understand the improvement of the correction terms of our first-order approximation, we follow the same calibration procedure outlined above for the Heston model (i.e. when the vol-vol is constant).

We present in Figure 5 the implied volatility and its calibrated approximation using the Heston and the Heston SVV models. One maturitiy of approximately 120 days was chosen in order to show the capacity of the model to capture both skews of the S&P 500 and the VIX implied volatilities. Very short maturities are refrained because of the nature of our first-order approximation, see, for instance, Fouque et al. (2011). As explained in Generalization 3 in Section 5, the term-structure of VIX could be captured by allowing time-dependence in some parameters of the model.

The calibrated parameters are shown in Table 2. The reader should notice that the $V^{\varepsilon,\delta}$’s parameters satisfy their basic assumption of being small and that under the Heston SVV model, the vol-of-vol decreases compared to the standard Heston model, a desirable aspect for this model. In Table 2 and Figure 6, one can observe a major improvement of the calibration when comparing the mean squared error, as in Equation (4.1). The enhancement of the calibration is better seen in S&P 500’s implied volatilities. Regarding VIX, Heston SVV model shows a better fit, presenting the concave shape of the curve.

\DTLloaddb

hestonparam_paper_20150821.csv \DTLloaddberrorparam_paper_error_20150821.csv

5 Generalizations

Our model could be generalized in many ways, some simpler than others. Below we briefly discuss some of these possibilities:

Generalization 1**.**

Notice that we could have explicitly considered the market prices of volatility risk as it is done in Fouque et al. (2011), and in doing so we would have an additional term of order $\varepsilon^{-1/2}$ and a term of order $\delta^{1/2}$ in the drifts of $Y^{\varepsilon}$ and $Z^{\delta}$ respectively, both depending on $Y^{\varepsilon}$ and $Z^{\delta}$. They could have been handled in the same way it is done in the aforesaid reference. For simplicity, we do not consider this generalization here.

Generalization 2**.**

It is fairly simple to generalize the model above to deal with time-dependent interest rate and dividend yield. In order to make the exposition clearer, we will not consider it in the derivation of the first-order approximation shown in the sections to follow. However, this feature is present in the numerical computation in Section 4.1.

Generalization 3**.**

In order to capture VIX’s term structure, one would need to introduce time dependent parameters. The most parsimonious choice would be the long-run mean, $m$. This would add a moderate computational difficulty to our model, along with a more cumbersome notation and derivation of the first-order approximation, and therefore outside the scope of the paper. For a similar generalization, we refer the reader to Sepp (2008) and Mikhailov and Nögel (2004). Another approach to deal with VIX’s term structure would be along the lines of the work Fouque et al. (2004).

Generalization 4**.**

It should be straight forward to adapt the machinery developed here to deal with the two-factor model of Christoffersen et al. (2014):

[TABLE]

with a simple correlation structure for the four-dimensional Brownian motion $(W_{t}^{S^{1}},W_{t}^{S^{2}},W_{t}^{V^{1}},W_{t}^{V^{2}})$.

Generalization 5**.**

A more complex generalization would be to consider that the mean-reverting rate, $\kappa$, is itself stochastic depending on $Y^{\varepsilon}$ and $Z^{\delta}$ and that the long-run mean, $m$, satisfies $\kappa(y,z)m(y,z)=a$:

[TABLE]

6 Perturbation Framework

The first-order approximation for option prices on the the stock and its volatility index will be developed next. The arguments shown here justify the results and formulas presented in Section 3. These arguments follow the ideas thoroughly explained in Fouque et al. (2011).

6.1 Options on the Stock

In this section, we derive the first-order approximation of the price of derivative contracts on $S$. Notice that, by the Feynman-Kac’s Formula, $P_{S}^{\varepsilon,\delta}$, that is given in Equation (3.1), satisfies the following PDE

[TABLE]

where the differential operator $\mathcal{L}^{\varepsilon,\delta}_{S}$ is given by

[TABLE]

with

[TABLE]

We now develop the singular and regular perturbation analysis for the option price $P_{S}^{\varepsilon,\delta}$ following the method outlined in Fouque et al. (2011). The reader will readily realize that the derivation is very similar to the Black–Scholes perturbation analysis performed in the aforesaid reference.

We formally write $P_{S}^{\varepsilon,\delta}$ in powers of $\sqrt{\delta}$,

[TABLE]

and then, by Equation (6.4), we choose $P^{\varepsilon}_{S_{0}}$ and $P^{\varepsilon}_{S_{1}}$ to satisfy

[TABLE]

6.1.1 Computing $P_{S_{0}}$

We formally expand $P^{\varepsilon}_{S_{0}}$ in powers of $\sqrt{\varepsilon}$,

[TABLE]

and denote $P_{S_{0,0}}$ simply by $P_{S_{0}}$, where we assume that, at maturity, P_{S_{0}}(T,s,v,y,z)=\raisebox{1.5pt}{\varphi}_{S}(s), $P_{S_{1,0}}(T,s,v,y,z)=0$ and $P_{S_{0,1}}(T,s,v,y,z)=0$. Substituting expansion (6.20) into Equation (6.15), we get the following PDEs:

[TABLE]

with the notation $(i,j)$ denoting the term of $i$th order in $\varepsilon$ and $j$th in $\delta$. Therefore, using the well-known arguments of Fouque et al. (2011), we might choose:

• •

$P_{S_{0}}=P_{S_{0}}(t,s,v,z)$ and $P_{S_{1,0}}=P_{S_{1,0}}(t,s,v,z)$ independent of $y$;

• •

$P_{S_{0}}$ satisfing $\langle\mathcal{L}_{2}^{S}\rangle P_{S_{0}}=0$;

• •

$P_{S_{1,0}}$ solving $\langle\mathcal{L}_{2}^{S}\rangle P_{S_{1,0}}=-v\langle\mathcal{L}_{1}^{S}P_{S_{2,0}}\rangle$.

Define then the Heston differential operator:

[TABLE]

where $\mathcal{L}_{BS}(\sigma)$ is the Black–Scholes differential operator with volatility $\sigma$:

[TABLE]

Hence, by Equation (6.7), \mathcal{L}_{2}^{S}=\mathcal{L}_{H}(\eta(y,z),\raisebox{1.5pt}{\rho}_{SV}). Define now the averaged coefficients

[TABLE]

so that $\langle\mathcal{L}_{2}^{S}\rangle=\mathcal{L}_{H}(\overline{\eta}(z),\overline{\rho}(z))$ and $P_{S_{0}}$ solves the PDE

[TABLE]

The function $P_{S_{0}}$ might be computed using the method developed in Heston (1993), which is precisely described in Appendix A.1.

6.1.2 Computing $P^{\varepsilon}_{S_{1,0}}$

By the 0-order equation (6.23), the following formula holds true

[TABLE]

for some function $c$ that does not depend on $y$. Notice

[TABLE]

Then, denote by $\phi(y,z)$ and $\psi(y,z)$ the solutions of the following Poisson equations

[TABLE]

Hence

[TABLE]

and thus, Equation (6.33) implies

[TABLE]

Therefore, $P_{S_{1,0}}^{\varepsilon}=\sqrt{\varepsilon}P_{S_{1,0}}$ will be chosen to satisfy

[TABLE]

where

[TABLE]

In Appendix A.2, using Fourier transform techniques, we derive a quasi-closed formula for $P^{\varepsilon}_{S_{1,0}}$.

6.1.3 Computing $P^{\delta}_{S_{0,1}}$

We now expand $P_{S_{1}}^{\varepsilon}$ in powers of $\sqrt{\varepsilon}$,

[TABLE]

and then substitute this and the expansion for $P_{S_{0}}^{\varepsilon}$ into Equation (6.19) to find

[TABLE]

Thus, we choose:

• •

$P_{S_{0,1}}=P_{S_{0,1}}(t,s,v,z)$ and $P_{S_{1,1}}=P_{S_{1,1}}(t,s,v,z)$ independent of $y$;

• •

$P_{S_{0,1}}$ satisfying $\langle\mathcal{L}_{2}^{S}\rangle P_{S_{0,1}}=-v\langle\mathcal{M}_{1}^{S}\rangle P_{S_{0}}$.

Notice that

[TABLE]

and therefore,

[TABLE]

where

[TABLE]

In Appendix A.3, using Fourier transform techniques, we find a quasi-closed formula for $P^{\delta}_{S_{0,1}}$.

6.2 Options on the Volatility Index

In this section, we will develop the first-order approximation for the price of options on VIX. Before continuing, it is necessary to study the dynamics of VIX under our model. Observe that, under mild conditions on $\eta(y,z)$, we have the well-known formula

[TABLE]

which implies

[TABLE]

where

[TABLE]

Moreover, we define

[TABLE]

and notice that $\mbox{VIX}_{t}=\gamma(V_{t})$. This implies that we may consider, in this model, derivative contracts on VIX as contracts on $V$ but with a more complicated payoff. In Appendix B, we present the Fourier method presented in Sepp (2008) to compute the price of options on VIX, under constant vol-vol, using this observation.

Additionally, notice that, given $\mbox{VIX}_{t}$, $m$ and $\kappa$, we can find the current value of $V$, $V_{t}$, using Equation (6.58).

Remark 6.1**.**

Under the more complex model described in Generalization 5 in Section 5, $\mbox{VIX}_{t}$ would not be independent of $\varepsilon$ and $\delta$. Indeed, it is fairly easy to show that

[TABLE]

for some functions $F_{i,j}$ that could be explicitly computed. The constant $V_{1,\kappa}^{\delta}(z)$ is related to the fast and slow time scales in $\kappa$. Therefore, in order to compute the first-order approximation for option on VIX, one should consider the same approach as in Fouque et al. (2014). In this paper, the authors used the first-order approximation of future prices and examined the problem of computing the first-order approximation on derivatives on futures as a singular and regular perturbation of an asset whose dynamics itself is only known up to its first-order approximation.

We will now derive the first-order approximation for derivatives contracts on VIX. Define the following differential operator

[TABLE]

where $\mathcal{L}_{0}$, $\mathcal{M}_{2}$ and $\mathcal{M}_{3}$ are the same as in the pricing PDE for derivatives on $S$, see Equations (6.5), (6.9) and (6.10), respectively, and $\mathcal{L}_{1}^{V}$, $\mathcal{L}_{2}^{V}$ and $\mathcal{M}_{1}^{V}$ are given by

[TABLE]

Hence, $P^{\varepsilon,\delta}_{V}$, defined in Equation (3.2), satisfies the following PDE:

[TABLE]

where $\gamma$ is given by Equation (6.60). Following Fouque et al. (2011), as we have done in the previous section, we conclude that the first-order approximation for $P^{\varepsilon,\delta}_{V}$ solves the following PDEs

[TABLE]

where

[TABLE]

We are using the notation $\mathcal{L}_{\!{}_{C\!I\!R}}(\eta)$ because it is related to the infinitesimal generator of a CIR process with constant volatility $\eta$. The quasi-closed formulas for $P_{V_{0}}$, $P^{\varepsilon}_{V_{1,0}}$ and $P^{\delta}_{V_{0,1}}$ are given in Appendix B.

6.3 Accuracy of the Approximation

We now state the precise accuracy result for the formal approximation determined in the previous sections. All the reasoning in Sections 6.1 and 6.2 were only formal arguments and well-thought choices for the proposed first-order approximations. The following theorem is the result that establishes the order of accuracy of this approximation and justifies, a posteriori, the choices made earlier. The proof is very similar to the ones presented in Fouque et al. (2011) and Fouque and Lorig (2011) and, therefore, omitted.

Theorem 6.2**.**

Under Assumption 2.2 and if the payoff functions \raisebox{1.5pt}{\varphi}_{S} and \raisebox{1.5pt}{\varphi}_{V} are continuous and piecewise smooth, then

[TABLE]

7 Conclusion

In this paper, we have proposed a continuous diffusion model for the stock price that is able to capture both skews in the stock’s and volatility index’s options data and that allows for quasi-closed formulas for the first-order approximation for option prices on the spot and its volatility index. These features were not achieve by any other continuous diffusion model. We have exemplified our calibration procedure with real data on S&P 500 and VIX.

Further research could be conducted in order to develop the generalizations outlined in Section 5. For instance, to be able to achieve the fit of VIX’s term structure, one could consider the method outlined in this paper to derive the first-order approximation under a time-dependent generalization of the model proposed here.

Appendix A Fourier Method to Compute the First-Order Approximation for Options on $S$

The computations presented here are based on the ideas shown in Fouque and Lorig (2011).

Let us first change variables to better apply the Fourier method.

[TABLE]

Therefore, one concludes

[TABLE]

A.1 A Quasi-Closed Formula for $P_{S_{0}}$

Define

[TABLE]

Hence, if we denote by $\widehat{P}_{S_{0}}(\tau,\xi,v,z)$ the Fourier transform of $\widetilde{P}_{S_{0}}(\tau,x,v,z)$ with respect to $x$, $\widehat{P}_{S_{0}}$ satisfies the following PDE

[TABLE]

Therefore, by the arguments presented in Heston (1993), we can write

[TABLE]

where

[TABLE]

For call options, we must set $\xi_{i}>1$.

A.2 A Quasi-Closed Formula for $P^{\varepsilon}_{S_{1,0}}$

If we denote by $\widehat{P}_{S_{1,0}}^{\varepsilon}(\tau,\xi,v,z)$ the Fourier transform of $\widetilde{P}_{S_{1,0}}^{\varepsilon}(\tau,x,v,z)$ with respect to $x$, $\widehat{P}_{S_{1,0}}^{\varepsilon}$ satisfies the following PDE

[TABLE]

where

[TABLE]

Consider now the following ansatz:

[TABLE]

So,

[TABLE]

and hence we should choose $f^{\varepsilon}_{0}$ and $f^{\varepsilon}_{1}$ to solve

[TABLE]

Therefore, once we solve the ODE system above, we need to compute

[TABLE]

A.3 A Quasi-Closed Formula for $P^{\delta}_{S_{0,1}}$

Now, if we denote by $\widehat{P}_{S_{0,1}}^{\delta}(\tau,\xi,v,z)$ the Fourier transform of $\widetilde{P}_{S_{0,1}}^{\delta}(\tau,x,v,z)$ with respect to $x$, $\widehat{P}_{S_{0,1}}^{\delta}$ satisfies the following PDE

[TABLE]

where

[TABLE]

Consider now the following ansatz:

[TABLE]

and notice

[TABLE]

Therefore,

[TABLE]

and hence we should choose $g^{\delta}_{0}$, $g^{\delta}_{1}$ and $g^{\delta}_{2}$ to solve

[TABLE]

Therefore, once we solve the ODE system above, we need to compute

[TABLE]

Appendix B Fourier Method to Compute the First-Order Approximation for Options on VIX

Firstly, notice that

[TABLE]

where $d\overline{V}_{t}=\kappa(m-\overline{V}_{t})dt+\overline{\eta}(z)\sqrt{\overline{V}_{t}}dW_{t}^{V}$. Hence,

[TABLE]

Using the Green’s function technique, as in Sepp (2008), the solution of this PDE can be written as

[TABLE]

where $\tau=T-t$, $\nu=\nu_{r}+i\nu_{i}$, $\nu_{r}\leq 0$ and

[TABLE]

Furthermore, the Fourier transform of some typical payoff functions are

[TABLE]

where $\mbox{erf}(z)$ is the complex error function, see Sepp (2008). For the VIX future case (\raisebox{1.5pt}{\varphi}_{V}(v)=v), one needs to remove the discount factor $e^{-r(T-t)}$ in $P_{V_{0}}$.

By equation (B.1),

[TABLE]

We then consider the following educated guess:

[TABLE]

where $h^{\varepsilon,\delta}_{V_{i}}$ are functions of $\tau$, $\nu$ and $z$. By Equations (6.73) and (6.77), one may conclude that $h^{\varepsilon,\delta}_{V_{i}}$ must solve the ODE system

[TABLE]

Funding

J.-P. Fouque was supported by NSF grant DMS-1409434.