2.1. Model formulation

Inspired by recent research on jump activity,Footnote ³ we consider a stochastic volatility model with both CJs and IJs in a time-varying jump intensity framework. The following dynamics are proposed to model the S&P 500 index, $S_t$, under the risk-neutral probability measure $\mathbb {Q}$:

(1) \begin{equation} \left\{\begin{aligned} \dfrac{\mathrm{d}S_t}{S_t} & = (r-q-\bar{\zeta}\lambda^{C}_t-\zeta\lambda^{S}_t)\,\mathrm{d}t +\sqrt{V_t}\,\mathrm{d}W^{S}_t + (e^{\bar{J}_S}-1)\,\mathrm{d}N^{C}_t +(e^{J_S}-1)\,\mathrm{d}N^{S}_t,\\ \mathrm{d}V_t & = \kappa(\theta-V_t)\,\mathrm{d}t + \sigma_V\sqrt{V_t}\,\mathrm{d}W^{V}_t + \bar{J}_V\,\mathrm{d}N^{C}_t +J_V\,\mathrm{d}N^{V}_t, \end{aligned} \right. \end{equation}

where $r$ is the risk-free rate and $q$ is the dividend. $W^{S}_t$ and $W^{V}_t$ are correlated Brownian motions with $\mathrm {d}W^{S}_t\,\mathrm {d}W^{V}_t=\rho \,\mathrm {d}t$. $\{N^{C}_t,N^{S}_t,N^{V}_t\}$ is an independent trivariate Poisson processes with time-varying rate $\lambda ^{C}_t$, $\lambda ^{S}_t$ and $\lambda ^{V}_t$, respectively. The jumps $\bar {J}_S$ and $\bar {J}_V$ occur simultaneously in the asset price and volatility and are driven by $N^{C}_t$. The independent price jump $J_S$ and volatility jump $J_V$ are driven by $N^{S}_t$ and $N^{V}_t$, respectively. The jump sizes are assumed to have the following distributions: $\bar {J}_V\sim \exp (\bar {\mu }_V)$ with mean $\bar {\mu }_V$, $\bar {J}_S\,|\,\bar {J}_V \sim N(\bar {\mu }_S+\rho _J\bar {J}_V,\bar {\sigma }_S^{2})$. In addition, $\bar {\zeta } \triangleq \mathbb {E}^{\mathbb {Q}}(e^{\bar {J}_S}-1) = e^{\bar {\mu }_S+{\bar {\sigma }_S^{2}}/{2}}/(1-\rho _J\bar {\mu }_V)-1$ and $\zeta \triangleq \mathbb {E}^{\mathbb {Q}}(e^{J_S}-1)=e^{\mu _S+{\sigma _S^{2}}/{2}}-1$.

Since the time-varying jump intensity is very popular in recent studies, for example, Bates [Reference Bates7] and Aït-Sahalia et al. [Reference Aït-Sahalia, Cacho-Diaz and Laeven1] report that more jumps occur during more volatile periods, we consider three cases as follows. (1) We first consider the stochastic volatility model with both CJs and IJs in Bandi and Reno [Reference Bandi and Reno3] and Chen and Ye [Reference Chen and Ye9], that is, $\lambda ^{i}_t=\lambda ^{i}$ is constant for $i\in \{C,S,V\}$. (2) The jump intensity is a linear function (e.g., [Reference Bardgett, Gourier and Leippold4,Reference Cheng, Ibraimi, Leippold and Zhang10]) of the instantaneous variance, and in this case, we name the model “SVCIJ-I” (note that the SVCIJ model is a special case of the SVCIJ-I model, where the dependent coefficient $\lambda ^{2,C}=\lambda ^{2,S}=\lambda ^{2,V}=0$). (3) The jump intensity ramps up in response to the occurrences of jumps, making further jumps more likely to follow (called “jump propagation”, see, e.g., [Reference Aït-Sahalia, Cacho-Diaz and Laeven1]), and in this case, we denote the model “SVCIJ-H”. A summary of the model specifications is given in Table 1. Our models embrace most of the important derivative pricing models in the literature (such as SVJ, SVIJ, SVCJ and SVCJ-H, i.e., [Reference Aït-Sahalia, Cacho-Diaz and Laeven1,Reference Bates5,Reference Eraker14,Reference Pan24]).

TABLE 1.

Summary of jump intensities.

Note. $J_{\lambda }^{i}\sim {\rm Exp}(\mu _{\lambda }^{i})$ with mean $\mu _{\lambda }^{i}$, $i\in \{C,S,V\}$.

The VIX is primarily related to spot volatility, and we are particularly interested in the effect of the flexible jump structure on spot volatility. In the next section, we study the distributional properties of volatility with respect to CJs and IJs.

2.2. Volatility distribution and jumps

To study how jumps affect volatility, we illustrate the probability density of volatility with CJs and IJs. The probability density function of volatility $\sqrt {V_T}$ can be obtained through Monte Carlo simulation with 500,000 sample paths. Then, we plot the probability density function of volatility at time $T$.

Figure 1 shows how the CJs and IJs of volatility pushes the distribution toward the equilibrium mean level. When the volatility process is examined separately, the CJs should yield the same effect strength as the independent jump. When the jump intensity increases, the volatility of volatility seems to be invariant and the shape of the distribution of volatility tends to be stable.

FIGURE 1.

Probability density function of volatility over a 1-year horizon with different jump intensities, based on the SVCIJ model. The parameter values are shown in Table 2; we simulate 500,000 sample paths of data with length of time $T$.

TABLE 2.

Summary of parameter settings.

Note. The superscript $i\in \{C,S,V\}$.

Naturally, the impact of contemporaneous and independent jumps on volatility cannot be ignored. Jumps in volatility can clearly improve understandings of the dynamics of spot volatility. In addition, Psychoyios et al. [Reference Psychoyios, Dotsis and Markellos27] argues that incorrectly omitting jumps may cause considerable problems in VIX futures and option pricing. Cheng et al. [Reference Cheng, Ibraimi, Leippold and Zhang10] argues that by accounting for return jumps and volatility jumps, the stochastic volatility model can further improve the performance of VIX option pricing formulas. In the next section, we derive formulas for VIX futures and option pricing.

3.1. Characteristic function

In this subsection, we present the conditional characteristic functions of $V_T$ and $(\lambda _T^{C},\lambda _T^{S},\lambda ^{V}_T)$ under the risk-neutral measure. The characteristic function $\mathcal {L}:\mathbb {C}^{1}\rightarrow \mathbb {C}^{1}$ of some process $X_T$ is defined as

$$\mathcal{L}(\phi;T-t,X_t)\triangleq \mathbb{E}_t(e^{\phi X_T}),$$

where $T\geqslant t$ and $\phi \in \mathcal {D}_x$, and $\mathcal {D}_x$ denotes the complex domain of $\mathcal {L}(\phi ;T-t,X_t)$, that is

$$\mathcal{D}_x=\{\phi \in \mathbb{C}: \mathbb{E}[\mathbb{E}_t(e^{\phi X_T})]<\infty\}.$$

The advantage of affine models is that their characteristic function has an exponentially affine form (see [Reference Duffie, Pan and Singleton13]), and thus, various transformations can be quasi-analytically calculated by solving a system of ordinary differential equations (ODEs).

Proposition 1 (Characteristic function of $V_T$ under the SVCIJ-I model)

If the spot variance $V_t$ follows the dynamics given by Eq. (1) and the intensities are represented as linear functions of $V_t$, Duffie et al. [Reference Duffie, Pan and Singleton13] shows that for any $\phi \in \mathcal {D}_v$, where $\mathcal {D}_v$ denotes the complex domain of $\mathcal {L}(\phi ;\tau,V_t)$, the exponentially affine expression below holds:

(6) \begin{equation} \mathcal{L}(\phi;\tau,V_t)\triangleq \mathbb{E}_t(e^{\phi V_T}) = e^{h_1(\phi,\tau)V_t + h_2(\phi,\tau)}, \end{equation}

where $\tau =T-t$ and $h_1(\phi,\tau )$, $h_2(\phi,\tau )$, $h_3(\phi,\tau )$ are solutions of the system of ODEs:

$$\left\{\begin{aligned} \dfrac{\partial}{\partial\tau}h_1(\phi,\tau) & ={-}\kappa h_1(\phi,\tau)+\frac{1}{2}\sigma_V^{2} h_1^{2}(\phi,t)+\lambda^{2,C}\left(\frac{1}{1-\bar{\mu}_Vh_1(\phi,\tau)}-1\right)\\ & \quad +\lambda^{2,V}\left(\frac{1}{1-\mu_Vh_1(\phi,\tau)}-1\right),\\ \dfrac{\partial}{\partial\tau}h_2(\phi,\tau) & = \kappa\theta h_1(\phi,\tau)+\lambda^{1,C}\left(\frac{1}{1-\bar{\mu}_Vh_1(\phi,\tau)}-1\right) +\lambda^{1,V}\left(\frac{1}{1-\mu_Vh_1(\phi,\tau)}-1\right), \end{aligned} \right.$$

with the initial conditions $h_1(\phi,0)=\phi$, $h_2(\phi,0)=0$.

Notably, the SVCIJ model is a special case of the SVCIJ-I model. When $\lambda ^{2,C}=\lambda ^{2,S}=\lambda ^{2,V}=0$, the SVCIJ-I model degenerates into the SVCIJ model. Although it is only the fading of the $\lambda ^{2,C}$, $\lambda ^{2,S}$ and $\lambda ^{2,V}$ parameters, the analytical solution exists for the characteristic function under the SVCIJ model.

Corollary 1 (Characteristic function of $V_T$ under the SVCIJ model)

The conditional characteristic function of the instantaneous variance $V_T$ is given by the following exponentially affine expression:

(7) \begin{equation} \mathcal{L}(\phi;\tau,V_t)\triangleq \mathbb{E}_t(e^{\phi V_T}) = e^{h_1(\phi,\tau)V_t + h_2(\phi,\tau)+h_3(\phi,\tau)}, \end{equation}

where $\tau =T-t$ and $h_1(\phi,\tau )$, $h_2(\phi,\tau )$, $h_3(\phi,\tau )$ are defined as follows:

$$\left\{ \begin{aligned} h_1(\phi,\tau) & = \frac{2\kappa \phi}{\sigma_V^{2}\phi+(2\kappa-\sigma_V^{2}\phi)e^{\kappa \tau}},\\ h_2(\phi,\tau) & =\frac{-2\kappa \theta}{\sigma_V^{2}}\ln{\left(1+\frac{\sigma_V^{2}\phi}{2\kappa}(e^{-\kappa\tau}-1)\right)},\\ h_3(\phi,\tau) & =\frac{2\lambda_{C}\bar{\mu}_{V}}{2\kappa\bar{\mu}_{V}-\sigma_V^{2}}\ln \left(1+\frac{(\sigma_V^{2}-2 \kappa \bar{\mu}_{V}) \phi}{2 \kappa(1-\bar{\mu}_{V} \phi)}(e^{-\kappa \tau}-1)\right)\\ & \quad +\frac{2 \lambda_{V} \mu_{V}}{2 \kappa \mu_{V}-\sigma_V^{2}} \ln \left(1+\frac{(\sigma_V^{2}-2 \kappa \mu_{V}) \phi}{2 \kappa(1-\mu_{V} \phi)}(e^{-\kappa \tau}-1)\right), \end{aligned} \right.$$

with the initial conditions $h_1(\phi,0)=\phi$, $h_2(\phi,0)=h_3(\phi,0)=0$, with $\phi \in \mathcal {D}_v$. In addition, under this model, $\phi _R$ satisfies

$$\phi_R<\min \left(\frac{2\kappa}{\sigma_V^{2}(1-e^{-\kappa\tau})},\frac 1\mu_V,\frac{2\kappa}{\sigma_V^{2}(1-e^{-\kappa\tau})+2\mu_V\kappa e^{-\kappa\tau}},\frac {1}{\bar{\mu}_V},\frac{2\kappa}{\sigma_V^{2}(1-e^{-\kappa\tau})+2\bar{\mu}_V\kappa e^{-\kappa\tau}}\right).$$

Proposition 2 (Characteristic function of $(V_T, \lambda ^{C}_T, \lambda ^{S}_T, \lambda ^{V}_T)$ under the SVCIJ-H model)

If the spot variance $V_t$ follows the dynamics given by Eq. (1) and the intensity satisfies the bottom row of Table 1, then the conditional joint characteristic function of the instantaneous variance and intensities, for any $\boldsymbol {\phi }$ that belongs to their complex domain of $\mathcal {L}(\boldsymbol {\phi };\tau,V_t,\lambda ^{C}_t,\lambda ^{S}_t,\lambda ^{V}_t)$, is given by the following exponential affine form:

\begin{align*} \mathcal{L}(\boldsymbol{\phi};\tau,V_t,\lambda^{C}_t,\lambda^{S}_t,\lambda^{V}_t)& \triangleq \mathbb{E}_t(e^{\phi_1 V_T+\phi_2\lambda_T^{C}+\phi_3\lambda_T^{S}+\phi_4\lambda_T^{V}}) \\ & =e^{h_1(\boldsymbol{\phi};\tau)V_t +h_2(\boldsymbol{\phi};\tau)\lambda^{C}_t+h_3(\boldsymbol{\phi};\tau)\lambda^{S}_t +h_4(\boldsymbol{\phi};\tau)\lambda^{V}_t+h_5(\boldsymbol{\phi};\tau)}, \end{align*}

where $\boldsymbol {\phi }=(\phi _1,\phi _2,\phi _3,\phi _4,\phi _5)$ and $h_1(\boldsymbol {\phi },\tau )$, $h_2(\boldsymbol {\phi },\tau )$, $h_3(\boldsymbol {\phi },\tau )$, $h_4(\boldsymbol {\phi },\tau )$, $h_5(\boldsymbol {\phi },\tau )$ are solutions of the system of ODEs:

$$\left\{ \begin{aligned} h_1(\boldsymbol{\phi},\tau) & = \frac{2\kappa \phi_1}{\sigma_V^{2}\phi_1+(2\kappa-\sigma_V^{2}\phi_1)e^{\kappa \tau}},\\ \frac{\partial}{\partial\tau}h_2(\boldsymbol{\phi},\tau) & ={-}\alpha^{C} h_2(\boldsymbol{\phi},\tau) +\frac{1}{1-\bar{\mu}_Vh_1(\boldsymbol{\phi},\tau)}\times\frac{1}{1-\mu_{\lambda}^{C}h_2(\boldsymbol{\phi},\tau)}-1,\\ \frac{\partial}{\partial\tau}h_3(\boldsymbol{\phi},\tau) & ={-}\alpha^{S} h_3(\boldsymbol{\phi},\tau) +\frac{1}{1-\mu_{\lambda}^{S}h_3(\boldsymbol{\phi},\tau)}-1,\\ \frac{\partial}{\partial\tau}h_4(\boldsymbol{\phi},\tau) & ={-}\alpha^{V} h_4(\boldsymbol{\phi},\tau) +\frac{1}{1-\mu_Vh_1(\boldsymbol{\phi},\tau)}\times\frac{1}{1-\mu_{\lambda}^{V}h_4(\boldsymbol{\phi},\tau)}-1,\\ \frac{\partial}{\partial\tau}h_5(\boldsymbol{\phi},\tau) & = \kappa\theta h_1(\boldsymbol{\phi},\tau)+\alpha^{C}\lambda^{C}_{\infty} h_2(\boldsymbol{\phi},\tau) + \alpha^{S}\lambda^{S}_{\infty} h_3(\boldsymbol{\phi},\tau)+ \alpha^{V}\lambda^{V}_{\infty} h_4(\boldsymbol{\phi},\tau), \end{aligned}\right.$$

with the initial conditions $h_1(\boldsymbol {\phi },0)=\phi _1$, $h_2(\boldsymbol {\phi },0)=\phi _2$, $h_3(\boldsymbol {\phi },0)=\phi _3$, $h_4(\boldsymbol {\phi },0)=\phi _4$, $h_5(\boldsymbol {\phi },0)=0$, with $\phi _1\sim \phi _4 \in \mathcal {D}_v$.

Since the above ODEs of $h_2$, $h_3$ and $h_4$ of the SVCIJ-I and SVCIJ-H modelsFootnote ⁴ have no analytical solutions, we use a numerical method, such as the Runge-Kutta algorithm, to solve these ODEs.

3.2. Pricing VIX futures and options

In this section, we present our model framework for pricing VIX futures and options, which is a general case of the previous affine pricing frameworks developed in the previous literature. First, simple affine state dynamics are employed to price the VIX in explicit form in general equilibrium; then, the generalized Laplace payoff transform analysis is conducted and performed to derive a pricing formula for VIX futures and options as a single integral.

The characteristic functions of variance and intensities have been discussed above. To find the characteristic functions of the VIX, a natural idea is to explore VIX squared as a linear function of instantaneous variance and intensities. Fortunately, Duan and Yeh [Reference Duan and Yeh12] provide a linear expression for the relationship between the VIX and variance in a stochastic volatility model. Therefore, we generalize and derive this expression under the SVCIJ-I and SVCIJ-H models. In Lemma 1, we derive the linear relation between the state variable (i.e., the instantaneous variance and time-varying intensity) and the 30-calendar-day VIX under these three models.

Lemma 1 (Linear property of VIX)

The $\text {VIX}^{2}$ is simply given by an affine transformation of the spot variance and time-varying intensities. The 30-day VIX index is a measure of the risk-neutral expected volatility of the SPX over the next 30 calendar days. Specifically, the calculation formula of VIX squared at a given time $t$ is given by

$$\text{VIX}_t^{2}(\tilde{\tau})=\frac{2}{\tilde{\tau}}\mathbb{E}_t \left[\int_t^{t+\tilde{\tau}}\left(\frac{\mathrm{d}S_u}{S_u}-\mathrm{d}\log S_u\right)\right],$$

where $\tilde {\tau } = \frac {30}{365}$. To express the VIX as an explicit function in state variables, we follow Duan and Yeh [Reference Duan and Yeh12] in using the affine property of VIX squared. Lemma 1 shows that we can transform VIX squared into state variables $V_t$ and $\lambda _t$.

1. (1) For the SVCIJ-I model,

   (8) \begin{equation} \text{VIX}_t^{2} = aV_t+b, \end{equation}
   where
   (9) \begin{equation} \left\{ \begin{aligned} a & = \eta(1+2\lambda^{2,C}(\bar{\zeta} - \bar{\mu}_S-\rho_J \bar{\mu}_V)+2\lambda^{2,S}(\zeta - \mu_S)),\\ b & = \frac{B}{A}(1-\eta)+2(\bar{\zeta} - \bar{\mu}_S-\rho_J \bar{\mu}_V)\left(\lambda^{1,C}+\lambda^{2,C}\frac{B}{A}(1-\eta)\right)\\ & \quad + 2(\zeta - \mu_S )\left(\lambda^{1,S}+\lambda^{2,S}\frac{B}{A}(1-\eta)\right). \end{aligned} \right. \end{equation}

   Here, $\bar {\mu }_V\triangleq \mathbb {E}(\bar {J}^{V})$, $\mu _V\triangleq \mathbb {E}(J^{V})$, $\bar {\zeta } \triangleq \mathbb {E}(e^{\bar {J}^{S}}-1)$, $\zeta \triangleq \mathbb {E}(e^{J^{S}}-1)$ and $\eta ={(1-e^{-A\tilde {\tau }})}/{A\tilde {\tau }}$, $A=\kappa -\lambda ^{2,V}\mu _V-\lambda ^{2,C}\bar {\mu }_V$ and $B=\theta \kappa +\lambda ^{1,V}\mu _V+\lambda ^{1,C}\bar {\mu }_V$. For convenience, we abbreviate $\text {VIX}_t(\tilde {\tau })$ as $\text {VIX}_t$. When $\lambda ^{2,C}=\lambda ^{2,S}=\lambda ^{2,V}=0$, this linear relationship can be applied to the SVCIJ model. This linear function can also be expressed as a linear function of $V_t$ and $\lambda _t$ if we substitute $\lambda _t=\lambda ^{1}+\lambda ^{2}V_t$ into the linear result above.

2. (2) For the SVCIJ-H model,

   (10) \begin{equation} \text{VIX}_t^{2} = aV_t+b\lambda_t^{C}+c\lambda_t^{S}+d\lambda_t^{V}+e, \end{equation}
   where
   (11) \begin{equation} \left\{ \begin{aligned} a & = \frac{1-e^{-\kappa\tilde{\tau}}}{\kappa\tilde{\tau}},\\ b & = \frac{\bar{\mu}_V(\varphi(\beta^{C})-a)}{\kappa-\beta^{C}}+2\varphi(\beta^{C})(\bar{\zeta} - \bar{\mu}_S-\rho_J \bar{\mu}_V),\\ c & = 2\varphi(\beta^{S})(\zeta - \mu_S ),\\ d & = \frac{\mu_V(\varphi(\beta^{V})-a)}{\kappa-\beta^{V}},\\ e & = \left(\theta+\frac{\tilde{\lambda}^{C}\bar{\mu}_V + \tilde{\lambda}^{V} \mu_V}{\kappa}\right)(1-a)+2 \tilde{\lambda}^{C} (\bar{\zeta} - \bar{\mu}_S-\rho_J \bar{\mu}_V)(1-\varphi(\beta^{C})) \\ & \quad + 2\tilde{\lambda}^{S} (\zeta - \mu_S )(1-\varphi(\beta^{S}))- \frac{\bar{\mu}_V\tilde{\lambda}^{C}}{\kappa-\beta^{C}}(\varphi(\beta^{C})-a)\\ & \quad-\frac{\mu_V\tilde{\lambda}^{V}}{\kappa-\beta^{V}}(\varphi(\beta^{V})-a), \end{aligned} \right. \end{equation}
   where $a\sim e$ are constants. Here, $\beta ^{i}=\alpha ^{i}-\mu _{\lambda }^{i}$, $\tilde {\lambda }^{i}=\lambda ^{i}_{\infty }\alpha ^{i}/\beta ^{i}$ and $\varphi (x)={(1-e^{-x\tilde {\tau }})}/{x\tilde {\tau }}$, $i\in \{C,S,V\}$.

Many studies make use of this linear relationship. Indeed, this linear relationship plays an important role in connecting the latent volatility to the observable VIX in the market. Lemma 1 derives the linear relationship between the $\text {VIX}_T$ and the two state variables: $V_T$ and $\lambda _T$. Thus, we can exploit this linear function relationship to obtain the explicit analytical formula for the prices of VIX futures and options under different models. As stated at the beginning of this section, we derive the formulas by the direct integral transformations to calculate the prices of VIX derivatives. See Theorem 1 for the calculation of $\text {VIX}$ futuresFootnote ⁵ and Theorem 2 for the calculation of $\text {VIX}$ options.

Theorem 1 (VIX future pricing)

Using the transformation of Eq. (4) and the characteristic function of VIX squared, the price of a VIX future, $F(T)$, is given by the following formula:

(12) \begin{equation} F(T)=\mathbb{E}[\text{VIX}_T]=\frac{1}{2\sqrt{\pi}}\int_0^{\infty} \frac{1-\mathbb{E}(e^{{-}s\text{VIX}_T^{2}})}{s^{3/2}}\,\mathrm{d} s. \end{equation}

1. (1) For the SVCIJ-I model,

   $$\mathbb{E}(e^{{-}s \text{VIX}_T^{2}})=e^{{-}bs}\mathcal{L}({-}as;T,V_t),$$
   where $s\in \mathbb {R}^{+}$, and the initial value $\text {VIX}^{2}_0=aV_0+b$ and $a, b$ are defined in Eq. (9).

2. (2) For the SVCIJ-H model,

   $$\mathbb{E}(e^{{-}s \text{VIX}_T^{2}})=e^{{-}es}\mathcal{L}(({-}as,-bs,-cs,-ds);T,V_t,\lambda_t^{C},\lambda^{S}_t,\lambda^{V}_t),$$
   and the initial value $\text {VIX}_0^{2} = aV_t+b\lambda ^{C}_0+c\lambda ^{S}_0+d\lambda _0^{V}+e$ and $a\sim e$ are defined in Eq. (11).

Notably, the integrals in this paper are integrable, the conditions of the Fubini theorem are satisfied, and the order of integration in Theorems 1 and 2 can be interchanged. Others, such as the Feller condition $2\kappa \theta >\sigma _V^{2}$, will not be elaborated in the rest of this paper, as the authors have verified its rationality.

Remark 1 The VIX futures prices have a single one-dimensional numerical integration, but for a hypothetical futures contract written on VIX squared (introduced by [Reference Eraker and Wu15]), which has a closed-form expression, since VIX squared is a linear function of $V_T$ and $\lambda _T$, the VIX squared futures prices under the SVCIJ-H model are

$$F_{\text{VIX}^{2}}(T)=\mathbb{E}(\text{VIX}_t^{2}) = a\mathbb{E}(V_t)+b\mathbb{E}(\lambda_t^{C})+c\mathbb{E}(\lambda_t^{S})+d\mathbb{E}(\lambda_t^{V})+e.$$

A European VIX call renders its holder the right, but not obligation, to obtain the difference between the VIX index at an expiration date $T$ and a pre-specified strike $K$. Compared to an equity option, it is more convenient to directly calculate price of a VIX call without normalization. The no-arbitrage condition implies that the VIX call price is represented by the following theorem.

Theorem 2 (VIX call option pricing)

Using the transformation of Eq. (5) and the characteristic function of VIX squared, the price of a VIX call option, $C(T,K)$, is given by the following formula:

(13) \begin{equation} C(T,K) =\mathbb{E}[(\text{VIX}_T-K)^{+}]=\frac{e^{{-}rT}}{2\sqrt{\pi}}\int_0^{\infty} \operatorname{Re}\left[\frac{1-\text{erf}(K\sqrt{\phi})}{\phi^{3/2}}\mathbb{E}(e^{\phi \text{VIX}_T^{2}})\right]\mathrm{d}\phi_I, \end{equation}

where $\phi =\phi _R + i \phi _I$ is the complex Laplace transform variable with $\phi _I\in \mathbb {R}$ and $\phi _R>0$ and where the complex error function $\text {erf}(Z)$ is defined by

\begin{equation}\text{erf}(Z) = \frac{2}{\sqrt{\pi}}\int_0^{Z} e^{{-}s^{2}}\,\mathrm{d}s.\end{equation}

1. (1) For the SVCIJ-I model,

   $$\mathbb{E}(e^{\phi \text{VIX}_T^{2}})=e^{b\phi}\mathcal{L}(a\phi;T,V_t),$$
   with the initial value $\text {VIX}^{2}_0=aV_0+b$.

2. (2) For the SVCIJ-H model,

   $$\mathbb{E}(e^{\phi \text{VIX}_T^{2}})=e^{e\phi}\mathcal{L}((a\phi,b\phi,c\phi,d\phi);T,V_t,\lambda_t^{C},\lambda^{S}_t,\lambda^{V}_t),$$
   with the initial value $\text {VIX}_0^{2} = aV_t+b\lambda ^{C}_0+c\lambda ^{S}_0+d\lambda _0^{V}+e$.

Our formulas are extremely efficient and accurate for calculating the futures and options prices of $\text {VIX}_T$, regardless of whether the model considers both CJs and IJs or time-varying intensities. We use clever integral transformations, and the formula involves only a single integral and thus exhibits more efficiency. In the next section, we perform some numerical experiments to verify the correctness of our formulas by comparing the results calculated by the Monte Carlo simulation with the results obtained by our formulas.

4.1. Monte Carlo simulations

Euler discretization can be used to approximate the path of the variance and intensity processes over a discrete interval; in this case, we choose a discrete time small enough for $\Delta t=1/1,000$. The discretization of the variance and intensity processes are

\begin{align*} V_t& =V_{t-1}+\kappa(\theta-V_{t-1})\Delta t+\sigma_V\sqrt{|V_{t-1}|}\sqrt{\Delta t}\epsilon_t+\bar{J}_VN_t^{C}+J_VN^{V}_t,\\ \lambda_t^{i}& =\lambda^{1,i}+\lambda^{2,i}V_t,\quad (\text{SVCIJ-I model}, i\in\{C,V\}),\\ \lambda_t^{i}& =\lambda_{t-1}^{i}+\alpha^{i}(\lambda_{\infty}^{i}-\lambda_{t-1}^{i})\Delta t+J_{\lambda}^{i}N^{i}_t,\quad (\text{SVCIJ-H model},i\in\{C,S,V\}), \end{align*}

where $\epsilon _t\sim \mathcal {N}(0,1)$, $\bar {J}_V\sim \exp {(\bar {\mu }_V)}$, $J_V\sim \exp {(\mu _V)}$, $J^{i}_{\lambda }\sim \exp {(\mu _{\lambda }^{i})}$, $N_t^{i}\sim {\rm Bernoulli}(\lambda _{t-1}^{i}\Delta t)$ ($\lambda _{t-1}^{i}$ is constant for the SVCIJ model).

With the sample path of variance and intensities (if any), we calculate the sampling VIX path according to $\text {VIX}_T^{2}=aV_T+b\lambda ^{C}_T+c\lambda ^{S}_T+d\lambda ^{V}_T+e$ or $\text {VIX}_T^{2}=aV_T+b$. Then, we obtain the VIX futuresFootnote ⁶ by $({1}/{N})\sum _{i=1}^{N}\text {VIX}_{T,i}$ and VIX options by $({1}/{N})\sum _{i=1}^{N}(\text {VIX}_{T,i}-K)^{+}$, where $N$ is the number of trajectories and $\text {VIX}_{T,i}$ is the value in $i$th trajectory at time $T$.

To measure the error between the Monte Carlo method and the results calculated by our formulas, we provide two measures: absolute error and relative error, which are defined as follows:

\begin{align*} \text{absolute error}& =|V-V_{\text{MC}}|,\\ \text{relative error}& =\frac{|V-V_{\text{MC}}|}{V_{\text{MC}}}, \end{align*}

where $V=F$, $C$ represent futures and options, respectively.

In our experiments, we supposed all parameters (see Table 2). In fact, suitable arbitrary parameters have no effect on the accuracy of our formula calculation.

4.2. Accuracy check

Table 3 shows the comparison of VIX futures prices obtained from our formula and Monte Carlo simulations, Table 4 shows the comparison of VIX options prices obtained from our formula and Monte Carlo simulations. The Monte Carlo simulation uses 200,000 sample paths and each path with time steps of 1/1,000. We list the results of the three models for given parameters and find that all three models match perfectly with those obtained from the Monte Carlo simulations.

It is clear that the differences between the results of our futures pricing formula and the Monte Carlo simulation results (the 200,000 Monte Carlo simulated trajectories) are very small, and the relative error range is between 0% and 0.13%. This confirms the accuracy of our futures calculation formula.

The differences between our formula's option prices and the simulation results are small, in the range of 0% to 1.75%. If we regard the Monte Carlo price as the benchmark, then this adequately demonstrates that both our option pricing formula and futures pricing formula are accurate. We have reason, theoretically, to believe that even with different calculation formulas, as long as the given parameters and model are the same, the pricing results are consistent. For the SVCIJ-H model, since the SVCIJ model is a special case of the SVCIJ-I model, to more reliably verify the formula of the SVCIJ-I model, we take $\lambda ^{2,C}$, $\lambda ^{2,S}$ and $\lambda ^{2,V}$ to be the large value of 20.

Compared with futures contracts written on commodities and equities, VIX futures are quite unusual. As time to maturity increases substantially, the prices of VIX futures flatten out and become less sensitive to the time to maturity. In the limiting case, VIX futures prices tend toward a constant, this phenomenon of the term structure of VIX futures prices is in line with the observed traded VIX futures prices in the CBOE (see, e.g., [Reference Eraker and Yang16]). This phenomenon is due to volatility's mean-reverting nature, which means that over the long run, volatility tends to revert its long-run average. For the SVCIJ-I model, $\mathbb {E}(V_T)\stackrel {T\rightarrow \infty }{\longrightarrow }(\theta \kappa +\lambda ^{1,C} \bar {\mu }_V+\lambda ^{1,V}\mu _V)/(\kappa -\lambda ^{2,C}\bar {\mu }_V-\lambda ^{2,V}\mu _V)$ is the long-run mean of variance. For the SVCIJ-H model, $\mathbb {E}(V_T)\stackrel {T\rightarrow \infty }{\longrightarrow } (\kappa \theta +\tilde {\lambda }^{C}\bar {\mu }_V+\tilde {\lambda }^{V}\mu _V)/\kappa$ is the long-run mean of variance, and $\lim _{T\rightarrow \infty }\mathbb {E}(\lambda ^{i}_T)\stackrel {T\rightarrow \infty }{\longrightarrow }\tilde {\lambda }^{i}$ ($i\in \{C,S,V\}$) are the long-run means of intensities. All of these limits above mean that as time-to-maturity increases, prices are becoming less sensitive to initial volatility or jump intensities changes. The existence of these limits confirms that as the time-to-maturity increases, the futures price tends toward a constant.

(14) \begin{equation} \lim_{T\rightarrow\infty} F(T) =\text{Constant}, \end{equation}

where the constant is a value independent of the VIX. The left panel of Figure 2 illustrates the prices of VIX futures approaching a constant (based on the SVCIJ-H model).

FIGURE 2.

The prices of VIX futures and options over long maturities, based on the SVCIJ-H model.

For this property of volatility, it is natural that the price of its options tends to be close to zero at long maturities. We can also derive this from our formula for options pricing (Eq. (13)), where the integral is bounded, and $e^{-rT}$ goes to zero as $T$ increases.

(15) \begin{equation} \lim_{T\rightarrow\infty} C(T,K) =0. \end{equation}

The right panel of Figure 2 illustrates the prices of VIX options are approaching a constant. Therefore, at long maturities, the VIX call option will become worthless and less sensitive (which is in line with the results of the models of [Reference Eraker and Yang16,Reference Lian and Zhu21]).

4.3. Sensitivity analysis

To further study the jump structure, we investigate the sensitivity of the VIX options with respect to the parameters of the jump parameters in Figures 3–5. To identify the ceteris paribus effects on the VIX options, we change one or more parameters while holding all remaining parameters fixed. For each parameter set, we plot the 3D diagram of VIX options against time to maturity and the parameters to be studied.

FIGURE 3.

The prices of VIX options ($K=25$), based on the SVCIJ model.

Naturally, considering the upward jump risk in the model, the call price would increase. Eraker et al. [Reference Eraker, Johannes and Polson17] shows that the SVCJ and SVIJ models exhibit similar behavior, but the SVIJ model is more flexible than the SVCJ model. Therefore, we may also be concerned about how the jump factors affect VIX futures and options in detail under the more general jump-diffusion models. This section examines how jump factors impact VIX derivatives pricing. These prices of VIX derivatives are calculated from our pricing formulas. Figure 3 illustrates the sensitivity of a VIX call option to CJs and IJs in volatility and time to maturity under the SVCIJ model.

Figure 3 illustrates the scenarios for different types of jumps, showing that the call price is affected by various jumps (by setting different jump intensities), especially for CJs. In particular, the top panels of Figure 3 show that the VIX call price is more sensitive to CJs for long-term contracts. These results are in line with the markets because CJs are frequently accompanied by crises and systemic risks, while most IJs occur independently of variance and are relatively accidental. Therefore, from the perspective of the pricing of VIX derivatives, CJs deserve more attention.

Figure 4 illustrates the scenarios for different $\lambda ^{2,C}$, $\lambda ^{2,S}$ and $\lambda ^{2,V}$, showing that the call prices are relatively sensitive to the time-varying intensities, especially for large $\lambda ^{2,C}$, $\lambda ^{2,S}$ and $\lambda ^{2,V}$. Alternatively, by removing the time-varying component of intensities, which is achieved by setting $\lambda ^{2,C}=\lambda ^{2,S}=\lambda ^{2,V}=0$, we can obtain the SVCIJ model. In addition, a large $\lambda ^{2,C}$, $\lambda ^{2,S}$ and $\lambda ^{2,V}$ are also associated with higher derivative prices and is more sensitive to time to maturity. Therefore, when the coefficient $\lambda ^{2,C}$, $\lambda ^{2,S}$ and $\lambda ^{2,V}$of the time-varying component of jump intensity may be large, we may not be able to ignore its effect on pricing.

FIGURE 4.

The prices of VIX futures and options ($K=26$), based on the SVCIJ-I model.

Figure 5 shows the sensitivity of a VIX call price to the impacts of jumps in intensity and time to maturity in the short run and long run (by setting $\alpha ^{C}$, $\alpha ^{S}$ and $\alpha ^{V}$ to a large value or small value), which shows that the VIX call prices with $\alpha ^{C}=\alpha ^{S} =\alpha ^{V}=2$ (long run) are more sensitive than VIX call prices with $\alpha ^{C}=\alpha ^{S}=\alpha ^{V}=15$ (short run). This is because of the significant difference in the mean-reversion speed of intensity. For the SVCIJ-H model, VIX squared is not only related to spot variance but also linearly related to time-varying jump intensities. Thus, intensities can affect VIX squared not only by acting on spot variance but can also directly affect VIX squared.

FIGURE 5.

The prices of VIX options (K=26), based on the SVCIJ-H model.

In summary, although variance jumps in CJs and IJs are both important to the pricing of VIX options, CJs seem to be more important. The time-varying intensity as a linear function of variance seems to be important in pricing VIX derivatives when variance has a large effect on intensity. More important, the mean-reversion of intensity mitigates the value of jumps in intensity. In the case of fast mean-reversion speed, the jumps in intensity seem to have little effect on pricing VIX calls.

4.4. Impact of jumps in implied volatility

To characterize the pricing of VIX call options, we investigate implied volatility. Specifically, we perform numerical tests to investigate the behavior of the implied volatility of VIX options. Recall that the implied volatility $\sigma _t^{{\rm implied}}$ associated with $F(T)$ and $C(T,K)$ is defined by inverting the Black and Scholes [Reference Black and Scholes8] formula (see Appendix A for details). Figure 6 shows the implied volatility with the strike price under different maturities and $V_0$. Figures 7 and 8 illustrate the implied volatility patterns generated by varying jump parameters. Table 5 reports the implied volatility with respect to the maturity and strike under the SVCIJ model with different jump patterns.

FIGURE 6.

Implied volatility for VIX option in the SVCIJ model. Note. In the left panel, we set varying values of maturity in the SVCIJ model. In the right panel, we set varying values of spot variance $V_0$ and a fixed maturity at $\tau =0.5$ in the SVCIJ model.

FIGURE 7.

Implied volatility for VIX option with maturity $\tau =0.5$ with respect to the SVCIJ model and SVCIJ-I model. Note. In the left panel, we set varying values of the jump intensity of contemporaneous and volatility-independent jumps in the SVCIJ model. In the right panel, we set varying values of the dependent coefficient in the SVCIJ-I model.

FIGURE 8.

Implied volatility for VIX option with maturity $\tau =0.5$ in the SVCIJ-H model. Note. In the left panel, we set the mean reversion of both intensities very low ($\alpha ^{C}=\alpha ^{S}=\alpha ^{V}=2$, implying long-run intensities) and set varying values of jump size intensities in the SVCIJ-H model. In the right panel, we set the mean reversion of both intensities very high ($\alpha ^{C}=\alpha ^{S}=\alpha ^{V}=15$, implying short-run intensities) and set varying values of jump size intensities in the SVCIJ-H model.

TABLE 5.

Implied volatility for the SVCIJ model.

Note. The table reports implied volatility for VIX options with respect to maturity and strike in the SVCIJ model, as well as the jump intensities of contemporaneous and volatility-independent jumps. $\tau$ is the maturity, and $\lambda ^{C}$ and $\lambda ^{V}$ are jump intensities.

There are several noteworthy features of implied volatility.

First, in the SVCIJ model, for a given strike, the implied volatility is greater for short-maturity options (left panel of Figure 6). That is, the term structure of implied volatility is downward sloping irrespective of the level of VIX. To understand why this is, it is crucial to understand that implied volatility, in this case calculated from the Black–Scholes model for pricing an option on VIX futures, assumes that the underlying follows a random walk. If a given process is a random walk, its forecasted variance increases linearly with the forecast horizon. VIX options exhibit a downward-sloping term structure, and therefore, the VIX variance does not increase proportionally with the forecast horizon. Rather, it suggests that the VIX is not a martingale.

Second, the shapes of the implied volatility uniformly form a concave “frown” rather than the usual convex “smile” seen in most equity options, such as SPX options.

Third, low (high) strikes, over some strike range, generally have low (high) implied volatilities, suggesting a positively skewed underlying distribution. Unlike equity options, which exhibit a negative skew in implied volatility, implied volatility in VIX options exhibits positive skew. Notably, accounting for positive jumps in volatility can explain the positive implied volatility skew exhibited by options written on volatility such as CBOE's VIX options. For the SVCIJ-H model, VIX squared would be a linear function of $V_t$ and $\lambda _t$. Similarly, the positive jumps in $\lambda _t$ can also be explained by the positive VIX implied volatility skew.

Fourth, the implied volatility values decrease as the initial variance $V_0$ increases but not significantly.

In the rest of this section, we investigate the sensitivity of the implied volatility of VIX options to the jump-related parameters in the proposed models. To identify the ceteris paribus effects on the implied volatility, we change the specified jump-related parameters while holding all remaining parameters fixed. For each parameter set, we plot the implied volatility curve against the strike.

For the SVCIJ model, Figure 7 shows the implied volatility under different jump structures against the strike. As the figure shows, the CJs and IJs are evident: the structures under CJs and IJs have different implied curves, but for a given jump intensity, the implied volatility is greater for IJs. The implied volatility of IJs is indistinguishable, and the implied volatility curve is very close. However for CJs, there is a significant decreasing trend given increases in $\lambda ^{V}$ over some strike range. Moreover, the impact of these CJ and IJ jump structures on implied volatility is consistent at different maturity horizons. Table 5 shows the VIX implied volatility surface under different jump structures against strike and maturity. In addition, we know from the previous section on option sensitivity that IJs are not sensitive to option prices, while CJs are considerably more sensitive. This point of view reflects the significance of distinguishing CJs and IJs in our model. During the financial crisis, there were relatively more CJs, while IJs generally predominate (evidence from [Reference Chen and Ye9]). This is reflected by the fact that our stochastic volatility model contains both CJs and IJs.

For the SVCIJ-I model, the left panel of Figure 7 shows the implied volatility against the dependent coefficient $\lambda ^{2,C}$ and $\lambda ^{2,V}$. We observe that the yellow curve is close to the green curve and that the red curve is close to the purple curve, which means that if $\lambda ^{2,C}$ and $\lambda ^{2,V}$ are of the same value, then the impact on implied volatility is similar. However, note that, in the low strike, the implied volatility on $\lambda ^{2,C}$ is smaller, in the high strike, it is greater, relative to implied volatility on $\lambda ^{2,V}$ of the same value. In the low strike, the implied volatility in the non-zero dependent coefficient is close to the implied volatility for a zero dependent coefficient. In our SVCIJ-I model, therefore, when a strike is small, the effect of the dependent coefficient on implied volatility is relatively small.

In the previous section, we illustrated the option sensitivity in the long run ($\alpha ^{C}=\alpha ^{S}=\alpha ^{V}=2$) and short run ($\alpha ^{C}=\alpha ^{S}=\alpha ^{V}=15$). We are interested in the implied volatility patterns generated by different mean reversions of jump intensity. Figure 8 shows the result when we condition on high and low mean reversion intensity. We see that the implied volatility curves are almost always increasing with respect to jump size intensity increases, especially under long-run mean reversion. In the short-run case, jump size intensity has no significant effect on implied volatility. However, for a low strike, it does not appear to have much effect on implied volatility. Our model shows that as maturity increases, the futures price also increases and exhibits the concave “frown” pattern. This implies that VIX option-implied volatility is positively correlated with the VIX itself. In the SVCIJ-H model, a fluctuating $\lambda _t$ process is important in shaping the positive correlation between implied volatility and the VIX. This can be interpreted to mean that, first, as $\lambda _t$ increases, it drives up the VIX because the VIX loads positively on it; second, it also drives up implied volatility and thereby increases the prices of VIX calls. Part of this phenomenon is strikingly similar to the price sensitivity of options discussed previously.

4.5. Empirical studies

In this section, we introduce market data and investigate whether the flexible jump structure improves pricing performance of VIX options. The market data used in this paper include VIX futures, LIBOR and VIX options that provided by the CBOE and OptionMetrics. The sampling dates are June 21, 2017 (Wednesday) as in-sample data for model calibration, and June 22, 2017 (Thursday) as out-of-sample data for evaluating model performance. We filter the option data as follows: (1) we use the mid price (the average of the best bid and best ask) as the market VIX option price; (2) we eliminate any options with zero bid-quote and trading volumes; and (3) we consider only VIX call options because the trading volume of VIX calls is more than puts. Table 6 provides the summary statistics of VIX option data divided into three categories based on moneyness. Here, moneyness is defined as $m:=\ln (K/F)$ with the strike price of option $K$ and corresponding VIX future price $F$.

TABLE 6.

Summary statistics of VIX option data.

Note. The table is divided into three categories based on moneyness (defined as $m:=\ln (K/F)$ with the strike price of option $K$ and corresponding VIX future price $F$) on June 21, 2017 (in-sample) and June 22, 2017 (out-of-sample). The reported observation is the number of VIX call option contracts.

For ease of comparison among different models, four typical models are considered: SV model, SVCIJ model and two generalized forms of SVCIJ model by allowing independent return jump to be self-exciting (i.e., SVCIJ-I and SVCIJ-H). Table 7 presents the estimated results by minimizing the mean-square error between market prices and model option prices. The estimated value of $\lambda ^{2,S}$ is 139.2725 of SVCIJ-I model, which implies that state-depend jump intensity may not be ignored. The estimated value of $\alpha ^{S}$ and $\mu _{\lambda }^{S}$ of the SVCIJ-H model are 11.2192 and 11.6707, respectively, which implies that Hawkes-type propagation also plays a very important role.

TABLE 7.

Calibrated model parameters.

Note. The table exhibits the calibration results of VIX option for SV, SVCIJ, SVCIJ-I, and SVCIJ-H modes. The parameters are backed out by minimizing the mean-square error between market and model option prices by using the in-sample data (June 21, 2017).

Next, we compare the performance among different models based on the mean absolute percentage errors (MAPE), which is defined as:

$$\text{MAPE}=\frac{1}{N_t}\sum_{i=1}^{N_t}\frac{|C_i^{{\rm Model}}-C_i^{{\rm Market}}|}{C_i^{{\rm Market}}},$$

where $N_t$ is the number of VIX call option contracts on day $t$, $C_i^{{\rm Market}}$ is the market price of contract $i$ and $C_i^{{\rm Model}}$ represents the model-implied price of contract $i$.

Figure 9 presents a cross-sectional VIX option pricing performance on in-sample data for different maturities and moneyness. It is clear that our proposed models can match market VIX option prices well for the long maturity date and OTM options. Furthermore, the advantage of our models over the benchmark model are most significant for short maturity date and ITM options. Figure 10 shows a cross-sectional VIX option pricing performance on out-of-sample data for different maturities and moneyness and have a similar behavior with in-sample result. Table 8 gives the in-sample and out-of-sample pricing errors of VIX option across different categories of moneyness and maturity, respectively. Clearly, both SVCIJ-I and SVCIJ-H have better performance than other models in all categories. Additionally, we find that the SVCIJ-I and SVCIJ-H have similar performance. Note that the long-term maturity and OTM options show the lowest price errors. Overall, our models provide the best in-sample fitting and out-of-sample prediction across all moneyness and maturities.

FIGURE 9.

In-sample pricing performance of VIX options against moneyness on June 21, 2017. Note. The moneyness is defined as $\ln (K/F)$ with the strike price of option $K$ and corresponding VIX future price $F$.

FIGURE 10.

Out-of-sample pricing performance of VIX options against moneyness on June 21, 2017. Note. The moneyness is defined as $\ln (K/F)$ with the strike price of option $K$ and corresponding VIX future price $F$.

TABLE 8.

VIX option pricing errors.

Note. It presents in-sample (June 21, 2017) and out-of-sample (June 22, 2017) option pricing errors for different ranges of moneyness and maturity for the SV, SVCIJ, SVCIJ-I and SVCIJ-H models.

Abbreviations: MAPE, mean absolute percentage errors; m, moneyness; DTM, days to maturity.