JEL Classification: G12, G13
Introduction
Since the global financial crisis of 2008, the VIX has become a central gauge of financial-market uncertainty and the reference index for a major market in volatility derivatives (Whaley [28 ] , Luo et al. [23 ] ). The index itself is not directly tradable: it is calculated from S&P 500 option prices and summarizes the market’s risk-neutral expectation of S&P 500 variation over the next 30 days (Luo and Zhang [22 ] ). Tradable exposure is instead provided by VIX futures and options, which allow investors to hedge volatility risk and take directional positions on changes in market uncertainty (Luo et al. [23 ] ). At the same time, evidence that volatility exhibits rough, non-Markovian dynamics has prompted growing interest in the implications of roughness for VIX futures and options (Gatheral et al. [14 ] , Jacquier et al. [18 ] , Rømer [27 ] ). The prominence of the VIX market, together with this evidence, creates a need for models that can capture the dynamics of the VIX while remaining sufficiently tractable for derivative pricing.1
Existing models for VIX derivatives can be organized into two broad approaches. The first derives VIX from the joint risk-neutral dynamics of the S&P 500 and its latent variance process, thereby maintaining a structural link between SPX and VIX derivatives, but often requiring a relatively rich specification and joint calibration to both markets [30 , 31 , 9 , 20 , 4 , 13 ] . The second approach models variance, volatility, VIX, or VIX futures directly [21 , 22 , 26 , 17 , 23 ] (Dupoyet et al. 2011). By avoiding an explicit model for the SPX spot price, direct specifications are generally more parsimonious and yield more tractable VIX derivative prices, although they do not impose a structural link between the equity and volatility markets.
This paper follows the instantaneous squared VIX (ISVIX) approach developed by Luo and Zhang [22 ] , Luo et al. [23 ] , which models directly the short-horizon limiting object whose conditional average determines squared VIX. The resulting Markovian specification is parsimonious and tractable, but it restricts the transmission of current volatility conditions across maturities to the pattern implied by exponential mean reversion. A natural way to relax this restriction is to introduce memory. Rough-volatility models use singular fractional kernels to capture the low local regularity and history dependence observed in volatility dynamics [5 , 14 , 12 ] . Most applications to VIX derivatives, however, introduce roughness into the latent variance or forward-variance process of the equity index and then derive VIX as a nonlinear functional of that process. This route is well suited to joint SPX–VIX modeling, but it makes the VIX term structure and VIX option prices depend on a comparatively rich latent-state system.
We instead propose a rough instantaneous squared VIX (RISVIX) model, by introducing a fractional Volterra kernel directly into the ISVIX framework of Luo et al. [23 ] . We derive an explicit VIX term structure governed by a Mittag–Leffler decay function and a semi-closed pricing formula for VIX options based on an affine Fourier–Laplace transform. These results lead to a two-step procedure in which the VIX term structure is used to estimate the memory and mean-reversion parameters, while VIX options are used to calibrate the remaining distributional parameters. Empirically, we estimate RISVIX using VIX term-structure indices and VIX option prices from 2018 to 2025. The gains from fractional memory vary substantially with market conditions. RISVIX improves the fit of the VIX term structure, with larger gains for its slope and curvature and during periods of elevated volatility. It also improves out-of-sample VIX option pricing, reducing MSE and MAE by 13.06% and 6.33%, respectively, with larger gains at shorter maturities and in high-VIX states.
Our contributions are threefold. First, we introduce fractional memory into an ISVIX model, yielding a parsimonious rough framework for modeling the VIX term structure and pricing VIX options. Second, we retain the tractability of direct VIX modeling by deriving an explicit term-structure formula and a semi-closed-form option-pricing formula, which together support a transparent two-step estimation and calibration procedure. Third, we show empirically that the gains from fractional memory vary substantially across volatility regimes: RISVIX improves the fit of the VIX term structure and the out-of-sample pricing of VIX options more strongly when volatility is elevated.
Related literature. Our paper is related to two strands of literature. First, it contributes to the literature on direct modeling of VIX term structures and derivatives. Lin [21 ] price VIX futures by integrating information under the physical and risk-neutral measures, while Dupoyet et al. (2011) develop a direct volatility-futures model and compare CEV and CIR specifications with and without jumps. Huskaj and Nossman [17 ] model the term structure of VIX futures directly rather than deriving it from an equity-variance process, and Mencía and Sentana [26 ] conduct an extensive empirical evaluation of VIX derivative valuation models across different market conditions. More closely related to our work, Luo and Zhang [22 ] extend VIX to arbitrary horizons and introduce instantaneous and long-term squared VIX as the short- and long-maturity limits of the VIX term structure. Building on this representation, Luo et al. [23 ] develop a unified Markovian ISVIX model with a stochastic mean level and jumps for pricing VIX derivatives across horizons. RISVIX extends their framework by replacing the propagation pattern implied by Markovian mean reversion with a fractional one.
Second, our paper contributes to the literature on rough volatility and VIX derivatives. Rough-volatility models use fractional kernels to capture the low local regularity and history dependence of volatility [5 , 14 , 12 ] . Existing applications to VIX derivatives generally specify rough dynamics for the latent variance or forward-variance curve of the equity index and then derive VIX from that process [18 , 16 , 27 , 8 , 19 ] . The closest study is Ye et al. (2026), who derive the VIX term structure from rough Heston using a Markovian approximation. Our RISVIX differs from their approach in two respects. First, it retains the fractional Volterra kernel directly in the ISVIX dynamics, thereby preserving the non-Markovian memory structure without relying on a Markovian approximation. Second, its affine Volterra structure yields a semi-closed pricing formula for VIX options, allowing the VIX term structure and option prices to be studied within the same framework.
The remainder of the paper is organized as follows. § 2 reviews the ISVIX framework, fractional Brownian motion and affine Volterra processes. § 3 introduces RISVIX, derives the VIX term structure and option-pricing formula, and presents the estimation procedure. § 4 reports the empirical results for the VIX term structure and VIX options. § 5 concludes. The Appendix contains the proofs and auxiliary results.
Preliminaries
This section reviews instantaneous squared VIX, fractional Brownian motion and affine Volterra models, which will be used throughout this paper. We work on a filtered probability space ( Ω , F , ( F t ) t ≥ 0 , Q ) (\Omega,\mathcal{F},(\mathcal{F}_t)_{t\ge0},\mathbb{Q}) ( Ω , F , ( F t ) t ≥ 0 , Q ) satisfying the usual conditions, where Q \mathbb{Q} Q is a risk-neutral measure. Write E t Q [ ⋅ ] = E Q [ ⋅ ∣ F t ] . \mathbb{E}^{\mathbb{Q}}_t[\cdot] = \mathbb{E}^{\mathbb{Q}}[\cdot\mid\mathcal{F}_t]. E t Q [ ⋅ ] = E Q [ ⋅ ∣ F t ] .
Instantaneous Squared VIX
Let S t S_t S t be a strictly positive semimartingale representing the equity index. We use the continuous-strike idealisation of VIX, whose squared value over a horizon τ > 0 \tau> 0 τ > 0 is represented by the log contract ([21 , 22 ] ):
V I X t , τ 2 = 2 τ E t Q [ ∫ t t + τ d S s S s − − log ( S t + τ S t ) ] . (1) \mathrm{VIX}_{t,\tau}^{2} = \frac{2}{\tau}\mathbb{E}^{\mathbb{Q}}_t\left[\int_t^{t+\tau}\frac{\mathrm{d}S_s}{S_{s-}}-\log\left(\frac{S_{t+\tau}}{S_t}\right)\right]. \tag*{(1)} VIX t , τ 2 = τ 2 E t Q [ ∫ t t + τ S s − d S s − log ( S t S t + τ ) ] . (1)
Suppose that the continuous component of the index return has instantaneous variance v t v_t v t , while log-price jumps have intensity λ t x \lambda_t^x λ t x and sizes distributed as a random variable x x x with a state-independent law. When the integrated variance and jump correction have finite expectations, Equation (1) admits the instantaneous squared VIX (ISVIX) representation of Luo and Zhang [22 ] :
V t = v t + 2 λ t x E Q [ e x − 1 − x ] , V I X t , τ 2 = 1 τ E t Q [ ∫ t t + τ V s d s ] . (2) V_t = v_t + 2\lambda_t^x\mathbb{E}^{\mathbb{Q}}\left[e^x-1-x\right], \qquad\mathrm{VIX}_{t,\tau}^{2} = \frac{1}{\tau}\mathbb{E}^{\mathbb{Q}}_t\left[\int_t^{t+\tau}V_s\,\mathrm{d}s\right]. \tag*{(2)} V t = v t + 2 λ t x E Q [ e x − 1 − x ] , VIX t , τ 2 = τ 1 E t Q [ ∫ t t + τ V s d s ] . (2)
If V V V is right-continuous in L 1 ( Q ) L^1(\mathbb{Q}) L 1 ( Q ) , then V I X t , τ 2 → V t \mathrm{VIX}_{t,\tau}^{2}\to V_t VIX t , τ 2 → V t in L 1 ( Q ) L^1(\mathbb{Q}) L 1 ( Q ) as τ ↓ 0 \tau\downarrow0 τ ↓ 0 . ISVIX coincides with the instantaneous diffusive variance in the absence of price jumps and otherwise includes a jump correction ([23 ] ).
Fractional Brownian Motion and Fractional Kernels
A standard fractional Brownian motion (fBM) with Hurst parameter H ∈ ( 0 , 1 ) H\in(0,1) H ∈ ( 0 , 1 ) is a centred Gaussian process B H B^H B H with covariance ([25 ] )
E Q [ B s H B t H ] = 1 2 ( s 2 H + t 2 H − ∣ t − s ∣ 2 H ) , s , t ≥ 0. (3) \mathbb{E}^{\mathbb{Q}}\left[B_s^H B_t^H\right] = \frac{1}{2}\left(s^{2H}+t^{2H}-|t-s|^{2H}\right), \qquad s,t\ge0. \tag*{(3)} E Q [ B s H B t H ] = 2 1 ( s 2 H + t 2 H − ∣ t − s ∣ 2 H ) , s , t ≥ 0. (3)
It has stationary increments and admits sample paths that are locally Hölder continuous of every order strictly less than H H H . The case H = 1 / 2 H=1/2 H = 1/2 is standard Brownian motion, whereas H < 1 / 2 H<1/2 H < 1/2 gives rougher paths.
We restrict attention to H ∈ ( 0 , 1 / 2 ] H\in(0,1/2] H ∈ ( 0 , 1/2 ] and set α = H + 1 / 2 \alpha=H+1/2 α = H + 1/2 . The corresponding fractional kernel is K ( t ) = t α − 1 / Γ ( α ) K(t)=t^{\alpha-1}/\Gamma(\alpha) K ( t ) = t α − 1 /Γ ( α ) for t > 0 t>0 t > 0 , where Γ \Gamma Γ is the gamma function. For a standard Brownian motion B B B , define
B ~ t H = ∫ 0 t K ( t − s ) d B s . (4) \widetilde{B}_t^H = \int_0^t K(t-s)\,\mathrm{d}B_s. \tag*{(4)} B t H = ∫ 0 t K ( t − s ) d B s . (4)
For H < 1 / 2 H<1/2 H < 1/2 , this Riemann–Liouville process differs from standard fBM: its increments are not stationary, although its paths have the same local Hölder regularity ([5 ] ). Since K K K is square-integrable on every bounded interval for α > 1 / 2 \alpha>1/2 α > 1/2 , the convolution is a Brownian Itô integral. At α = 1 \alpha=1 α = 1 , K ≡ 1 K\equiv1 K ≡ 1 and B ~ 1 / 2 = B \widetilde{B}^{1/2}=B B 1/2 = B .
Affine Volterra Models
We consider a class of affine Volterra models with diffusion and jump innovations ([3 , 6 ] ):
X t = X 0 + ∫ 0 t K ( t − s ) ∘ { μ ( X s ) d s + σ ( X s ) d B s + d J s } , (5) X_t = X_0 + \int_0^t \mathcal{K}(t-s) \circ\{\mu(X_s)\,\mathrm{d}s + \sigma(X_s)\,\mathrm{d}B_s + \mathrm{d}J_s\}, \tag*{(5)} X t = X 0 + ∫ 0 t K ( t − s ) ∘ { μ ( X s ) d s + σ ( X s ) d B s + d J s } , (5)
where X X X is a predictable d d d -dimensional process with locally integrable paths, X 0 X_0 X 0 is deterministic, B B B is an m m m -dimensional Brownian motion, K = ( K 1 , … , K d ) ⊤ \mathcal{K} = (K_1,\ldots,K_d)^\top K = ( K 1 , … , K d ) ⊤ consists of deterministic, locally square-integrable kernels, and ∘ \circ ∘ denotes componentwise multiplication. Here J J J is a compensated jump process whose jumps have conditional intensity λ y ( X s ) \lambda^y(X_s) λ y ( X s ) and sizes distributed as a random vector y ∈ R d y \in\mathbb{R}^d y ∈ R d with a fixed, state-independent law satisfying E Q [ ∥ y ∥ 2 ] < ∞ \mathbb{E}^{\mathbb{Q}}[\lVert y\rVert^2] < \infty E Q [∥ y ∥ 2 ] < ∞ ; its compensator subtracts λ y ( X s ) E Q [ y ] d s \lambda^y(X_s)\mathbb{E}^{\mathbb{Q}}[y]\,\mathrm{d}s λ y ( X s ) E Q [ y ] d s from the uncompensated jump increments. For singular kernels, Equation (5) is understood to hold Q ⊗ d t \mathbb{Q} \otimes\mathrm{d}t Q ⊗ d t -almost everywhere ([6 ] ).
The drift, diffusion covariance and jump intensity are affine on the state space:
μ ( x ) = b 0 + ∑ k = 1 d b k x k , σ ( x ) σ ( x ) ⊤ = H 0 + ∑ k = 1 d H k x k , λ y ( x ) = λ 0 y + ∑ k = 1 d λ k y x k , (6) \mu(x) = b_0 + \sum_{k=1}^{d} b_kx_k,\quad\sigma(x)\sigma(x)^\top= H_0 + \sum_{k=1}^{d} H_kx_k,\quad\lambda^y(x) = \lambda^y_0 + \sum_{k=1}^{d}\lambda^y_kx_k, \tag*{(6)} μ ( x ) = b 0 + k = 1 ∑ d b k x k , σ ( x ) σ ( x ) ⊤ = H 0 + k = 1 ∑ d H k x k , λ y ( x ) = λ 0 y + k = 1 ∑ d λ k y x k , (6)
where b k ∈ R d b_k \in\mathbb{R}^d b k ∈ R d , each H k ∈ R d × d H_k \in\mathbb{R}^{d \times d} H k ∈ R d × d is symmetric, and λ k y ∈ R \lambda^y_k \in\mathbb{R} λ k y ∈ R for k = 0 , … , d k = 0,\ldots,d k = 0 , … , d . We assume that the diffusion covariance is positive semidefinite and the jump intensity is nonnegative on the state space. Taking K i ≡ 1 K_i \equiv1 K i ≡ 1 for every component recovers an affine jump-diffusion, while setting λ y ≡ 0 \lambda^y \equiv0 λ y ≡ 0 gives the diffusion setting of Abi Jaber et al. [3 ] .
Whenever the expectation is absolutely convergent, write ϱ ( z ) = E Q [ e z ⊤ y ] \varrho(z) = \mathbb{E}^{\mathbb{Q}}[e^{z^\top y}] ϱ ( z ) = E Q [ e z ⊤ y ] and define the affine cumulant functions by
R k ( z ) = b k ⊤ z + 1 2 z ⊤ H k z + λ k y [ ϱ ( z ) − 1 − z ⊤ E Q [ y ] ] , k = 0 , … , d . (7) R_k(z) = b_k^\top z + \frac{1}{2}z^\top H_kz + \lambda^y_k\left[\varrho(z) - 1 - z^\top\mathbb{E}^{\mathbb{Q}}[y]\right],\quad k = 0,\ldots,d. \tag*{(7)} R k ( z ) = b k ⊤ z + 2 1 z ⊤ H k z + λ k y [ ϱ ( z ) − 1 − z ⊤ E Q [ y ] ] , k = 0 , … , d . (7)
For a horizon T > 0 T > 0 T > 0 and transform argument u ∈ C d u \in\mathbb{C}^d u ∈ C d , the associated Riccati–Volterra equation is
Ψ ( t , u ) = u ∘ K ( t ) + ∫ 0 t K ( t − s ) ∘ R ( Ψ ( s , u ) ) d s , 0 < t ≤ T , (8) \Psi(t,u) = u \circ\mathcal{K}(t) + \int_0^t \mathcal{K}(t-s) \circ R(\Psi(s,u))\,\mathrm{d}s,\quad0 < t \le T, \tag*{(8)} Ψ ( t , u ) = u ∘ K ( t ) + ∫ 0 t K ( t − s ) ∘ R ( Ψ ( s , u )) d s , 0 < t ≤ T , (8)
where R ( z ) = ( R 1 ( z ) , … , R d ( z ) ) ⊤ R(z) = (R_1(z),\ldots,R_d(z))^\top R ( z ) = ( R 1 ( z ) , … , R d ( z ) ) ⊤ .
Under suitable integrability and martingale conditions, solutions of Equation (8) yield initial-time exponential-affine transforms; see Abi Jaber et al. [3 ] and Bondi et al. [6 ] . The functions R 1 , … , R d R_1,\ldots,R_d R 1 , … , R d determine the Riccati–Volterra equation, whereas R 0 R_0 R 0 enters only the state-independent part of the transform exponent. At later valuation times, conditional transforms generally depend on the past through an adjusted forward curve rather than on X t X_t X t alone ([3 , 6 ] ).
Framework
In this section, we propose a rough ISVIX (RISVIX) model and derive the corresponding VIX term structure and VIX option pricing formulas. Throughout the theoretical analysis, time zero denotes the valuation origin.
Model Specification
We extend the ISVIX framework of Luo et al. [23 ] using the fractional kernel K ( t ) = t α − 1 / Γ ( α ) K(t) = t^{\alpha-1}/\Gamma(\alpha) K ( t ) = t α − 1 /Γ ( α ) introduced in Subsection 2.2. Fix a finite horizon T > 0 T > 0 T > 0 . Under the risk-neutral measure Q \mathbb{Q} Q , the RISVIX model is
V t = V 0 + ∫ 0 t K ( t − s ) { κ ( θ s − V s ) d s + σ V V s d B 1 , s + d J s y } , θ t = θ 0 + σ θ ∫ 0 t θ s d B 2 , s , 0 ≤ t ≤ T . (9) \begin{aligned}
V_t &= V_0 + \int_0^t K(t-s)\left\{\kappa(\theta_s - V_s)\,\mathrm{d}s + \sigma_V\sqrt{V_s}\,\mathrm{d}B_{1,s} + \mathrm{d}J_s^y\right\},\\
\theta_t &= \theta_0 + \sigma_\theta\int_0^t \sqrt{\theta_s}\,\mathrm{d}B_{2,s}, \qquad0 \le t \le T. \tag*{(9)}
\end{aligned} V t θ t = V 0 + ∫ 0 t K ( t − s ) { κ ( θ s − V s ) d s + σ V V s d B 1 , s + d J s y } , = θ 0 + σ θ ∫ 0 t θ s d B 2 , s , 0 ≤ t ≤ T . (9)
Here V 0 , θ 0 ≥ 0 V_0, \theta_0 \ge0 V 0 , θ 0 ≥ 0 are F 0 \mathcal{F}_0 F 0 -measurable initial states satisfying E Q [ V 0 2 + θ 0 2 ] < ∞ \mathbb{E}^{\mathbb{Q}}[V_0^2 + \theta_0^2] < \infty E Q [ V 0 2 + θ 0 2 ] < ∞ , and B 1 B_1 B 1 and B 2 B_2 B 2 are independent Brownian motions. Let N y N^y N y be a Poisson process with intensity λ y \lambda^y λ y , and let { Y n } n ≥ 1 \{Y_n\}_{n\ge1} { Y n } n ≥ 1 be an independent sequence of exponentially distributed random variables with mean μ y \mu_y μ y . We assume that B 1 B_1 B 1 , B 2 B_2 B 2 , N y N^y N y and { Y n } n ≥ 1 \{Y_n\}_{n\ge1} { Y n } n ≥ 1 are mutually independent and independent of F 0 \mathcal{F}_0 F 0 , and define J t y = ∑ n = 1 N t y Y n − λ y μ y t J_t^y = \sum_{n=1}^{N_t^y}Y_n - \lambda^y\mu_y t J t y = ∑ n = 1 N t y Y n − λ y μ y t . Thus, J y J^y J y is a square-integrable martingale, and the jump-size moment generating function is ϱ ( z ) = ( 1 − μ y z ) − 1 \varrho(z) = (1-\mu_y z)^{-1} ϱ ( z ) = ( 1 − μ y z ) − 1 for Re ( z ) < μ y − 1 \operatorname{Re}(z) < \mu_y^{-1} Re ( z ) < μ y − 1 .
The parameters satisfy κ , μ y > 0 \kappa, \mu_y > 0 κ , μ y > 0 , σ V , σ θ , λ y ≥ 0 \sigma_V, \sigma_\theta, \lambda^y \ge0 σ V , σ θ , λ y ≥ 0 and α ∈ ( 1 / 2 , 1 ) \alpha\in(1/2,1) α ∈ ( 1/2 , 1 ) . The parameter κ \kappa κ governs mean reversion towards the stochastic level θ \theta θ , while σ V \sigma_V σ V and σ θ \sigma_\theta σ θ govern innovations to ISVIX and its mean level, respectively. The fractional kernel introduces memory into the ISVIX dynamics; formally setting α = 1 \alpha= 1 α = 1 gives K ≡ 1 K \equiv1 K ≡ 1 and a Markovian jump-diffusion.
On [ 0 , T ] [0,T] [ 0 , T ] , we assume that Equation (9) admits a nonnegative weak solution ( V , θ ) (V,\theta) ( V , θ ) . The process θ \theta θ is continuous and nonnegative, while V V V is predictable with ∫ 0 T V s d s < ∞ \int_0^T V_s\,\mathrm{d}s < \infty ∫ 0 T V s d s < ∞ almost surely; both V ≥ 0 V \ge0 V ≥ 0 and the first equation in Equation (9) hold Q ⊗ d t \mathbb{Q} \otimes\mathrm{d}t Q ⊗ d t -almost everywhere. Under the initial second-moment condition above, Lemma A.2 gives E Q [ ∫ 0 T ( V t 2 + θ t 2 ) d t ] < ∞ \mathbb{E}^{\mathbb{Q}}\left[\int_0^T (V_t^2+\theta_t^2)\,\mathrm{d}t\right] < \infty E Q [ ∫ 0 T ( V t 2 + θ t 2 ) d t ] < ∞ and establishes that θ \theta θ is a square-integrable martingale. General weak-existence results for affine stochastic Volterra equations can be found in Abi Jaber [1 ] . The following proposition gives the conditional mean of ISVIX.
Proposition 3.1. Suppose that Equation (9) admits the finite-horizon solution specified in Appendix § 8 .§ 2 on [ 0 , T ] [0,T] [ 0 , T ] . Then, the following identity holds almost surely:
E Q [ V s ] = θ 0 + ( V 0 − θ 0 ) E α , 1 ( − κ s α ) , (10) \mathbb{E}^{\mathbb{Q}}[V_s] = \theta_0 + (V_0-\theta_0)E_{\alpha,1}(-\kappa s^\alpha), \tag*{(10)} E Q [ V s ] = θ 0 + ( V 0 − θ 0 ) E α , 1 ( − κ s α ) , (10) where E a , b ( z ) = ∑ n = 0 ∞ z n / Γ ( a n + b ) E_{a,b}(z) = \sum_{n=0}^{\infty} z^n/\Gamma(an+b) E a , b ( z ) = ∑ n = 0 ∞ z n /Γ ( an + b ) is the two-parameter Mittag–Leffler function. The right-hand side gives the unique continuous representative of the conditional mean on [ 0 , T ] [0,T] [ 0 , T ] .
Remark. The conditional mean in Equation (10) depends only on ( α , κ , V 0 , θ 0 ) (\alpha,\kappa,V_0,\theta_0) ( α , κ , V 0 , θ 0 ) . The Brownian and compensated jump terms have zero conditional mean, while E 0 Q [ θ s ] = θ 0 \mathbb{E}^{\mathbb{Q}}_0[\theta_s] = \theta_0 E 0 Q [ θ s ] = θ 0 . Hence, although ( σ V , σ θ , λ y , μ y ) (\sigma_V, \sigma_\theta, \lambda^y, \mu_y) ( σ V , σ θ , λ y , μ y ) do not enter the first moment, they shape the conditional distribution of future ISVIX and therefore affect VIX option prices.
VIX Term Structure
Recall from Equation (2) that, for t ≥ 0 t \ge0 t ≥ 0 and τ > 0 \tau> 0 τ > 0 , the squared VIX is V I X t , τ 2 = τ − 1 E t Q [ ∫ t t + τ V s d s ] \mathrm{VIX}_{t,\tau}^2 = \tau^{-1}\mathbb{E}^{\mathbb{Q}}_t\left[\int_t^{t+\tau} V_s\,\mathrm{d}s\right] VIX t , τ 2 = τ − 1 E t Q [ ∫ t t + τ V s d s ] . The following theorem gives its representation under the RISVIX model.
Theorem 3.1. Suppose that Equation (9) admits the finite-horizon solution specified in Appendix § 8 .§ 2 on [ 0 , T ] [0,T] [ 0 , T ] . Then, for each fixed t ∈ [ 0 , T ) t \in[0,T) t ∈ [ 0 , T ) and τ ∈ ( 0 , T − t ] \tau\in(0,T-t] τ ∈ ( 0 , T − t ] , the following identity holds almost surely:
V I X t , τ 2 = w ( τ ) V 0 + [ 1 − w ( τ ) ] θ t + ∫ 0 t a τ ( t − s ) { κ ( θ s − V 0 ) d s + σ V V s d B 1 , s + d J s y } . (11) \mathrm{VIX}_{t,\tau}^2 = w(\tau)V_0 + [1-w(\tau)]\theta_t + \int_0^t a_\tau(t-s)\left\{\kappa(\theta_s-V_0)\,\mathrm{d}s+\sigma_V\sqrt{V_s}\,\mathrm{d}B_{1,s}+\mathrm{d}J_s^y\right\}. \tag*{(11)} VIX t , τ 2 = w ( τ ) V 0 + [ 1 − w ( τ )] θ t + ∫ 0 t a τ ( t − s ) { κ ( θ s − V 0 ) d s + σ V V s d B 1 , s + d J s y } . (11) where w ( τ ) = E α , 2 ( − κ τ α ) w(\tau)=E_{\alpha,2}(-\kappa\tau^\alpha) w ( τ ) = E α , 2 ( − κ τ α ) and a τ ( s ) = [ E α , 1 ( − κ s α ) − E α , 1 ( − κ ( s + τ ) α ) ] / ( κ τ ) a_\tau(s)=[E_{\alpha,1}(-\kappa s^\alpha)-E_{\alpha,1}(-\kappa(s+\tau)^\alpha)]/(\kappa\tau) a τ ( s ) = [ E α , 1 ( − κ s α ) − E α , 1 ( − κ ( s + τ ) α )] / ( κ τ ) .
Remark. At the valuation origin, Equation (11) reduces to V I X 0 , τ 2 = V 0 w ( τ ) + θ 0 [ 1 − w ( τ ) ] \mathrm{VIX}_{0,\tau}^2=V_0w(\tau)+\theta_0[1-w(\tau)] VIX 0 , τ 2 = V 0 w ( τ ) + θ 0 [ 1 − w ( τ )] . The weight w ( τ ) w(\tau) w ( τ ) decreases from one to zero as τ \tau τ increases. Hence, over any horizon on which the finite-horizon solution is defined, the time-zero squared VIX curve moves monotonically from V 0 V_0 V 0 towards θ 0 \theta_0 θ 0 : it is upward sloping when V 0 < θ 0 V_0<\theta_0 V 0 < θ 0 , downward sloping when V 0 > θ 0 V_0>\theta_0 V 0 > θ 0 , and flat when the two levels coincide. The time-zero term structure depends only on ( α , κ , V 0 , θ 0 ) (\alpha,\kappa,V_0,\theta_0) ( α , κ , V 0 , θ 0 ) . The diffusion and jump parameters do not enter this curve, but they affect the distribution of future VIX through the history term in Equation (11) . This separation motivates the two-stage calibration used below.
VIX Option Pricing
Consider a European VIX call valued at time zero, with expiry T > 0 T>0 T > 0 , strike K > 0 K>0 K > 0 , and variance horizon τ 0 > 0 \tau_0>0 τ 0 > 0 . By Theorem 3.1, its terminal squared VIX is given by Equation (11) evaluated at t = T t=T t = T and τ = τ 0 \tau=\tau_0 τ = τ 0 . We price the option through the Fourier–Laplace transform
f ( u ; 0 , T ) = E 0 Q [ exp { u V I X T , τ 0 2 } ] . (12) f(u;0,T)=\mathbb{E}^{\mathbb{Q}}_0\left[\exp\left\{u\mathrm{VIX}_{T,\tau_0}^2\right\}\right]. \tag*{(12)} f ( u ; 0 , T ) = E 0 Q [ exp { u VIX T , τ 0 2 } ] . (12)
Let R κ R_\kappa R κ denote the resolvent of κ K \kappa K κ K , satisfying R κ + κ K ∗ R κ = κ K R_\kappa+\kappa K*R_\kappa=\kappa K R κ + κ K ∗ R κ = κ K . By Lemma A.1 ,
R κ ( s ) = κ s α − 1 E α , α ( − κ s α ) , a τ 0 ( s ) = 1 κ τ 0 ∫ s s + τ 0 R κ ( v ) d v . (13) R_\kappa(s)=\kappa s^{\alpha-1}E_{\alpha,\alpha}(-\kappa s^\alpha),\qquad a_{\tau_0}(s)=\frac{1}{\kappa\tau_0}\int_s^{s+\tau_0}R_\kappa(v)\,\mathrm{d}v. \tag*{(13)} R κ ( s ) = κ s α − 1 E α , α ( − κ s α ) , a τ 0 ( s ) = κ τ 0 1 ∫ s s + τ 0 R κ ( v ) d v . (13)
Suppose that Equation (9) admits the finite-horizon solution specified in Appendix § 8 .§ 2 on [ 0 , T + τ 0 ] [0,T+\tau_0] [ 0 , T + τ 0 ] .
Theorem 3.2. *There exists η ˉ T , τ 0 > 0 \bar{\eta}_{T,\tau_0}>0 η ˉ T , τ 0 > 0 such that, for every u ∈ C u\in\mathbb{C} u ∈ C with 0 ≤ Re ( u ) < η ˉ T , τ 0 0\le\operatorname{Re}(u)<\bar{\eta}_{T,\tau_0} 0 ≤ Re ( u ) < η ˉ T , τ 0 , the
transform satisfies almost surely
f ( u ; 0 , T ) = exp { ϕ ( T ) + χ 1 ( T ) V 0 + ψ 2 ( T ) θ 0 } , (14) f(u;0,T)=\exp\{\phi(T)+\chi_1(T)V_0+\psi_2(T)\theta_0\}, \tag*{(14)} f ( u ; 0 , T ) = exp { ϕ ( T ) + χ 1 ( T ) V 0 + ψ 2 ( T ) θ 0 } , (14)
where the coefficients ψ 1 \psi_1 ψ 1 , ψ 2 \psi_2 ψ 2 , χ 1 \chi_1 χ 1 and ϕ \phi ϕ are the unique continuous functions on [ 0 , T ] [0,T] [ 0 , T ] satisfying
ψ 1 ( s ) = u a τ 0 ( s ) + σ V 2 2 κ ∫ 0 s R κ ( s − v ) ψ 1 ( v ) 2 d v , ψ 2 ( s ) = u [ 1 − w ( τ 0 ) ] + ∫ 0 s [ κ ψ 1 ( v ) + σ θ 2 2 ψ 2 ( v ) 2 ] d v , χ 1 ( s ) = u w ( τ 0 ) + ∫ 0 s [ − κ ψ 1 ( v ) + σ V 2 2 ψ 1 ( v ) 2 ] d v , ϕ ( s ) = λ y ∫ 0 s [ ρ ( ψ 1 ( v ) ) − 1 − μ y ψ 1 ( v ) ] d v . (15) \begin{aligned}
\psi_1(s)&=ua_{\tau_0}(s)+\frac{\sigma_V^2}{2\kappa}\int_0^s R_\kappa(s-v)\psi_1(v)^2\,\mathrm{d}v,\\
\psi_2(s)&=u[1-w(\tau_0)]+\int_0^s\left[\kappa\psi_1(v)+\frac{\sigma_\theta^2}{2}\psi_2(v)^2\right]\,\mathrm{d}v,\\
\chi_1(s)&=uw(\tau_0)+\int_0^s\left[-\kappa\psi_1(v)+\frac{\sigma_V^2}{2}\psi_1(v)^2\right]\,\mathrm{d}v,\\
\phi(s)&=\lambda^y\int_0^s[\rho(\psi_1(v))-1-\mu_y\psi_1(v)]\,\mathrm{d}v. \tag*{(15)}
\end{aligned} ψ 1 ( s ) ψ 2 ( s ) χ 1 ( s ) ϕ ( s ) = u a τ 0 ( s ) + 2 κ σ V 2 ∫ 0 s R κ ( s − v ) ψ 1 ( v ) 2 d v , = u [ 1 − w ( τ 0 )] + ∫ 0 s [ κ ψ 1 ( v ) + 2 σ θ 2 ψ 2 ( v ) 2 ] d v , = u w ( τ 0 ) + ∫ 0 s [ − κ ψ 1 ( v ) + 2 σ V 2 ψ 1 ( v ) 2 ] d v , = λ y ∫ 0 s [ ρ ( ψ 1 ( v )) − 1 − μ y ψ 1 ( v )] d v . (15)
The solution satisfies sup 0 ≤ s ≤ T Re ψ 1 ( s ) < μ y − 1 \sup_{0\le s\le T}\operatorname{Re}\psi_1(s)<\mu_y^{-1} sup 0 ≤ s ≤ T Re ψ 1 ( s ) < μ y − 1 , so the jump transform ρ \rho ρ is well defined.
Remark. For each transform argument u u u , computing the coefficients requires solving a nonlinear Volterra equation for ψ 1 \psi_1 ψ 1 and a scalar Riccati ODE for ψ 2 \psi_2 ψ 2 ; χ 1 \chi_1 χ 1 and ϕ \phi ϕ are then obtained by integration. The variance-window averaging enters through a τ 0 a_{\tau_0} a τ 0 , while fractional memory enters through both a τ 0 a_{\tau_0} a τ 0 and R κ R_\kappa R κ . Formally setting α = 1 \alpha=1 α = 1 gives a τ 0 ( s ) = w ( τ 0 ) e − κ s a_{\tau_0}(s)=w(\tau_0)e^{-\kappa s} a τ 0 ( s ) = w ( τ 0 ) e − κ s and χ 1 = ψ 1 \chi_1=\psi_1 χ 1 = ψ 1 , reducing Equation (15) to the corresponding Markovian Riccati system.
To recover option prices from Equation (14) , we view the call payoff as a function of squared VIX, g ( x ) = ( x − K ) + g(x)=(\sqrt{x}-K)^+ g ( x ) = ( x − K ) + for x ≥ 0 x\ge0 x ≥ 0 . Its Laplace transform is ∫ 0 ∞ e − z x g ( x ) d x = π erfc ( K z ) / ( 2 z 3 / 2 ) \int_0^\infty e^{-zx}g(x)\,\mathrm{d}x=\sqrt{\pi}\,\operatorname{erfc}(K\sqrt{z})/(2z^{3/2}) ∫ 0 ∞ e − z x g ( x ) d x = π erfc ( K z ) / ( 2 z 3/2 ) for Re ( z ) > 0 \operatorname{Re}(z)>0 Re ( z ) > 0 . Theorem 3.2 guarantees that the VIX transform is finite on every vertical contour u = η + i ξ u=\eta+i\xi u = η + i ξ with η ∈ ( 0 , η ˉ T , τ 0 ) \eta\in(0,\bar{\eta}_{T,\tau_0}) η ∈ ( 0 , η ˉ T , τ 0 ) . Hence any such η \eta η can be used as the damping parameter in the inversion below; the resulting option value does not depend on its particular choice.
Corollary 3.1. Assume a constant risk-free rate r r r and let η ∈ ( 0 , η ˉ T , τ 0 ) \eta\in(0,\bar{\eta}_{T,\tau_0}) η ∈ ( 0 , η ˉ T , τ 0 ) . Then the time-zero price of a European VIX call with expiry T T T and strike K K K is
c ( 0 , T ; K ) = e − r T E Q [ ( V I X T , τ 0 − K ) + ] = e − r T 2 π ∫ 0 ∞ Re [ erfc ( K η + i ξ ) ( η + i ξ ) 3 / 2 f ( η + i ξ ; 0 , T ) ] d ξ . (16) \begin{aligned}
c(0,T;K)&=e^{-rT}\mathbb{E}_{\mathbb{Q}}\left[(\mathrm{VIX}_{T,\tau_0}-K)^+\right]\\
&=\frac{e^{-rT}}{2\sqrt{\pi}}\int_0^\infty\operatorname{Re}\left[\frac{\operatorname{erfc}(K\sqrt{\eta+i\xi})}{(\eta+i\xi)^{3/2}}f(\eta+i\xi;0,T)\right]\,\mathrm{d}\xi. \tag*{(16)}
\end{aligned} c ( 0 , T ; K ) = e − r T E Q [ ( VIX T , τ 0 − K ) + ] = 2 π e − r T ∫ 0 ∞ Re [ ( η + i ξ ) 3/2 erfc ( K η + i ξ ) f ( η + i ξ ; 0 , T ) ] d ξ . (16)
Here erfc ( z ) = 1 − 2 π ∫ 0 z e − v 2 d v \operatorname{erfc}(z)=1-\frac{2}{\sqrt{\pi}}\int_0^z e^{-v^2}\,\mathrm{d}v erfc ( z ) = 1 − π 2 ∫ 0 z e − v 2 d v , and the square root and complex power are taken on the principal branch.
Thus, VIX option pricing under RISVIX reduces to solving one nonlinear Volterra equation for each transform argument and evaluating a one-dimensional numerical integral.
Calibration Strategy
The parameter separation leads naturally to a two-step procedure, following the estimation–calibration structure of Luo et al. [23 ] . The first step estimates the common term-structure parameters and date-specific initial states from VIX indices, while the second step calibrates the remaining distributional parameters to VIX option prices.
Step 1: Estimation from the VIX term structure Let T \mathcal{T} T denote the set of calibration dates and let τ 1 , … , τ J \tau_1,\ldots,\tau_J τ 1 , … , τ J denote the observed VIX horizons. At each calendar date t 0 ∈ T t_0 \in\mathcal{T} t 0 ∈ T , we treat t 0 t_0 t 0 as the valuation origin of a local model and apply the time-zero relation in Theorem 3.1, ( V I X t 0 , τ j ) 2 = V t 0 w ( τ j ) + θ t 0 [ 1 − w ( τ j ) ] (\mathrm{VIX}_{t_0,\tau_j})^2 = V_{t_0}w(\tau_j) + \theta_{t_0}[1-w(\tau_j)] ( VIX t 0 , τ j ) 2 = V t 0 w ( τ j ) + θ t 0 [ 1 − w ( τ j )] . Here ( V t 0 , θ t 0 ) (V_{t_0},\theta_{t_0}) ( V t 0 , θ t 0 ) are the date-specific initial states of the local model; the subscript t 0 t_0 t 0 indexes the calendar valuation date. We estimate the common parameters ( α , κ ) (\alpha,\kappa) ( α , κ ) jointly with these initial states by
( α ^ , κ ^ , { V ^ t 0 , θ ^ t 0 } t 0 ∈ T ) = arg min α ∈ ( 1 / 2 , 1 ] , κ > 0 V t 0 ≥ 0 , θ t 0 ≥ 0 , t 0 ∈ T ∑ t 0 ∈ T ∑ j = 1 J [ V I X t 0 , τ j m k t − V I X t 0 , τ j m o d e l ( α , κ , V t 0 , θ t 0 ) ] 2 . (17) \left(\widehat{\alpha},\widehat{\kappa},\{\widehat{V}_{t_0},\widehat{\theta}_{t_0}\}_{t_0\in\mathcal{T}}\right)
=
\underset{\substack{\alpha\in(1/2,1],\ \kappa>0\\ V_{t_0}\geq0,\ \theta_{t_0}\geq0,\ t_0\in\mathcal{T}}}{\operatorname{arg\,min}}
\sum_{t_0\in\mathcal{T}}\sum_{j=1}^{J}
\left[
\mathrm{VIX}^{\mathrm{mkt}}_{t_0,\tau_j}
-
\mathrm{VIX}^{\mathrm{model}}_{t_0,\tau_j}(\alpha,\kappa,V_{t_0},\theta_{t_0})
\right]^2 . \tag*{(17)} ( α , κ , { V t 0 , θ t 0 } t 0 ∈ T ) = α ∈ ( 1/2 , 1 ] , κ > 0 V t 0 ≥ 0 , θ t 0 ≥ 0 , t 0 ∈ T arg min t 0 ∈ T ∑ j = 1 ∑ J [ VIX t 0 , τ j mkt − VIX t 0 , τ j model ( α , κ , V t 0 , θ t 0 ) ] 2 . (17)
For any fixed ( α , κ ) (\alpha,\kappa) ( α , κ ) , the date-specific initial-state estimates are separable across dates and are profiled out subject to nonnegativity. Our LZZ benchmark applies the same procedure while imposing α = 1 \alpha=1 α = 1 .
Step 2: Calibration to VIX option prices Conditional on the Step 1 estimates, we calibrate ζ = ( σ V , σ θ , λ y , μ y ) ⊤ \zeta=(\sigma_V,\sigma_\theta,\lambda^y,\mu_y)^\top ζ = ( σ V , σ θ , λ y , μ y ) ⊤ using VIX call prices. Let T t r a i n \mathcal{T}_{\mathrm{train}} T train denote the set of option-calibration dates. For each t 0 ∈ T t r a i n t_0\in\mathcal{T}_{\mathrm{train}} t 0 ∈ T train , let M t 0 \mathcal{M}_{t_0} M t 0 denote the observed times to expiration and let K t 0 , T \mathcal{K}_{t_0,T} K t 0 , T denote the strikes available for T ∈ M t 0 T\in\mathcal{M}_{t_0} T ∈ M t 0 . The local model at t 0 t_0 t 0 is initialized with ( V ^ t 0 , θ ^ t 0 ) (\widehat{V}_{t_0},\widehat{\theta}_{t_0}) ( V t 0 , θ t 0 ) . The option maturity T T T is measured from t 0 t_0 t 0 , and the VIX variance horizon τ 0 \tau_0 τ 0 is fixed at 30 days. The second-step estimator is
ζ ^ = arg min ζ ∈ P ∑ t 0 ∈ T t r a i n ∑ T ∈ M t 0 ∑ K ∈ K t 0 , T [ C t 0 m k t ( T , K ) − C t 0 m o d e l ( T , K ; ζ , α ^ , κ ^ , V ^ t 0 , θ ^ t 0 ) ] 2 . (18) \widehat{\zeta}
=
\underset{\zeta\in\mathcal{P}}{\operatorname{arg\,min}}
\sum_{t_0\in\mathcal{T}_{\mathrm{train}}}
\sum_{T\in\mathcal{M}_{t_0}}
\sum_{K\in\mathcal{K}_{t_0,T}}
\left[
C^{\mathrm{mkt}}_{t_0}(T,K)
-
C^{\mathrm{model}}_{t_0}\left(T,K;\zeta,\widehat{\alpha},\widehat{\kappa},\widehat{V}_{t_0},\widehat{\theta}_{t_0}\right)
\right]^2 . \tag*{(18)} ζ = ζ ∈ P arg min t 0 ∈ T train ∑ T ∈ M t 0 ∑ K ∈ K t 0 , T ∑ [ C t 0 mkt ( T , K ) − C t 0 model ( T , K ; ζ , α , κ , V t 0 , θ t 0 ) ] 2 . (18)
Here C t 0 m k t ( T , K ) C^{\mathrm{mkt}}_{t_0}(T,K) C t 0 mkt ( T , K ) is the bid–ask midpoint, and C t 0 m o d e l ( T , K ; ⋅ ) C^{\mathrm{model}}_{t_0}(T,K;\cdot) C t 0 model ( T , K ; ⋅ ) is the price c ( 0 , T ; K ) c(0,T;K) c ( 0 , T ; K ) in Equation (16) , evaluated using the Step 1 estimates and the candidate value of ζ \zeta ζ . The parameters in ζ \zeta ζ are common across all valuation dates, expirations and strikes.
For each valuation date t 0 t_0 t 0 , we re-index that date as local time zero. Under this re-indexing, the relevant local model interval for an option with time to expiration T T T is [ 0 , T + τ 0 ] [0,T+\tau_0] [ 0 , T + τ 0 ] , corresponding to the calendar interval [ t 0 , t 0 + T + τ 0 ] [t_0,t_0+T+\tau_0] [ t 0 , t 0 + T + τ 0 ] . Such re-indexing is standard in rough-volatility pricing; depending on the model, conditioning at the new origin is encoded through an initial forward curve or another model-specific initial input ([11 , 16 ] ). In our empirical implementation, each t 0 t_0 t 0 defines a local RISVIX valuation problem initialized by the date-specific initial states ( V ^ t 0 , θ ^ t 0 ) (\widehat{V}_{t_0},\widehat{\theta}_{t_0}) ( V t 0 , θ t 0 ) inferred from the contemporaneous VIX term structure.
Empirical Study
Data
We obtain daily VIX term-structure observations from CBOE and VIX option quotes from OptionMetrics. The term structure consists of the closing values of VIX9D, VIX, VIX3M, VIX6M and VIX1Y, corresponding to volatility horizons of nine days, 30 days, three months, six months and one year, respectively. The sample runs from January 2, 2018 to June 30, 2025. Restricting attention to dates on which all five indices are observed gives a balanced panel of 1,883 trading days.
The VIX option data also run from January 2, 2018 to June 30, 2025 and contain trading dates, expiration dates, strikes, option types, closing bid and ask quotes, and trading volume. We use the bid–ask midpoint as the observed option price, with option prices and strikes expressed in VIX index points. We retain call options with 30–126 calendar days to expiration and spot moneyness K / V I X t ∈ [ 0.7 , 1.3 ] K/\mathrm{VIX}_{t} \in[0.7, 1.3] K / VIX t ∈ [ 0.7 , 1.3 ] , where K K K denotes the strike and V I X t \mathrm{VIX}_{t} VIX t is the contemporaneous closing value of the 30-day VIX index. We further require bid and ask quotes of at least 0.025 points, positive trading volume and open interest, a midquote of at least 0.50 points, a nonnegative bid–ask spread no greater than 10% of the midquote. These filters leave 61,994 quote–day observations: 47,451 in the training period from January 2018 to December 2023 and 14,543 in the test period from January 2024 to June 2025.
Table 1 reports summary statistics for the two datasets. Panel A shows an upward-sloping average term structure, with mean index levels increasing from 19.29 for VIX9D to 23.55 for VIX1Y, while the corresponding standard deviations decrease from 9.50 to 4.92. Panel B reports an average option midquote of 4.28 points. Longer-dated calls have higher average prices, whereas shorter-dated calls exhibit greater positive skewness and kurtosis.
Table 1. Summary statistics of VIX indices and options. Panel A reports daily closing index levels. Panel B reports VIX call-option midquotes, grouped by calendar days to expiration.
Figure 1 illustrates the historical evolution and term-structure variation of the VIX indices. Panel (a) shows substantial co-movement across horizons, with shorter-horizon indices exhibiting sharper fluctuations. Panel (b) shows that the term structure is typically upward sloping during low-volatility periods but can invert during volatility spikes, most notably in early 2020. Thus, the data exhibit substantial variation in both the level and shape of the VIX term structure.
Figure 1. Historical dynamics and term structure of VIX indices. Panel (a) plots the available daily closing values of VIX9D, VIX, VIX3M, VIX6M and VIX1Y on a common scale. Panel (b) displays the balanced term-structure panel used in the empirical analysis.
Term-Structure Calibration
We compare RISVIX with the benchmark of Luo et al. [23 ] , denoted LZZ (α = 1 \alpha= 1 α = 1 ), using the VIX-level least-squares criterion in § 3 . The full-sample comparison also includes the rough Heston term-structure model of Ye et al. (2026), denoted YFK. The full calibration sample in this subsection covers January 2018–June 2025 and contains 1,883 trading days and 9,415 horizon–date observations. We also examine a separate calibration of RISVIX and LZZ over January–June 2020. For the subsequent option pricing, the first-step inputs are estimated separately using January 2018–December 2023, and January 2024–June 2025 is reserved as the test period.
For RISVIX, the common parameters (α , κ \alpha, \kappa α , κ ) and date-specific nonnegative initial states (V t , θ t V_t, \theta_t V t , θ t ) are updated alternately. Given (α , κ \alpha, \kappa α , κ ), we fit (V t , θ t V_t, \theta_t V t , θ t ) separately for each date; holding those states fixed, we then update (α , κ \alpha, \kappa α , κ ) using all calibration dates. The states are re-estimated at the updated parameters, and this cycle is repeated until the stopping tolerance or iteration limit is reached. LZZ imposes α = 1 \alpha= 1 α = 1 and uses the same VIX-level least-squares criterion.
Full-Sample Fit
Table 2 reports the parameter estimates and calibration errors for January 2018–June 2025. We recalibrate YFK (Ye et al. 2026) to the same VIX term-structure observations using a 20-factor
Table 2. Full-sample term-structure calibration. Panel A reports parameter estimates for January 2018–June 2025 (1,883 trading days). LZZ imposes α = 1 \alpha= 1 α = 1 , and κ \kappa κ corresponds to the parameter denoted by λ \lambda λ in YFK. Panel B reports MSE for the three models on the same observations and the percentage reductions achieved by RISVIX relative to LZZ and YFK. Positive reductions favour RISVIX.
Markovian lift of the fractional volatility kernel.
In addition to VIX levels, we assess term-structure slope and curvature using divided differences. Let y t , j y_{t,j} y t , j denote the VIX index level on date t t t at horizon τ j \tau_j τ j , where τ 1 < ⋯ < τ 5 \tau_1 < \cdots< \tau_5 τ 1 < ⋯ < τ 5 are expressed in years. We approximate slope and curvature by
s t , j = y t , j + 1 − y t , j τ j + 1 − τ j , j = 1 , … , 4 , c t , j = 2 ( s t , j + 1 − s t , j ) τ j + 2 − τ j , j = 1 , … , 3. (19) \begin{aligned}
s_{t,j} &= \frac{y_{t,j+1} - y_{t,j}}{\tau_{j+1} - \tau_j}, && j = 1,\ldots,4,\\
c_{t,j} &= \frac{2(s_{t,j+1} - s_{t,j})}{\tau_{j+2} - \tau_j}, && j = 1,\ldots,3. \tag*{(19)}
\end{aligned} s t , j c t , j = τ j + 1 − τ j y t , j + 1 − y t , j , = τ j + 2 − τ j 2 ( s t , j + 1 − s t , j ) , j = 1 , … , 4 , j = 1 , … , 3. (19)
These definitions account for unequal horizon spacings and are applied to both observed and model-implied VIX levels.
For each measure, MSE is the average squared difference between observed and model-implied values over the relevant dates and horizons. Each model is evaluated on 9,415 level errors, 7,532 slope errors and 5,649 curvature errors. The percentage MSE reduction relative to benchmark B ∈ { L Z Z , Y F K } B \in\{\mathrm{LZZ},\mathrm{YFK}\} B ∈ { LZZ , YFK } is 100 ( 1 − M S E R I S V I X / M S E B ) 100(1-\mathrm{MSE}_{\mathrm{RISVIX}}/\mathrm{MSE}_B) 100 ( 1 − MSE RISVIX / MSE B ) ; positive values favour RISVIX.
RISVIX yields α ^ = 0.8942 \hat{\alpha} = 0.8942 α ^ = 0.8942 (H ^ = 0.3942 \hat{H} = 0.3942 H ^ = 0.3942 ), compared with α ^ = 0.6997 \hat{\alpha} = 0.6997 α ^ = 0.6997 (H ^ = 0.1997 \hat{H} = 0.1997 H ^ = 0.1997 ) for YFK, and achieves the lowest MSE across all three measures. Relative to LZZ, RISVIX reduces level MSE by 10.42%, with larger reductions of 17.21% and 13.18% for slope and curvature, respectively. The corresponding reductions relative to YFK are 88.39%, 57.86% and 27.34%.
State-Dependent Fit
The relative gains of RISVIX vary substantially across market conditions. We split the full calibration sample into dates with large VIX increases and the remaining dates, keeping the full-sample model parameters fixed. A shock date is defined as a daily VIX increase at or above the 90th percentile of the calibration-sample changes. Figure 2 reports percentage MSE reductions for levels, slopes and curvature in the full sample and in the two groups. On shock dates, the respective reductions are 28.73%, 26.34% and 20.60%, compared with 1.18%, 11.71% and 9.31% on the remaining dates. RISVIX therefore improves all three diagnostics in both groups, with larger relative improvements when VIX rises sharply.
Figure 2. Term-structure fit by VIX shocks. Bars report percentage MSE reductions for RISVIX relative to LZZ for each diagnostic and group. The sample runs from January 2, 2018 to June 30, 2025.
We next examine the initial COVID-19 period from January 2 to June 30, 2020, which contains 125 trading days and 625 horizon–date observations. Both models are separately estimated using this period.
Table 3 shows a larger departure from the Markovian restriction during the COVID-19 period, with α ^ = 0.7051 \hat{\alpha} = 0.7051 α ^ = 0.7051 (H ^ = 0.2051 \hat{H} = 0.2051 H ^ = 0.2051 ). The level MSE falls from 3.0105 for LZZ to 2.4808 for RISVIX, a reduction of 17.60%. The corresponding reductions in slope and curvature MSE are 41.14% and 46.59%. Relative to the full calibration sample, the improvement is particularly pronounced for the shape diagnostics.
Table 3. Term-structure calibration during the COVID-19 period. Panel A reports parameter values for January 2–June 30, 2020 (125 trading days); LZZ imposes α = 1 \alpha= 1 α = 1 . Panel B reports MSE on identical observations for the two models and its percentage reduction relative to LZZ. Positive reductions favour RISVIX.
Fit across Horizons
Figure 3 decomposes the percentage reduction in level MSE across VIX horizons. Panel (a) uses the full January 2018–June 2025 sample. The largest reduction is 33.30% at nine days, followed by 8.26% at three months and 7.34% at six months. At 30 days and one year, MSE is higher than LZZ, with reductions of −1.05% and −18.83%, respectively.
Figure 3. Term-structure fit across horizons. Bars report horizon-specific percentage reductions in level MSE for RISVIX relative to LZZ. Panel (a) uses the full January 2018–June 2025 sample; Panel (b) uses models separately calibrated over January–June 2020. Each panel has one bar per horizon. Green bars indicate positive reductions and red bars indicate negative reductions. Dashed lines show the corresponding aggregate level-MSE reductions. The panels share the same vertical scale.
Panel (b) uses the separate COVID-19 calibrations and displays one bar per horizon. The largest reduction, 47.69%, occurs at the 30-day horizon, followed by 23.47% at nine days and 18.64% at six months. The one-year MSE declines by 3.20%, while the three-month MSE increases by 3.33%. The dashed lines mark the overall level-MSE reductions of 10.42% and 17.60% for the two calibration samples.
VIX Option Pricing
We next evaluate whether the term-structure flexibility of RISVIX is associated with more accurate out-of-sample VIX option prices. The training period is January 2018–December 2023, and the test period is January 2024–June 2025. The global parameter sets reported below are common to all dates, expirations and strikes and are held fixed throughout the test period. For each test date t 0 t_0 t 0 , the local initial states ( V t 0 , θ t 0 ) (V_{t_0}, \theta_{t_0}) ( V t 0 , θ t 0 ) are inferred from the contemporaneous VIX term structure with ( α , κ ) (\alpha, \kappa) ( α , κ ) fixed at their training-sample values.
Parameter Estimates
The first-step inputs comprise the common term-structure parameters ( α , κ ) (\alpha, \kappa) ( α , κ ) and the date-specific nonnegative states ( V t 0 , θ t 0 ) (V_{t_0}, \theta_{t_0}) ( V t 0 , θ t 0 ) . RISVIX obtains these inputs through the alternating updates described above. Conditional on these inputs, the option-pricing formula uses the four distributional parameters ( σ V , σ θ , λ y , μ y ) (\sigma_V, \sigma_\theta, \lambda_y, \mu_y) ( σ V , σ θ , λ y , μ y ) . Table 4 reports the global parameter values used to generate the option prices in this subsection. RISVIX uses α ^ = 0.8937 \hat{\alpha}=0.8937 α ^ = 0.8937 (H ^ = 0.3937 \hat{H}=0.3937 H ^ = 0.3937 ) and κ ^ = 7.0357 \hat{\kappa}=7.0357 κ ^ = 7.0357 , while LZZ imposes α = 1 \alpha=1 α = 1 and uses κ ^ = 6.2450 \hat{\kappa}=6.2450 κ ^ = 6.2450 . Relative to LZZ, RISVIX has higher σ ^ V \hat{\sigma}_V σ ^ V and σ ^ θ \hat{\sigma}_\theta σ ^ θ , a higher jump intensity λ ^ y \hat{\lambda}^y λ ^ y , and a smaller mean jump size μ ^ y \hat{\mu}_y μ ^ y .
Table 4 . Option-pricing parameters. The table reports the global parameter values used for the option-pricing results. The term-structure inputs (α \alpha α , H H H , κ \kappa κ ) are based on the January 2018–December 2023 term structure, with H = α − 1 / 2 H=\alpha-1/2 H = α − 1/2 . All global parameters are held fixed over the test period.
Table 4.
Pricing Performance
Table 5 reports mean squared error (MSE) and mean absolute error (MAE) for the complete option sample and for maturity and moneyness groups. Moneyness is measured as K / V I X t K/\mathrm{VIX}_t K / VIX t and grouped as < 0.9 <0.9 < 0.9 , ∈ [ 0.9 , 1.1 ] \in[0.9,1.1] ∈ [ 0.9 , 1.1 ] , and > 1.1 >1.1 > 1.1 . The two models are evaluated on the same quotes in every comparison.
Table 5. VIX option-pricing errors. Panel A uses 47,451 training-period quotes; Panel B uses 14,543 test-period quotes. MSE and MAE are computed from option prices in VIX index points. Red. is the percentage error reduction relative to LZZ; positive values favour RISVIX. The lower error in each model pair is in bold.
Table 5.
In sample, RISVIX lowers overall MSE from 0.5548 to 0.5015 and MAE from 0.5815 to 0.5528, corresponding to reductions of 9.61% and 4.93%. Both error measures are lower in every reported group, although the reductions are small for moneyness below 0.9. Out of sample, the overall MSE falls from 0.7238 to 0.6293 and MAE from 0.7537 to 0.7061, reductions of 13.06% and 6.33%. The shorter-maturity group has larger test-period gains, with MSE and MAE reductions of 15.10% and 6.83%, compared with 11.69% and 6.03% at 60–126 days. Across moneyness, RISVIX improves both measures in all three groups; the largest test-period reductions occur below 0.9, at 23.71% for MSE and 12.16% for MAE. The gains are smaller above 1.1, where MSE and MAE decline by 3.93% and 0.48%.
Figure 4 provides selected price cross sections for three training-period dates and one test-period date. Both models reproduce the decline in call prices across strikes. In these examples, LZZ tends to lie above market prices below spot VIX and below market prices at higher strikes. RISVIX reduces this tilt and has a lower mean squared error in each displayed cross section. The fitted curves remain close near the centre of the strike range, with more visible differences toward its edges. These selected examples illustrate cross-sectional shape; the complete-sample comparisons in Table 5 provide the aggregate model ranking.
Figure 4. Cross-sectional VIX option fits. Each panel compares market midquotes with RISVIX and LZZ prices across strikes for a fixed date and expiration. The vertical dotted line marks observed spot VIX. Panels (a)–(c) are in-sample; Panel (d) is out-of-sample. Model prices use the parameter values in Table 4 .
Pricing Performance across VIX States
Given the state dependence in term-structure fit, we next examine whether a similar pattern appears in VIX option prices. Figure 5 reports quote-weighted MSE reductions by the observed VIX level. Within each bin, MSE is computed by pooling all quote errors before taking the ratio between models.
Figure 5. Option-pricing fit across VIX states. Bars report quote-weighted percentage MSE reductions for RISVIX relative to LZZ across bins of observed VIX. Panel (a) is in-sample and Panel (b) is out-of-sample. Bins are left-closed and right-open, and the two panels share a vertical scale; positive values favour RISVIX.
In sample, the reductions are 4.4% below VIX 20 and 4.9% between 20 and 25, rising to 14.4% between 25 and 30, 16.1% between 30 and 40, and 17.9% at VIX 40 and above. Out of sample, RISVIX also improves every bin: MSE reductions are 10.0%, 21.5%, 37.4%, 45.9% and 39.6%, respectively. The largest test-period gain occurs in the 30–40 bin.
The broader split at VIX 25 supports the same pattern. For VIXt _t t ≥ \ge ≥ 25, RISVIX reduces MSE by 15.67% in sample and 40.93% out of sample, compared with overall reductions of 9.61% and 13.06%. The larger gains in elevated-volatility states are consistent with the stronger term-structure improvements during shock dates and the COVID-19 period.
Implied-Volatility Fit
We complement the price-based evaluation with an implied-volatility representation for March 17, 2023. Figure 6 compares market implied volatilities with those obtained from the RISVIX prices for this date. For the Black conversion of model prices, each expiration uses the median forward inferred from its market midquotes and reported implied volatilities, with zero discounting; these conversion inputs are not observed futures quotes.
Figure 6. VIX option implied-volatility surface. Panel (a) compares market and RISVIX implied volatilities across moneyness for four maturities on March 17, 2023. Panel (b) shows the corresponding market points and the interpolated RISVIX surface.
Panel (a) reports slices at 33, 61, 96 and 124 days to expiration. Market and RISVIX implied volatilities both increase with K / V I X t K/\mathrm{VIX}_t K / VIX t over the displayed ranges. The maturity dependence is visible in the different slopes and levels of the curves. RISVIX tracks the upward slope across all four maturities, although its implied volatilities tend to lie below the market values at higher moneyness. The plot uses 62 quotes for which the market and both models have valid implied volatilities under the conversion convention.
Panel (b) displays the same observations as market points and a piecewise-linear RISVIX surface across moneyness and maturity. The surface represents an interpolation of the model values at the observed contracts. The date is selected to illustrate four well-populated, upward-sloping maturity slices.
Comparison with Alternative VIX Option Pricing Models
We compare the following five models:
RISVIX: The proposed rough ISVIX model, which introduces a fractional Volterra kernel into the dynamics of instantaneous squared VIX.
LZZ: The Markovian ISVIX model of Luo et al. [23 ] , with α = 1 \alpha= 1 α = 1 .
HJT: The modulated Volterra stochastic volatility model of Horvath et al. [16 ] , which introduces an additional stochastic factor to capture the VIX implied-volatility smile.
AIL: The Gaussian polynomial volatility models of Abi Jaber et al. [2 ] , in which VIX is derived from conditional integrated variance.
BPS: The rough Hawkes–Heston model of Bondi et al. [7 ] , which combines fractional volatility dynamics with self-exciting jumps and uses affine transform methods for VIX option pricing.
RISVIX and LZZ model instantaneous squared VIX directly, whereas HJT, AIL and BPS derive VIX from the volatility dynamics of the underlying equity index. Table 6 reports their pricing errors on the same VIX call-option sample.
Table 6. Comparison with alternative VIX option pricing models. Calibration: January 2018–December 2023; evaluation: January 2024–June 2025.
RISVIX achieves the lowest MSE and MAE among the five models in both periods. The relatively high calibration MSE of BPS mainly reflects larger pricing errors during the COVID-19 period in 2020. In the evaluation period, BPS is the closest competitor: its MSE is 0.6757, compared with 0.6293 for RISVIX, while their MAEs are nearly identical at 0.7067 and 0.7061, respectively.
Conclusion
This paper develops RISVIX, a fractional extension of the ISVIX framework of Luo et al. [23 ] for modeling the VIX term structure and pricing VIX options. We introduce a fractional Volterra kernel directly into instantaneous squared VIX and derive an explicit term structure governed by Mittag–Leffler decay, together with a semi-closed-form option-pricing formula based on an affine Fourier–Laplace transform. The separation between term-structure and distributional parameters supports a two-step procedure that estimates memory, mean reversion and date-specific initial states from VIX indices before calibrating the remaining parameters to VIX option prices.
Using daily VIX term-structure indices and option quotes from January 2018 to June 2025, we find that RISVIX improves term-structure fit, particularly during periods of market stress. Relative to LZZ, the full-sample MSE reductions are 10.42%, 17.21% and 13.18% for VIX levels, slopes and curvature, respectively. Over January–June 2020, these reductions increase to 17.60%, 41.14% and 46.59%, with an estimated fractional parameter of α ^ = 0.7051 \hat{\alpha} = 0.7051 α ^ = 0.7051 .
The option-pricing results show a similar pattern. With global parameters fixed at their training-sample estimates and date-specific initial states updated from the contemporaneous VIX term structure, RISVIX reduces out-of-sample MSE and MAE by 13.06% and 6.33%, respectively, relative to LZZ. The gains are larger for shorter-maturity options and high-VIX states, with an out-of-sample MSE reduction of 40.93% when V I X t ≥ 25 \mathrm{VIX}_t \ge25 VIX t ≥ 25 . RISVIX also achieves lower aggregate errors than the additional term-structure and option-pricing benchmarks.
Future work could extend RISVIX to a joint SPX–VIX framework to impose consistency across equity and volatility markets ([9 , 27 ] ), or enrich the direct VIX specification with multiple memory factors and more general path-dependent feedback ([15 , 19 ] ).
Acknowledgements
Wen Su acknowledges the support from the EPSRC Centre for Doctoral Training in Mathematics of Random Systems: Analysis, Modelling and Simulation (EP/S023925/1).
Disclosure statement
No potential conflict of interest was reported by the author(s).
Technical Lemmas
Resolvent Identities
For scalar functions f f f and g g g , write ( f ∗ g ) ( t ) = ∫ 0 t f ( t − s ) g ( s ) d s (f \ast g)(t) = \int_{0}^{t} f(t-s)g(s)\,\mathrm{d}s ( f ∗ g ) ( t ) = ∫ 0 t f ( t − s ) g ( s ) d s whenever the integral exists. For κ > 0 \kappa> 0 κ > 0 , the resolvent R κ R_{\kappa} R κ of κ K \kappa K κ K satisfies
R κ + κ K ∗ R κ = κ K , K − R κ ∗ K = 1 κ R κ . (20) R_{\kappa}+\kappa K\ast R_{\kappa}=\kappa K,\qquad K-R_{\kappa}\ast K=\frac{1}{\kappa}R_{\kappa}. \tag*{(20)} R κ + κ K ∗ R κ = κ K , K − R κ ∗ K = κ 1 R κ . (20)
Consequently, for locally integrable f f f and g g g , the equation g = f − κ K ∗ g g=f-\kappa K\ast g g = f − κ K ∗ g is equivalent to g = f − R κ ∗ f g=f-R_{\kappa}\ast f g = f − R κ ∗ f . Both implications follow by substitution and associativity of convolution. This is the scalar resolvent identity underlying the variation-of-constants formula in Abi Jaber et al. [3 ] .
For a locally square-integrable martingale M M M , write ( K ∗ d M ) t = ∫ 0 t K ( t − s ) d M s (K\ast\mathrm{d}M)_t=\int_{0}^{t}K(t-s)\,\mathrm{d}M_s ( K ∗ d M ) t = ∫ 0 t K ( t − s ) d M s . When the convolutions are well defined and stochastic Fubini applies, the same calculation gives
g = f − κ K ∗ g + K ∗ d M ⟺ g = f − R κ ∗ f + 1 κ R κ ∗ d M . (21) g=f-\kappa K\ast g+K\ast\mathrm{d}M\quad\Longleftrightarrow\quad g=f-R_{\kappa}\ast f+\frac{1}{\kappa}R_{\kappa}\ast\mathrm{d}M. \tag*{(21)} g = f − κ K ∗ g + K ∗ d M ⟺ g = f − R κ ∗ f + κ 1 R κ ∗ d M . (21)
Indeed, stochastic Fubini gives R κ ∗ ( K ∗ d M ) = ( R κ ∗ K ) ∗ d M R_{\kappa}\ast(K\ast\mathrm{d}M)=(R_{\kappa}\ast K)\ast\mathrm{d}M R κ ∗ ( K ∗ d M ) = ( R κ ∗ K ) ∗ d M , after which Equation A.1 applies. This argument also permits jumps in M M M ; the identities are understood at times for which the indicated convolutions are defined.
Lemma A.1. Let K ( t ) = t α − 1 / Γ ( α ) K(t)=t^{\alpha-1}/\Gamma(\alpha) K ( t ) = t α − 1 /Γ ( α ) with α ∈ ( 1 / 2 , 1 ] \alpha\in(1/2,1] α ∈ ( 1/2 , 1 ] and κ > 0 \kappa>0 κ > 0 . Then, for t > 0 t>0 t > 0 ,
R κ ( t ) = κ t α − 1 E α , α ( − κ t α ) = − d d t E α , 1 ( − κ t α ) , R_{\kappa}(t)=\kappa t^{\alpha-1}E_{\alpha,\alpha}(-\kappa t^{\alpha})=-\frac{\mathrm{d}}{\mathrm{d}t}E_{\alpha,1}(-\kappa t^{\alpha}), R κ ( t ) = κ t α − 1 E α , α ( − κ t α ) = − d t d E α , 1 ( − κ t α ) , ∫ 0 t R κ ( s ) d s = 1 − E α , 1 ( − κ t α ) , (22) \int_{0}^{t}R_{\kappa}(s)\,\mathrm{d}s=1-E_{\alpha,1}(-\kappa t^{\alpha}), \tag*{(22)} ∫ 0 t R κ ( s ) d s = 1 − E α , 1 ( − κ t α ) , (22) ∫ 0 t E α , 1 ( − κ s α ) d s = t E α , 2 ( − κ t α ) = t [ 1 − κ t α E α , α + 2 ( − κ t α ) ] . \int_{0}^{t}E_{\alpha,1}(-\kappa s^{\alpha})\,\mathrm{d}s=tE_{\alpha,2}(-\kappa t^{\alpha})=t[1-\kappa t^{\alpha}E_{\alpha,\alpha+2}(-\kappa t^{\alpha})]. ∫ 0 t E α , 1 ( − κ s α ) d s = t E α , 2 ( − κ t α ) = t [ 1 − κ t α E α , α + 2 ( − κ t α )] . And the weight w w w in Equation (11) is strictly decreasing on ( 0 , ∞ ) (0,\infty) ( 0 , ∞ ) , with w ( 0 + ) = 1 w(0+)=1 w ( 0 + ) = 1 and w ( ∞ ) = 0 w(\infty)=0 w ( ∞ ) = 0 .
Proof. The beta integral gives K ∗ n ( t ) = t n α − 1 / Γ ( n α ) K^{\ast n}(t)=t^{n\alpha-1}/\Gamma(n\alpha) K ∗ n ( t ) = t n α − 1 /Γ ( n α ) , where K ∗ n K^{\ast n} K ∗ n denotes the n n n -fold convolution of K K K . Hence the series
R κ ( t ) = ∑ n = 1 ∞ ( − 1 ) n − 1 κ n K ∗ n ( t ) = κ t α − 1 E α , α ( − κ t α ) (23) R_{\kappa}(t)=\sum_{n=1}^{\infty}(-1)^{n-1}\kappa^{n}K^{\ast n}(t)=\kappa t^{\alpha-1}E_{\alpha,\alpha}(-\kappa t^{\alpha}) \tag*{(23)} R κ ( t ) = n = 1 ∑ ∞ ( − 1 ) n − 1 κ n K ∗ n ( t ) = κ t α − 1 E α , α ( − κ t α ) (23)
converges in L 1 ( 0 , T ) L^{1}(0,T) L 1 ( 0 , T ) for every finite T T T and satisfies Equation A.1 by cancellation of successive convolution powers. Termwise differentiation of the Mittag–Leffler series on compact subsets of ( 0 , ∞ ) (0,\infty) ( 0 , ∞ ) gives the derivative identity. Integrating this identity and using E α , 1 ( 0 ) = 1 E_{\alpha,1}(0)=1 E α , 1 ( 0 ) = 1 gives the second line of Equation A.3. Termwise integration of the series for E α , 1 E_{\alpha,1} E α , 1 gives the third line, with its second equality following from E α , 2 ( z ) = 1 + z E α , α + 2 ( z ) E_{\alpha,2}(z)=1+zE_{\alpha,\alpha+2}(z) E α , 2 ( z ) = 1 + z E α , α + 2 ( z ) .
For 0 < α < 1 0<\alpha<1 0 < α < 1 , the function t ↦ E α , 1 ( − κ t α ) t\mapsto E_{\alpha,1}(-\kappa t^{\alpha}) t ↦ E α , 1 ( − κ t α ) is strictly decreasing from one to zero; see Mainardi [24 ] . For α = 1 \alpha=1 α = 1 , it equals e − κ t e^{-\kappa t} e − κ t . Since w ( t ) = t − 1 ∫ 0 t E α , 1 ( − κ s α ) d s w(t)=t^{-1}\int_{0}^{t}E_{\alpha,1}(-\kappa s^{\alpha})\,\mathrm{d}s w ( t ) = t − 1 ∫ 0 t E α , 1 ( − κ s α ) d s , its limits follow by continuity at zero and averaging at infinity. Moreover, w ′ ( t ) = [ E α , 1 ( − κ t α ) − w ( t ) ] / t < 0 w'(t) = [E_{\alpha,1}(-\kappa t^\alpha) - w(t)]/t < 0 w ′ ( t ) = [ E α , 1 ( − κ t α ) − w ( t )] / t < 0 , which proves strict monotonicity.
Remark. Integrating the first identity in Equation A.3 over [ s , s + τ 0 ] [s, s+\tau_0] [ s , s + τ 0 ] gives the averaged resolvent in Equation (13) . This expression extends continuously to s = 0 s=0 s = 0 ; for α = 1 \alpha=1 α = 1 , it reduces to a τ 0 ( s ) = w ( τ 0 ) e − κ s a_{\tau_0}(s) = w(\tau_0)e^{-\kappa s} a τ 0 ( s ) = w ( τ 0 ) e − κ s .
Moment Conditions
Fix T > 0 T>0 T > 0 , α ∈ ( 1 / 2 , 1 ) \alpha\in(1/2,1) α ∈ ( 1/2 , 1 ) , κ , μ y > 0 \kappa,\mu_y>0 κ , μ y > 0 and σ V , σ θ , λ y ≥ 0 \sigma_V,\sigma_\theta,\lambda^y\geq0 σ V , σ θ , λ y ≥ 0 . Let ( V , θ ) (V,\theta) ( V , θ ) be a nonnegative weak solution of Equation (9) on [ 0 , T ] [0,T] [ 0 , T ] , with the driving processes, jump law and independence conditions specified immediately after that equation. The initial states are nonnegative and F 0 \mathcal{F}_0 F 0 -measurable, with E Q [ V 0 2 + θ 0 2 ] < ∞ \mathbb{E}^{\mathbb{Q}}[V_0^2+\theta_0^2]<\infty E Q [ V 0 2 + θ 0 2 ] < ∞ . We require θ \theta θ to be continuous and V V V to be predictable with ∫ 0 T V s d s < ∞ \int_0^T V_s\,\mathrm{d}s<\infty ∫ 0 T V s d s < ∞ almost surely; nonnegativity of V V V and its Volterra equation are understood Q ⊗ d t \mathbb{Q}\otimes\mathrm{d}t Q ⊗ d t -almost everywhere.
Lemma A.2. Under the preceding assumptions, θ \theta θ is a square-integrable martingale and, for 0 ≤ s ≤ t ≤ T 0\leq s\leq t\leq T 0 ≤ s ≤ t ≤ T , E s Q [ θ t ] = θ s \mathbb{E}^{\mathbb{Q}}_s[\theta_t]=\theta_s E s Q [ θ t ] = θ s and E s Q [ θ t 2 ] = θ s 2 + σ θ 2 ( t − s ) θ s \mathbb{E}^{\mathbb{Q}}_s[\theta_t^2]=\theta_s^2+\sigma_\theta^2(t-s)\theta_s E s Q [ θ t 2 ] = θ s 2 + σ θ 2 ( t − s ) θ s almost surely. Moreover,
E Q [ ∫ 0 T ( V t 2 + θ t 2 ) d t ] ≤ C T ( 1 + E Q [ V 0 2 + θ 0 2 ] ) < ∞ , (24) \mathbb{E}^{\mathbb{Q}}\left[\int_0^T (V_t^2+\theta_t^2)\,\mathrm{d}t\right]\leq C_T\left(1+\mathbb{E}^{\mathbb{Q}}[V_0^2+\theta_0^2]\right)<\infty, \tag*{(24)} E Q [ ∫ 0 T ( V t 2 + θ t 2 ) d t ] ≤ C T ( 1 + E Q [ V 0 2 + θ 0 2 ] ) < ∞ , (24) where C T C_T C T depends only on T , α , κ , σ V , σ θ , λ y , μ y T,\alpha,\kappa,\sigma_V,\sigma_\theta,\lambda^y,\mu_y T , α , κ , σ V , σ θ , λ y , μ y . In particular, E 0 Q [ ∫ 0 T ( V t 2 + θ t 2 ) d t ] < ∞ \mathbb{E}^{\mathbb{Q}}_0\left[\int_0^T(V_t^2+\theta_t^2)\,\mathrm{d}t\right]<\infty E 0 Q [ ∫ 0 T ( V t 2 + θ t 2 ) d t ] < ∞ almost surely.
Proof. Since θ \theta θ is a nonnegative local martingale, it is a supermartingale and E Q [ ∫ 0 T θ s d s ] ≤ T E Q [ θ 0 ] < ∞ \mathbb{E}^{\mathbb{Q}}[\int_0^T\theta_s\,\mathrm{d}s]\leq T\mathbb{E}^{\mathbb{Q}}[\theta_0]<\infty E Q [ ∫ 0 T θ s d s ] ≤ T E Q [ θ 0 ] < ∞ . Thus its Brownian integral is square integrable, and so is θ \theta θ because θ 0 ∈ L 2 ( Q ) \theta_0\in L^2(\mathbb{Q}) θ 0 ∈ L 2 ( Q ) . Conditional Itô isometry and the martingale property give E s Q [ ( θ t − θ s ) 2 ] = σ θ 2 ( t − s ) θ s \mathbb{E}^{\mathbb{Q}}_s[(\theta_t-\theta_s)^2]=\sigma_\theta^2(t-s)\theta_s E s Q [( θ t − θ s ) 2 ] = σ θ 2 ( t − s ) θ s , proving the conditional moment identities. Consequently,
E Q [ ∫ 0 T θ t 2 d t ] = T E Q [ θ 0 2 ] + σ θ 2 T 2 2 E Q [ θ 0 ] . (25) \mathbb{E}^{\mathbb{Q}}\left[\int_0^T\theta_t^2\,\mathrm{d}t\right]
=T\mathbb{E}^{\mathbb{Q}}[\theta_0^2]+\frac{\sigma_\theta^2T^2}{2}\mathbb{E}^{\mathbb{Q}}[\theta_0]. \tag*{(25)} E Q [ ∫ 0 T θ t 2 d t ] = T E Q [ θ 0 2 ] + 2 σ θ 2 T 2 E Q [ θ 0 ] . (25)
Let R κ R_\kappa R κ be the resolvent in Equation A.1 and set q κ ( t ) = E α , 1 ( − κ t α ) q_\kappa(t)=E_{\alpha,1}(-\kappa t^\alpha) q κ ( t ) = E α , 1 ( − κ t α ) . By Lemma A.1 , q κ q_\kappa q κ is bounded, R κ ∈ L 1 ( 0 , T ) ∩ L 2 ( 0 , T ) R_\kappa\in L^1(0,T)\cap L^2(0,T) R κ ∈ L 1 ( 0 , T ) ∩ L 2 ( 0 , T ) and K ∈ L 2 ( 0 , T ) K\in L^2(0,T) K ∈ L 2 ( 0 , T ) . Define
τ n = inf { t ∈ [ 0 , T ] : ∫ 0 t V s d s ≥ n } ∧ T , I s ( n ) = 1 { s ≤ τ n } . (26) \tau_n=\inf\left\{t\in[0,T]:\int_0^t V_s\,\mathrm{d}s\geq n\right\}\wedge T,\qquad I_s^{(n)}=\mathbf{1}_{\{s\leq\tau_n\}}. \tag*{(26)} τ n = inf { t ∈ [ 0 , T ] : ∫ 0 t V s d s ≥ n } ∧ T , I s ( n ) = 1 { s ≤ τ n } . (26)
The cumulative integral is continuous and adapted, so τ n \tau_n τ n is a stopping time and I ( n ) I^{(n)} I ( n ) is predictable. Pathwise integrability gives ∫ 0 T I s ( n ) V s d s ≤ n \int_0^T I_s^{(n)}V_s\,\mathrm{d}s\leq n ∫ 0 T I s ( n ) V s d s ≤ n and τ n ↑ T \tau_n\uparrow T τ n ↑ T almost surely. For any deterministic h ∈ L 2 ( 0 , T ) h \in L^2(0,T) h ∈ L 2 ( 0 , T ) , Itô isometry and Tonelli’s theorem yield
E Q [ ∫ 0 T ∣ ∫ 0 t h ( t − s ) I s ( n ) V s d B 1 , s ∣ 2 d t ] = E Q [ ∫ 0 T ∫ 0 t h ( t − s ) 2 I s ( n ) V s d s d t ] ≤ n ∥ h ∥ L 2 ( 0 , T ) 2 . (27) \mathbb{E}^{\mathbb{Q}}\left[\int_0^T\left|\int_0^t h(t-s)I_s^{(n)}\sqrt{V_s}\,\mathrm{d}B_{1,s}\right|^2\mathrm{d}t\right]
=
\mathbb{E}^{\mathbb{Q}}\left[\int_0^T\int_0^t h(t-s)^2I_s^{(n)}V_s\,\mathrm{d}s\,\mathrm{d}t\right]
\le n\|h\|_{L^2(0,T)}^2. \tag*{(27)} E Q [ ∫ 0 T ∫ 0 t h ( t − s ) I s ( n ) V s d B 1 , s 2 d t ] = E Q [ ∫ 0 T ∫ 0 t h ( t − s ) 2 I s ( n ) V s d s d t ] ≤ n ∥ h ∥ L 2 ( 0 , T ) 2 . (27)
In particular, r ↦ ( ∫ 0 T K ( r − s ) 2 I s ( n ) V s d s ) 1 / 2 r \mapsto\left(\int_0^T K(r-s)^2I_s^{(n)}V_s\,\mathrm{d}s\right)^{1/2} r ↦ ( ∫ 0 T K ( r − s ) 2 I s ( n ) V s d s ) 1/2 belongs to L 2 ( 0 , T ) L^2(0,T) L 2 ( 0 , T ) almost surely. Its convolution with ∣ R κ ∣ |\mathcal{R}_{\kappa}| ∣ R κ ∣ is therefore finite almost everywhere, verifying the stochastic-Fubini condition for the stopped Brownian term in Equation A.2. Ordinary Fubini applies pathwise to the drift and finite-activity jump terms. On t < τ n t<\tau_n t < τ n , the stopped Brownian integrals and their past time convolutions coincide with the unstopped ones. Thus Equation A.2 and Equation A.1 show that the process
V t ( n ) = q κ ( t ) V 0 + ∫ 0 t R κ ( t − s ) θ s d s + σ V κ ∫ 0 t R κ ( t − s ) I s ( n ) V s d B 1 , s + 1 κ ∫ 0 t R κ ( t − s ) d J s Y (28) V_t^{(n)}=q_{\kappa}(t)V_0+\int_0^t R_{\kappa}(t-s)\theta_s\,\mathrm{d}s
+\frac{\sigma_V}{\kappa}\int_0^t R_{\kappa}(t-s)I_s^{(n)}\sqrt{V_s}\,\mathrm{d}B_{1,s}
+\frac{1}{\kappa}\int_0^t R_{\kappa}(t-s)\,\mathrm{d}J_s^Y \tag*{(28)} V t ( n ) = q κ ( t ) V 0 + ∫ 0 t R κ ( t − s ) θ s d s + κ σ V ∫ 0 t R κ ( t − s ) I s ( n ) V s d B 1 , s + κ 1 ∫ 0 t R κ ( t − s ) d J s Y (28)
agrees with V t V_t V t on { t < τ n } \{t<\tau_n\} { t < τ n } , Q ⊗ d t \mathbb{Q}\otimes\mathrm{d}t Q ⊗ d t -almost everywhere. Set F n = E Q [ ∫ 0 T ∣ V t ( n ) ∣ 2 d t ] F_n=\mathbb{E}^{\mathbb{Q}}\left[\int_0^T|V_t^{(n)}|^2\mathrm{d}t\right] F n = E Q [ ∫ 0 T ∣ V t ( n ) ∣ 2 d t ] . Young’s convolution inequality, Equation A.8 with h = R κ h=R_{\kappa} h = R κ , and the compensated-Poisson isometry with E Q [ Y 1 2 ] = 2 μ y 2 \mathbb{E}^{\mathbb{Q}}[Y_1^2]=2\mu_y^2 E Q [ Y 1 2 ] = 2 μ y 2 give
F n ≤ 4 T ∥ q κ ∥ L ∞ ( 0 , T ) 2 E Q [ V 0 2 ] + 4 ∥ R κ ∥ L 1 ( 0 , T ) 2 E Q [ ∫ 0 T θ s 2 d s ] + 8 λ μ y 2 T κ 2 ∥ R κ ∥ L 2 ( 0 , T ) 2 + 4 σ V 2 κ 2 ∥ R κ ∥ L 2 ( 0 , T ) 2 E Q [ ∫ 0 T I s ( n ) V s d s ] . (29) \begin{aligned}
F_n \le{}&4T\|q_{\kappa}\|_{L^\infty(0,T)}^2\mathbb{E}^{\mathbb{Q}}[V_0^2]
+4\|R_{\kappa}\|_{L^1(0,T)}^2\mathbb{E}^{\mathbb{Q}}\left[\int_0^T\theta_s^2\,\mathrm{d}s\right]\\
&+\frac{8\lambda\mu_y^2T}{\kappa^2}\|R_{\kappa}\|_{L^2(0,T)}^2
+\frac{4\sigma_V^2}{\kappa^2}\|R_{\kappa}\|_{L^2(0,T)}^2\mathbb{E}^{\mathbb{Q}}\left[\int_0^T I_s^{(n)}V_s\,\mathrm{d}s\right]. \tag*{(29)}
\end{aligned} F n ≤ 4 T ∥ q κ ∥ L ∞ ( 0 , T ) 2 E Q [ V 0 2 ] + 4∥ R κ ∥ L 1 ( 0 , T ) 2 E Q [ ∫ 0 T θ s 2 d s ] + κ 2 8 λ μ y 2 T ∥ R κ ∥ L 2 ( 0 , T ) 2 + κ 2 4 σ V 2 ∥ R κ ∥ L 2 ( 0 , T ) 2 E Q [ ∫ 0 T I s ( n ) V s d s ] . (29)
The right-hand side is finite for each n n n , since its last expectation is at most n n n ; hence F n < ∞ F_n<\infty F n < ∞ before any absorption argument. Agreement before τ n \tau_n τ n , nonnegativity of V V V , and Cauchy–Schwarz imply
E Q [ ∫ 0 T I s ( n ) V s d s ] ≤ E Q [ ∫ 0 T V s ( n ) d s ] ≤ T F n . (30) \mathbb{E}^{\mathbb{Q}}\left[\int_0^T I_s^{(n)}V_s\,\mathrm{d}s\right]
\le\mathbb{E}^{\mathbb{Q}}\left[\int_0^T V_s^{(n)}\,\mathrm{d}s\right]
\le\sqrt{TF_n}. \tag*{(30)} E Q [ ∫ 0 T I s ( n ) V s d s ] ≤ E Q [ ∫ 0 T V s ( n ) d s ] ≤ T F n . (30)
Substituting this bound into Equation A.10, using the moment formula for θ \theta θ and E Q [ θ 0 ] ≤ ( 1 + E Q [ θ 0 2 ] ) / 2 \mathbb{E}^{\mathbb{Q}}[\theta_0]\le(1+\mathbb{E}^{\mathbb{Q}}[\theta_0^2])/2 E Q [ θ 0 ] ≤ ( 1 + E Q [ θ 0 2 ]) /2 , and applying Young’s inequality give F n ≤ C T ( 1 + E Q [ V 0 2 + θ 0 2 ] ) + F n / 2 F_n\le C_T(1+\mathbb{E}^{\mathbb{Q}}[V_0^2+\theta_0^2])+F_n/2 F n ≤ C T ( 1 + E Q [ V 0 2 + θ 0 2 ]) + F n /2 , where C T C_T C T depends only on the stated quantities and not on n n n . Since τ n ↑ T \tau_n\uparrow T τ n ↑ T , the agreement before τ n \tau_n τ n gives V ( n ) → V V^{(n)}\to V V ( n ) → V almost everywhere under Q ⊗ d t \mathbb{Q}\otimes\mathrm{d}t Q ⊗ d t . Fatou’s lemma and the bound for θ \theta θ prove Equation A.5. Conditioning the resulting integrable random variable on F 0 \mathcal{F}_0 F 0 proves the final assertion.
Lemma A.3. Suppose that Equation (9) admits the solution specified above on [ 0 , T ] [0,T] [ 0 , T ] . Then there exists c T > 0 c_T>0 c T > 0 , depending only on T T T and the model parameters, such that
E 0 Q [ exp { c ∫ 0 T V t d t } ] < ∞ (31) \mathbb{E}^{\mathbb{Q}}_0\left[\exp\left\{c\int_0^T V_t\,\mathrm{d}t\right\}\right]<\infty \tag*{(31)} E 0 Q [ exp { c ∫ 0 T V t d t } ] < ∞ (31)
almost surely for every c ∈ [ 0 , c T ) c \in[0,c_T) c ∈ [ 0 , c T ) .
Proof. Let X T = ∫ 0 T V t d t X_T=\int_0^T V_t\,\mathrm{d}t X T = ∫ 0 T V t d t and q κ ( t ) = E α , 1 ( − κ t α ) q_\kappa(t)=E_{\alpha,1}(-\kappa t^\alpha) q κ ( t ) = E α , 1 ( − κ t α ) . The resolvent calculation in the proof of Lemma A.2 , followed by integration in time, gives
X T = V 0 ∫ 0 T q κ ( t ) d t + ∫ 0 T b T ( s ) θ s d s + σ V κ ∫ 0 T b T ( s ) V s d B 1 , s + 1 κ ∫ 0 T b T ( s ) d J s Y . (32) X_T=V_0\int_0^T q_\kappa(t)\,\mathrm{d}t+\int_0^T b_T(s)\theta_s\,\mathrm{d}s
+\frac{\sigma_V}{\kappa}\int_0^T b_T(s)\sqrt{V_s}\,\mathrm{d}B_{1,s}+\frac{1}{\kappa}\int_0^T b_T(s)\,\mathrm{d}J_s^Y. \tag*{(32)} X T = V 0 ∫ 0 T q κ ( t ) d t + ∫ 0 T b T ( s ) θ s d s + κ σ V ∫ 0 T b T ( s ) V s d B 1 , s + κ 1 ∫ 0 T b T ( s ) d J s Y . (32)
where b T ( s ) = 1 − q κ ( T − s ) ∈ [ 0 , 1 ] b_T(s)=1-q_\kappa(T-s)\in[0,1] b T ( s ) = 1 − q κ ( T − s ) ∈ [ 0 , 1 ] by Lemma A.1 . Here R κ ∈ L 2 ( 0 , T ) R_\kappa\in L^2(0,T) R κ ∈ L 2 ( 0 , T ) and E Q [ X T ] < ∞ \mathbb{E}^{\mathbb{Q}}[X_T]<\infty E Q [ X T ] < ∞ justify stochastic Fubini for the Brownian term, while ordinary Fubini applies to the drift and finite-activity jump terms. Throughout the proof, C T C_T C T denotes a positive deterministic constant depending only on T T T and the model parameters, which may increase from line to line.
Choose a > 0 a>0 a > 0 such that a σ θ 2 T ≤ 1 a\sigma_\theta^2T\leq1 a σ θ 2 T ≤ 1 . For fixed s ∈ [ 0 , T ] s\in[0,T] s ∈ [ 0 , T ] , applying Itô’s formula to exp { a θ v / [ 1 − a σ θ 2 ( s − v ) / 2 ] } \exp\{a\theta_v/[1-a\sigma_\theta^2(s-v)/2]\} exp { a θ v / [ 1 − a σ θ 2 ( s − v ) /2 ]} , 0 ≤ v ≤ s 0\leq v\leq s 0 ≤ v ≤ s , gives zero drift. Restricting first to { θ 0 ≤ m } \{\theta_0\leq m\} { θ 0 ≤ m } and stopping this exponential at increasing levels of θ \theta θ , conditional Fatou’s lemma, followed by m ↑ ∞ m\uparrow\infty m ↑ ∞ , yields
E 0 Q [ e a θ s ] ≤ exp { a θ 0 1 − a σ θ 2 s / 2 } ≤ e 2 a θ 0 . (33) \mathbb{E}^{\mathbb{Q}}_0[e^{a\theta_s}]
\leq\exp\left\{\frac{a\theta_0}{1-a\sigma_\theta^2s/2}\right\}
\leq e^{2a\theta_0}. \tag*{(33)} E 0 Q [ e a θ s ] ≤ exp { 1 − a σ θ 2 s /2 a θ 0 } ≤ e 2 a θ 0 . (33)
Conditional Markov’s inequality gives Q ( θ s > 2 θ 0 + y ∣ F 0 ) ≤ e − a y \mathbb{Q}(\theta_s>2\theta_0+y\mid\mathcal{F}_0)\leq e^{-ay} Q ( θ s > 2 θ 0 + y ∣ F 0 ) ≤ e − a y for y ≥ 0 y\geq0 y ≥ 0 . Integrating this tail bound shows that ( E 0 Q [ θ s p ] ) 1 / p ≤ 2 θ 0 + ( p ! ) 1 / p / a ≤ 2 θ 0 + p / a (\mathbb{E}^{\mathbb{Q}}_0[\theta_s^p])^{1/p}\leq2\theta_0+(p!)^{1/p}/a\leq2\theta_0+p/a ( E 0 Q [ θ s p ] ) 1/ p ≤ 2 θ 0 + ( p ! ) 1/ p / a ≤ 2 θ 0 + p / a for every integer p ≥ 2 p\geq2 p ≥ 2 . Conditional Minkowski’s inequality then bounds the conditional L p L^p L p norm of ∫ 0 T b T ( s ) θ s d s \int_0^T b_T(s)\theta_s\,\mathrm{d}s ∫ 0 T b T ( s ) θ s d s by C T ( θ 0 + p ) C_T(\theta_0+p) C T ( θ 0 + p ) .
For the jump term, ∣ ∫ 0 T b T ( s ) d J s Y ∣ ≤ ∑ j = 1 N T Y Y j + λ Y μ Y T \left|\int_0^T b_T(s)\,\mathrm{d}J_s^Y\right|\leq\sum_{j=1}^{N_T^Y}Y_j+\lambda^Y\mu_YT ∫ 0 T b T ( s ) d J s Y ≤ ∑ j = 1 N T Y Y j + λ Y μ Y T . Independence from F 0 \mathcal{F}_0 F 0 and the exponential jump law give E 0 Q [ exp { ( 2 μ Y ) − 1 ∑ j = 1 N T Y Y j } ] = e λ Y T \mathbb{E}^{\mathbb{Q}}_0[\exp\{(2\mu_Y)^{-1}\sum_{j=1}^{N_T^Y}Y_j\}]=e^{\lambda^YT} E 0 Q [ exp {( 2 μ Y ) − 1 ∑ j = 1 N T Y Y j }] = e λ Y T . Thus the last term in Equation A.13 has conditional L p L^p L p norm at most C T p C_Tp C T p for every integer p ≥ 2 p\geq2 p ≥ 2 .
Set x p = ( E 0 Q [ X T p ] ) 1 / p x_p=(\mathbb{E}^{\mathbb{Q}}_0[X_T^p])^{1/p} x p = ( E 0 Q [ X T p ] ) 1/ p . By Lemma A.2 and Cauchy–Schwarz, x 2 < ∞ x_2<\infty x 2 < ∞ almost surely. The continuous-martingale Burkholder–Davis–Gundy inequality, with L p L^p L p constant 2 2 p 2\sqrt{2p} 2 2 p , and b T 2 V ≤ V b_T^2V\leq V b T 2 V ≤ V give, for dyadic integers p ≥ 4 p\geq4 p ≥ 4 with x p / 2 < ∞ x_{p/2}<\infty x p /2 < ∞ ,
x p ≤ C T ( 1 + V 0 + θ 0 + p ) + C T p x p / 2 1 / 2 . (34) x_p\leq C_T(1+V_0+\theta_0+p)+C_T\sqrt{p}\,x_{p/2}^{1/2}. \tag*{(34)} x p ≤ C T ( 1 + V 0 + θ 0 + p ) + C T p x p /2 1/2 . (34)
For the conditional inequality, first stop the Brownian integral and restrict to F 0 \mathcal{F}_0 F 0 -measurable events on which x p / 2 x_{p/2} x p /2 is bounded; conditional Fatou removes these localisations. Thus the right-hand side is finite before any absorption argument, and dyadic induction gives x p < ∞ x_p<\infty x p < ∞ for every dyadic p ≥ 2 p\geq2 p ≥ 2 .
For dyadic p ≥ 4 p\geq4 p ≥ 4 , use x p / 2 ≤ x p x_{p/2}\leq x_p x p /2 ≤ x p and Young’s inequality in Equation A.15 to absorb half of x p x_p x p . Monotonicity of conditional L p L^p L p norms between successive dyadic integers then gives x p ≤ A T + C T p x_p\leq A_T+C_Tp x p ≤ A T + C T p for every integer p ≥ 2 p\geq2 p ≥ 2 , where A T = C T ( 1 + V 0 + θ 0 ) < ∞ A_T=C_T(1+V_0+\theta_0)<\infty A T = C T ( 1 + V 0 + θ 0 ) < ∞ almost surely and C T C_T C T is deterministic. Hence E 0 Q [ X T n ] ≤ 2 n − 1 ( A T n + C T n n n ) \mathbb{E}^{\mathbb{Q}}_0[X_T^n]\leq2^{n-1}(A_T^n+C_T^nn^n) E 0 Q [ X T n ] ≤ 2 n − 1 ( A T n + C T n n n ) for every integer n ≥ 2 n\geq2 n ≥ 2 . Conditional monotone convergence and n ! ≥ ( n / e ) n n!\geq(n/e)^n n ! ≥ ( n / e ) n in the exponential series show that E 0 Q [ e c X T ] < ∞ \mathbb{E}^{\mathbb{Q}}_0[e^{cX_T}]<\infty E 0 Q [ e c X T ] < ∞ whenever 2 e C T c < 1 2eC_Tc<1 2 e C T c < 1 . Taking c T = ( 2 e C T ) − 1 c_T=(2eC_T)^{-1} c T = ( 2 e C T ) − 1 proves Equation A.12.
Fix T , τ 0 > 0 T,\tau_0>0 T , τ 0 > 0 . For η ≥ 0 \eta\ge0 η ≥ 0 , consider the real Riccati–Volterra equation
p η ( s ) = η a τ 0 ( s ) + σ V 2 2 κ ∫ 0 s R κ ( s − v ) p η ( v ) 2 d v , 0 ≤ s ≤ T , (35) p_\eta(s)=\eta a_{\tau_0}(s)+\frac{\sigma_V^2}{2\kappa}\int_0^s R_\kappa(s-v)p_\eta(v)^2\,\mathrm{d}v,\qquad0\le s\le T, \tag*{(35)} p η ( s ) = η a τ 0 ( s ) + 2 κ σ V 2 ∫ 0 s R κ ( s − v ) p η ( v ) 2 d v , 0 ≤ s ≤ T , (35)
together with the Riccati equation
ψ 2 ′ ( s ; η ) = κ p η ( s ) + σ θ 2 2 ψ 2 ( s ; η ) 2 , ψ 2 ( 0 ; η ) = η [ 1 − w ( τ 0 ) ] . (36) \psi_2'(s;\eta)=\kappa p_\eta(s)+\frac{\sigma_\theta^2}{2}\psi_2(s;\eta)^2,\qquad\psi_2(0;\eta)=\eta[1-w(\tau_0)]. \tag*{(36)} ψ 2 ′ ( s ; η ) = κ p η ( s ) + 2 σ θ 2 ψ 2 ( s ; η ) 2 , ψ 2 ( 0 ; η ) = η [ 1 − w ( τ 0 )] . (36)
Lemma A.4. There exists η ˉ T , τ 0 > 0 \bar{\eta}_{T,\tau_0}>0 η ˉ T , τ 0 > 0 such that, for every η ∈ [ 0 , η ˉ T , τ 0 ) \eta\in[0,\bar{\eta}_{T,\tau_0}) η ∈ [ 0 , η ˉ T , τ 0 ) , the preceding equations admit a unique nonnegative continuous solution ( p η , ψ 2 ( ⋅ ; η ) ) (p_\eta,\psi_2(\cdot;\eta)) ( p η , ψ 2 ( ⋅ ; η )) on [ 0 , T ] [0,T] [ 0 , T ] , with
sup 0 ≤ s ≤ T p η ( s ) ≤ 2 η ∥ a τ 0 ∥ L ∞ ( 0 , T ) < μ y − 1 , \sup_{0\le s\le T}p_\eta(s)\le2\eta\|a_{\tau_0}\|_{L^\infty(0,T)}<\mu_y^{-1}, 0 ≤ s ≤ T sup p η ( s ) ≤ 2 η ∥ a τ 0 ∥ L ∞ ( 0 , T ) < μ y − 1 , sup 0 ≤ s ≤ T ψ 2 ( s ; η ) ≤ 2 η [ 1 − w ( τ 0 ) + 2 κ T ∥ a τ 0 ∥ L ∞ ( 0 , T ) ] . (37) \sup_{0\le s\le T}\psi_2(s;\eta)\le2\eta\left[1-w(\tau_0)+2\kappa T\|a_{\tau_0}\|_{L^\infty(0,T)}\right]. \tag*{(37)} 0 ≤ s ≤ T sup ψ 2 ( s ; η ) ≤ 2 η [ 1 − w ( τ 0 ) + 2 κ T ∥ a τ 0 ∥ L ∞ ( 0 , T ) ] . (37)
In particular, both functions converge uniformly to zero as η ↓ 0 \eta\downarrow0 η ↓ 0 .
Proof. By Lemma A.1, R κ R_\kappa R κ is nonnegative and belongs to L 1 ( 0 , T ) L^1(0,T) L 1 ( 0 , T ) , a τ 0 a_{\tau_0} a τ 0 is nonnegative and continuous on [ 0 , T ] [0,T] [ 0 , T ] , and 0 < w ( τ 0 ) < 1 0<w(\tau_0)<1 0 < w ( τ 0 ) < 1 . Let A = ∥ a τ 0 ∥ L ∞ ( 0 , T ) A=\|a_{\tau_0}\|_{L^\infty(0,T)} A = ∥ a τ 0 ∥ L ∞ ( 0 , T ) and L = ∥ R κ ∥ L 1 ( 0 , T ) L=\|R_\kappa\|_{L^1(0,T)} L = ∥ R κ ∥ L 1 ( 0 , T ) . For η > 0 \eta>0 η > 0 , define Γ η p = η a τ 0 + ( σ V 2 / ( 2 κ ) ) R κ ∗ p 2 \Gamma_\eta p=\eta a_{\tau_0}+(\sigma_V^2/(2\kappa))R_\kappa*p^2 Γ η p = η a τ 0 + ( σ V 2 / ( 2 κ )) R κ ∗ p 2 . On the closed ball { p ∈ C ( [ 0 , T ] ; R ) : ∥ p ∥ ∞ ≤ 2 η A } \{p\in C([0,T];\mathbb{R}):\|p\|_\infty\le2\eta A\} { p ∈ C ([ 0 , T ] ; R ) : ∥ p ∥ ∞ ≤ 2 η A } ,
∥ Γ η p ∥ ∞ ≤ η A + 2 σ V 2 κ L η 2 A 2 , ∥ Γ η p − Γ η q ∥ ∞ ≤ 2 σ V 2 κ L η A ∥ p − q ∥ ∞ . (38) \|\Gamma_\eta p\|_\infty\le\eta A+\frac{2\sigma_V^2}{\kappa}L\eta^2A^2,\qquad\|\Gamma_\eta p-\Gamma_\eta q\|_\infty\le\frac{2\sigma_V^2}{\kappa}L\eta A\|p-q\|_\infty. \tag*{(38)} ∥ Γ η p ∥ ∞ ≤ η A + κ 2 σ V 2 L η 2 A 2 , ∥ Γ η p − Γ η q ∥ ∞ ≤ κ 2 σ V 2 L η A ∥ p − q ∥ ∞ . (38)
For sufficiently small η \eta η , this map preserves the ball and is a contraction. Iteration from zero gives a nonnegative fixed point p η p_\eta p η . Reducing η ˉ T , τ 0 \bar{\eta}_{T,\tau_0} η ˉ T , τ 0 if necessary ensures 2 η A < μ y − 1 2\eta A<\mu_y^{-1} 2 η A < μ y − 1 . If p ~ η \widetilde p_\eta p η is another nonnegative continuous solution, both are bounded and their difference is bounded in absolute value by ( σ V 2 / ( 2 κ ) ) R κ ∗ [ ( p η + p ~ η ) ∣ p η − p ~ η ∣ ] (\sigma_V^2/(2\kappa))R_\kappa*[(p_\eta+\widetilde p_\eta)|p_\eta-\widetilde p_\eta|] ( σ V 2 / ( 2 κ )) R κ ∗ [( p η + p η ) ∣ p η − p η ∣ ] . Since ∫ 0 δ R κ ( s ) d s → 0 \int_0^\delta R_\kappa(s)\,\mathrm{d}s\to0 ∫ 0 δ R κ ( s ) d s → 0 as δ ↓ 0 \delta\downarrow0 δ ↓ 0 , contraction on successive intervals of a common sufficiently short length gives uniqueness on [ 0 , T ] [0,T] [ 0 , T ] .
Given p η p_\eta p η , the right-hand side of Equation A.17 is continuous in time and locally Lipschitz in ψ 2 \psi_2 ψ 2 , so it has a unique maximal solution. Its derivative is nonnegative, and hence ψ 2 ( s ; η ) ≥ 0 \psi_2(s;\eta)\ge0 ψ 2 ( s ; η ) ≥ 0 on its interval of existence. Put C = 1 − w ( τ 0 ) + 2 κ T A > 0 C=1-w(\tau_0)+2\kappa T A>0 C = 1 − w ( τ 0 ) + 2 κ T A > 0 and reduce η ˉ T , τ 0 \bar{\eta}_{T,\tau_0} η ˉ T , τ 0 further so that 2 σ θ 2 T η C < 1 2\sigma_\theta^2T\eta C<1 2 σ θ 2 T η C < 1 . If ψ 2 \psi_2 ψ 2 first reached 2 η C 2\eta C 2 η C at some s ≤ T s\le T s ≤ T , its integral equation would give ψ 2 ( s ; η ) ≤ η C + 2 σ θ 2 T η 2 C 2 < 2 η C \psi_2(s;\eta)\le\eta C+2\sigma_\theta^2T\eta^2C^2<2\eta C ψ 2 ( s ; η ) ≤ η C + 2 σ θ 2 T η 2 C 2 < 2 η C , a contradiction. Thus the solution and its derivative remain bounded. It has a finite limit at any endpoint up to T T T , and local existence prevents termination before T T T . It therefore extends to [ 0 , T ] [0,T] [ 0 , T ] and satisfies the second bound in Equation A.18. At η = 0 \eta=0 η = 0 , both solutions are identically zero, with uniqueness following from the same arguments. The two bounds prove uniform convergence as η ↓ 0 \eta\downarrow0 η ↓ 0 . □
For u ∈ C u\in\mathbb{C} u ∈ C , consider
ψ 1 ( s ; u ) = u a τ 0 ( s ) + σ V 2 2 κ ∫ 0 s R κ ( s − v ) ψ 1 ( v ; u ) 2 d v . (39) \psi_1(s;u)=u a_{\tau_0}(s)+\frac{\sigma_V^2}{2\kappa}\int_0^s R_\kappa(s-v)\psi_1(v;u)^2\,\mathrm{d}v. \tag*{(39)} ψ 1 ( s ; u ) = u a τ 0 ( s ) + 2 κ σ V 2 ∫ 0 s R κ ( s − v ) ψ 1 ( v ; u ) 2 d v . (39)
Let η ∈ [ 0 , η ˉ T , τ 0 ) \eta\in[0,\bar{\eta}_{T,\tau_0}) η ∈ [ 0 , η ˉ T , τ 0 ) and ξ ∈ R \xi\in\mathbb{R} ξ ∈ R , and let p η p_\eta p η be the solution supplied by Lemma A.4. Suppose that Equation A.20 with u = η + i ξ u=\eta+i\xi u = η + i ξ has a continuous solution ψ 1 ( ⋅ ; η + i ξ ) \psi_1(\mathord{\cdot};\eta+i\xi) ψ 1 ( ⋅ ; η + i ξ ) on [ 0 , T ] [0,T] [ 0 , T ] .
Lemma A.5. Under the preceding assumptions, this solution is unique and satisfies
Re ψ 1 ( s ; η + i ξ ) ≤ p η ( s ) < μ y − 1 , 0 ≤ s ≤ T . (40) \operatorname{Re}\psi_1(s;\eta+i\xi)\leq p_\eta(s)<\mu_y^{-1},\qquad0\leq s\leq T. \tag*{(40)} Re ψ 1 ( s ; η + i ξ ) ≤ p η ( s ) < μ y − 1 , 0 ≤ s ≤ T . (40)
Proof. Write y = Im ψ 1 ( ⋅ ; η + i ξ ) y=\operatorname{Im}\psi_1(\mathord{\cdot};\eta+i\xi) y = Im ψ 1 ( ⋅ ; η + i ξ ) and d = p η − Re ψ 1 ( ⋅ ; η + i ξ ) d=p_\eta-\operatorname{Re}\psi_1(\mathord{\cdot};\eta+i\xi) d = p η − Re ψ 1 ( ⋅ ; η + i ξ ) . Subtracting the real part of Equation A.20 from Equation A.16 gives
d = σ V 2 2 κ R κ ∗ [ ( 2 p η − d ) d + y 2 ] . (41) d=\frac{\sigma_V^2}{2\kappa}R_\kappa*\left[(2p_\eta-d)d+y^2\right]. \tag*{(41)} d = 2 κ σ V 2 R κ ∗ [ ( 2 p η − d ) d + y 2 ] . (41)
Since p η p_\eta p η and d d d are continuous on [ 0 , T ] [0,T] [ 0 , T ] , choose λ ≥ 0 \lambda\geq0 λ ≥ 0 such that λ + ( σ V 2 / ( 2 κ ) ) ( 2 p η − d ) ≥ 0 \lambda+(\sigma_V^2/(2\kappa))(2p_\eta-d)\geq0 λ + ( σ V 2 / ( 2 κ )) ( 2 p η − d ) ≥ 0 throughout this interval. The resolvent identity in Equation A.1 gives ( I + λ R κ ∗ ) − 1 R κ ∗ = ( 1 + λ ) − 1 R κ ( 1 + λ ) ∗ (I+\lambda R_\kappa*)^{-1}R_\kappa*=(1+\lambda)^{-1}R_{\kappa(1+\lambda)}* ( I + λ R κ ∗ ) − 1 R κ ∗ = ( 1 + λ ) − 1 R κ ( 1 + λ ) ∗ . Adding λ R κ ∗ d \lambda R_\kappa*d λ R κ ∗ d to both sides of Equation A.22 therefore yields
d = 1 1 + λ R κ ( 1 + λ ) ∗ [ { λ + σ V 2 2 κ ( 2 p η − d ) } d + σ V 2 2 κ y 2 ] . (42) d=\frac{1}{1+\lambda}R_{\kappa(1+\lambda)}*\left[\left\{\lambda+\frac{\sigma_V^2}{2\kappa}(2p_\eta-d)\right\}d+\frac{\sigma_V^2}{2\kappa}y^2\right]. \tag*{(42)} d = 1 + λ 1 R κ ( 1 + λ ) ∗ [ { λ + 2 κ σ V 2 ( 2 p η − d ) } d + 2 κ σ V 2 y 2 ] . (42)
Treat the coefficient in braces as the fixed continuous function determined by the given solution. The resulting linear Volterra equation has a nonnegative kernel and nonnegative forcing. Boundedness of this coefficient and ∫ 0 δ R κ ( 1 + λ ) ( s ) d s → 0 \int_0^\delta R_{\kappa(1+\lambda)}(s)\,\mathrm{d}s\to0 ∫ 0 δ R κ ( 1 + λ ) ( s ) d s → 0 as δ ↓ 0 \delta\downarrow0 δ ↓ 0 make Picard iteration contractive on sufficiently short intervals of a common length. The iterates from zero are nonnegative on the first interval; on each subsequent interval, the already determined history contributes another nonnegative forcing term. Uniqueness on each interval identifies these nonnegative solutions with d d d , so d ≥ 0 d\geq0 d ≥ 0 on [ 0 , T ] [0,T] [ 0 , T ] . Together with Equation A.18, this proves Equation A.21.
Finally, subtracting Equation A.20 for two continuous solutions factors the difference of their squares into their sum times their difference. Both solutions are bounded and R κ ∈ L 1 ( 0 , T ) R_\kappa\in L^1(0,T) R κ ∈ L 1 ( 0 , T ) , so the same short-interval contraction argument proves uniqueness on [ 0 , T ] [0,T] [ 0 , T ] .
Proofs
Proof of Proposition 3.1
Proof. Fix T > 0 T>0 T > 0 and let m ( t ) = E Q [ V t ] m(t)=\mathbb{E}^{\mathbb{Q}}[V_t] m ( t ) = E Q [ V t ] for Lebesgue-almost every t ∈ [ 0 , T ] t\in[0,T] t ∈ [ 0 , T ] , taking a jointly measurable version. By Lemma A.2, conditional Jensen’s inequality and Fubini’s theorem give m ∈ L 2 ( 0 , T ) m\in L^2(0,T) m ∈ L 2 ( 0 , T ) almost surely. Since K ∈ L 2 ( 0 , T ) K\in L^2(0,T) K ∈ L 2 ( 0 , T ) and E Q [ ∫ 0 T V s d s ] < ∞ \mathbb{E}^{\mathbb{Q}}[\int_0^T V_s\,\mathrm{d}s]<\infty E Q [ ∫ 0 T V s d s ] < ∞ , Itô isometry shows that the Brownian convolution is square integrable for Lebesgue-almost every t t t . The compensated jump convolution is also square integrable by the compensated-Poisson isometry and E Q [ Y 1 2 ] = 2 μ y 2 \mathbb{E}^{\mathbb{Q}}[Y_1^2]=2\mu_y^2 E Q [ Y 1 2 ] = 2 μ y 2 . Both convolutions therefore have zero F 0 \mathcal{F}_0 F 0 -conditional mean. Conditional Fubini and E 0 Q [ θ s ] = θ 0 \mathbb{E}^{\mathbb{Q}}_0[\theta_s]=\theta_0 E 0 Q [ θ s ] = θ 0 give, Q ⊗ d t \mathbb{Q}\otimes dt Q ⊗ d t -almost everywhere,
m ( t ) − θ 0 = V 0 − θ 0 − κ ∫ 0 t K ( t − s ) [ m ( s ) − θ 0 ] d s . (43) m(t)-\theta_0=V_0-\theta_0-\kappa\int_0^t K(t-s)[m(s)-\theta_0]\,ds. \tag*{(43)} m ( t ) − θ 0 = V 0 − θ 0 − κ ∫ 0 t K ( t − s ) [ m ( s ) − θ 0 ] d s . (43)
Applying the resolvent identity in Equation A.1 to this pathwise L 1 ( 0 , T ) L^1(0,T) L 1 ( 0 , T ) equation yields, in the same almost-everywhere sense,
m ( t ) − θ 0 = ( V 0 − θ 0 ) [ 1 − ∫ 0 t R κ ( s ) d s ] = ( V 0 − θ 0 ) E α , 1 ( − κ t α ) , (44) m(t)-\theta_0=(V_0-\theta_0)\left[1-\int_0^t R_\kappa(s)\,ds\right]=(V_0-\theta_0)E_{\alpha,1}(-\kappa t^\alpha), \tag*{(44)} m ( t ) − θ 0 = ( V 0 − θ 0 ) [ 1 − ∫ 0 t R κ ( s ) d s ] = ( V 0 − θ 0 ) E α , 1 ( − κ t α ) , (44)
where the second equality follows from Lemma A.1 . Adding θ 0 \theta_0 θ 0 gives the right-hand side of Equation (10) , which is continuous on [ 0 , T ] [0,T] [ 0 , T ] and hence is the unique continuous representative of m m m .
Proof of Theorem 3.1
Proof. Fix t ∈ [ 0 , T ) t\in[0,T) t ∈ [ 0 , T ) and τ ∈ ( 0 , T − t ] \tau\in(0,T-t] τ ∈ ( 0 , T − t ] . Let q κ ( s ) = E α , 1 ( − κ s α ) q_\kappa(s)=E_{\alpha,1}(-\kappa s^\alpha) q κ ( s ) = E α , 1 ( − κ s α ) and M u = ∫ 0 u σ V V s d B 1 , s + J u y M_u=\int_0^u \sigma_V\sqrt{V_s}\,dB_{1,s}+J_u^y M u = ∫ 0 u σ V V s d B 1 , s + J u y . By Lemma A.2 , both θ \theta θ and M M M are square-integrable martingales. The resolvent calculation in the proof of that lemma, using Equation A.2, gives, Q ⊗ d u \mathbb{Q}\otimes du Q ⊗ d u -almost everywhere on [ 0 , T ] [0,T] [ 0 , T ] ,
V u = q κ ( u ) V 0 + ∫ 0 u R κ ( u − s ) θ s d s + 1 κ ∫ 0 u R κ ( u − s ) d M s . (45) V_u=q_\kappa(u)V_0+\int_0^u R_\kappa(u-s)\theta_s\,ds+\frac{1}{\kappa}\int_0^u R_\kappa(u-s)\,dM_s. \tag*{(45)} V u = q κ ( u ) V 0 + ∫ 0 u R κ ( u − s ) θ s d s + κ 1 ∫ 0 u R κ ( u − s ) d M s . (45)
For the Brownian term, Tonelli’s theorem and Cauchy–Schwarz give
∫ 0 τ ( E Q [ ∫ 0 t + r R κ ( t + r − s ) 2 V s d s ] ) 1 / 2 d r ≤ τ ∥ R κ ∥ L 2 ( 0 , T ) ( E Q [ ∫ 0 T V s d s ] ) 1 / 2 < ∞ . (46) \int_0^\tau\left(\mathbb{E}^{\mathbb{Q}}\left[\int_0^{t+r}R_\kappa(t+r-s)^2V_s\,ds\right]\right)^{1/2}dr
\leq\sqrt{\tau}\,\lVert R_\kappa\rVert_{L^2(0,T)}
\left(\mathbb{E}^{\mathbb{Q}}\left[\int_0^T V_s\,ds\right]\right)^{1/2}<\infty. \tag*{(46)} ∫ 0 τ ( E Q [ ∫ 0 t + r R κ ( t + r − s ) 2 V s d s ] ) 1/2 d r ≤ τ ∥ R κ ∥ L 2 ( 0 , T ) ( E Q [ ∫ 0 T V s d s ] ) 1/2 < ∞. (46)
This verifies the stochastic-Fubini condition for the window average. The finite-activity jump term admits pathwise Fubini, while R κ ∈ L 1 ( 0 , T ) R_\kappa\in L^1(0,T) R κ ∈ L 1 ( 0 , T ) and the first moments supplied by Lemma A.2 justify the deterministic and conditional Fubini interchanges. For Lebesgue-almost every r ∈ ( 0 , τ ) r\in(0,\tau) r ∈ ( 0 , τ ) , the stochastic integrals in Equation B.3 at u = t + r u=t+r u = t + r are square integrable. Splitting them at t t t , using E t Q [ θ s ] = θ t \mathbb{E}^{\mathbb{Q}}_t[\theta_s]=\theta_t E t Q [ θ s ] = θ t for s ≥ t s\geq t s ≥ t and the zero conditional mean of the future martingale integrals, gives almost surely
E t Q [ V t + r ] = q κ ( t + r ) V 0 + [ 1 − q κ ( r ) ] θ t + ∫ 0 t R κ ( t + r − s ) θ s d s + 1 κ ∫ 0 t R κ ( t + r − s ) d M s , (47) \mathbb{E}^{\mathbb{Q}}_t[V_{t+r}]=q_\kappa(t+r)V_0+[1-q_\kappa(r)]\theta_t+\int_0^tR_\kappa(t+r-s)\theta_s\,ds+\frac{1}{\kappa}\int_0^tR_\kappa(t+r-s)\,dM_s, \tag*{(47)} E t Q [ V t + r ] = q κ ( t + r ) V 0 + [ 1 − q κ ( r )] θ t + ∫ 0 t R κ ( t + r − s ) θ s d s + κ 1 ∫ 0 t R κ ( t + r − s ) d M s , (47)
where Lemma A.1 gives ∫ 0 r R κ ( s ) d s = 1 − q κ ( r ) \int_0^rR_\kappa(s)\,ds=1-q_\kappa(r) ∫ 0 r R κ ( s ) d s = 1 − q κ ( r ) . Averaging this almost-everywhere identity over r ∈ ( 0 , τ ) r\in(0,\tau) r ∈ ( 0 , τ ) and applying conditional and stochastic Fubini yields
V I X t , τ 2 = V 0 τ ∫ t t + τ q κ ( s ) d s + [ 1 − w ( τ ) ] θ t + ∫ 0 t a τ ( t − s ) { κ θ s d s + d M s } , (48) \mathrm{VIX}_{t,\tau}^{2}=\frac{V_0}{\tau}\int_t^{t+\tau}q_\kappa(s)\,\mathrm{d}s+[1-w(\tau)]\theta_t+\int_0^t a_\tau(t-s)\{\kappa\theta_s\,\mathrm{d}s+\mathrm{d}M_s\}, \tag*{(48)} VIX t , τ 2 = τ V 0 ∫ t t + τ q κ ( s ) d s + [ 1 − w ( τ )] θ t + ∫ 0 t a τ ( t − s ) { κ θ s d s + d M s } , (48)
because
w ( τ ) = 1 τ ∫ 0 τ q κ ( s ) d s = E α , 2 ( − κ τ α ) , w(\tau)=\frac{1}{\tau}\int_0^\tau q_\kappa(s)\,\mathrm{d}s=E_{\alpha,2}(-\kappa\tau^\alpha), w ( τ ) = τ 1 ∫ 0 τ q κ ( s ) d s = E α , 2 ( − κ τ α ) ,
a τ ( s ) = 1 κ τ ∫ s s + τ R κ ( v ) d v = q κ ( s ) − q κ ( s + τ ) κ τ . (49) a_\tau(s)=\frac{1}{\kappa\tau}\int_s^{s+\tau}R_\kappa(v)\,\mathrm{d}v=\frac{q_\kappa(s)-q_\kappa(s+\tau)}{\kappa\tau}. \tag*{(49)} a τ ( s ) = κ τ 1 ∫ s s + τ R κ ( v ) d v = κ τ q κ ( s ) − q κ ( s + τ ) . (49)
Moreover, the second identity in Equation B.7 gives
1 τ ∫ t t + τ q κ ( s ) d s = w ( τ ) − κ ∫ 0 t a τ ( s ) d s . (50) \frac{1}{\tau}\int_t^{t+\tau}q_\kappa(s)\,\mathrm{d}s=w(\tau)-\kappa\int_0^t a_\tau(s)\,\mathrm{d}s. \tag*{(50)} τ 1 ∫ t t + τ q κ ( s ) d s = w ( τ ) − κ ∫ 0 t a τ ( s ) d s . (50)
Since a τ a_\tau a τ is bounded and continuous on [ 0 , t ] [0,t] [ 0 , t ] , the averaged martingale integral is square integrable for the fixed pair ( t , τ ) (t,\tau) ( t , τ ) . Substituting Equation B.8 into Equation B.6 and expanding M M M therefore gives Equation (11) almost surely.
Proof of Theorem 3.2
Proof. Fix a real parameter u = η ∈ [ 0 , η ˉ T , τ 0 ) u=\eta\in[0,\bar{\eta}_{T,\tau_0}) u = η ∈ [ 0 , η ˉ T , τ 0 ) and take ψ 1 = p η \psi_1=p_\eta ψ 1 = p η and ψ 2 = ψ 2 ( ⋅ ; η ) \psi_2=\psi_2(\cdot;\eta) ψ 2 = ψ 2 ( ⋅ ; η ) from Lemma A.4 . Write d Z s = κ ( θ s − V s ) d s + σ V V s d B 1 , s + d J s y \mathrm{d}Z_s=\kappa(\theta_s-V_s)\,\mathrm{d}s+\sigma_V\sqrt{V_s}\,\mathrm{d}B_{1,s}+\mathrm{d}J_s^y d Z s = κ ( θ s − V s ) d s + σ V V s d B 1 , s + d J s y , so that V − V 0 = K ∗ d Z V-V_0=K*\mathrm{d}Z V − V 0 = K ∗ d Z almost everywhere in time. By Equation (11) ,
V I X T , τ 0 2 = w ( τ 0 ) V 0 + [ 1 − w ( τ 0 ) ] θ T + ∫ 0 T ( a τ 0 + κ K ∗ a τ 0 ) ( T − s ) d Z s . (51) \mathrm{VIX}_{T,\tau_0}^{2}=w(\tau_0)V_0+[1-w(\tau_0)]\theta_T+\int_0^T\left(a_{\tau_0}+\kappa K*a_{\tau_0}\right)(T-s)\,\mathrm{d}Z_s. \tag*{(51)} VIX T , τ 0 2 = w ( τ 0 ) V 0 + [ 1 − w ( τ 0 )] θ T + ∫ 0 T ( a τ 0 + κ K ∗ a τ 0 ) ( T − s ) d Z s . (51)
Here and below, boundedness of the deterministic coefficient functions, K ∈ L 1 ( 0 , T ) ∩ L 2 ( 0 , T ) K\in L^1(0,T)\cap L^2(0,T) K ∈ L 1 ( 0 , T ) ∩ L 2 ( 0 , T ) and Lemma A.2 justify the drift and Brownian Fubini interchanges; the finite-activity jump terms admit pathwise Fubini.
Applying I + κ K ∗ I+\kappa K* I + κ K ∗ to Equation A.16 and using R κ + κ K ∗ R κ = κ K R_\kappa+\kappa K*R_\kappa=\kappa K R κ + κ K ∗ R κ = κ K gives
ψ 1 = u ( a τ 0 + κ K ∗ a τ 0 ) + K ∗ ( − κ ψ 1 + σ V 2 2 ψ 1 2 ) . (52) \psi_1=u\left(a_{\tau_0}+\kappa K*a_{\tau_0}\right)+K*\left(-\kappa\psi_1+\frac{\sigma_V^2}{2}\psi_1^2\right). \tag*{(52)} ψ 1 = u ( a τ 0 + κ K ∗ a τ 0 ) + K ∗ ( − κ ψ 1 + 2 σ V 2 ψ 1 2 ) . (52)
Define the remaining coefficients by
χ 1 ( s ) = u w ( τ 0 ) + ∫ 0 s [ − κ ψ 1 ( v ) + σ V 2 2 ψ 1 ( v ) 2 ] d v , \chi_1(s)=uw(\tau_0)+\int_0^s\left[-\kappa\psi_1(v)+\frac{\sigma_V^2}{2}\psi_1(v)^2\right]\,\mathrm{d}v, χ 1 ( s ) = u w ( τ 0 ) + ∫ 0 s [ − κ ψ 1 ( v ) + 2 σ V 2 ψ 1 ( v ) 2 ] d v ,
ϕ ( s ) = λ y ∫ 0 s [ ϱ ( ψ 1 ( v ) ) − 1 − μ y ψ 1 ( v ) ] d v . (53) \phi(s)=\lambda^y\int_0^s[\varrho(\psi_1(v))-1-\mu_y\psi_1(v)]\,\mathrm{d}v. \tag*{(53)} ϕ ( s ) = λ y ∫ 0 s [ ϱ ( ψ 1 ( v )) − 1 − μ y ψ 1 ( v )] d v . (53)
Substituting Equation B.10 into Equation B.9 and using V − V 0 = K ∗ d Z V-V_0=K*\mathrm{d}Z V − V 0 = K ∗ d Z under stochastic Fubini yields
u VIX T , τ 0 2 = χ 1 ( T ) V 0 + u [ 1 − w ( τ 0 ) ] θ T + κ ∫ 0 T ψ 1 ( T − s ) θ s d s + ∫ 0 T σ V ψ 1 ( T − s ) V s d B 1 , s + ∫ 0 T ψ 1 ( T − s ) d J s y − σ V 2 2 ∫ 0 T ψ 1 ( T − s ) 2 V s d s . (54) u\operatorname{VIX}_{T,\tau_0}^{2}=\chi_1(T)V_0+u[1-w(\tau_0)]\theta_T+\kappa\int_0^T\psi_1(T-s)\theta_s\,\mathrm{d}s
+\int_0^T\sigma_V\psi_1(T-s)\sqrt{V_s}\,\mathrm{d}B_{1,s}+\int_0^T\psi_1(T-s)\,\mathrm{d}J_s^y
-\frac{\sigma_V^2}{2}\int_0^T\psi_1(T-s)^2V_s\,\mathrm{d}s. \tag*{(54)} u VIX T , τ 0 2 = χ 1 ( T ) V 0 + u [ 1 − w ( τ 0 )] θ T + κ ∫ 0 T ψ 1 ( T − s ) θ s d s + ∫ 0 T σ V ψ 1 ( T − s ) V s d B 1 , s + ∫ 0 T ψ 1 ( T − s ) d J s y − 2 σ V 2 ∫ 0 T ψ 1 ( T − s ) 2 V s d s . (54)
For the mean factor, apply Itô’s formula to ψ 2 ( T − s ) θ s \psi_2(T-s)\theta_s ψ 2 ( T − s ) θ s . Since d θ s = σ θ θ s d B 2 , s \mathrm{d}\theta_s=\sigma_\theta\sqrt{\theta_s}\,\mathrm{d}B_{2,s} d θ s = σ θ θ s d B 2 , s and ψ 2 \psi_2 ψ 2 satisfies Equation A.17,
u [ 1 − w ( τ 0 ) ] θ T + κ ∫ 0 T ψ 1 ( T − s ) θ s d s = ψ 2 ( T ) θ 0 + ∫ 0 T σ θ ψ 2 ( T − s ) θ s d B 2 , s − σ θ 2 2 ∫ 0 T ψ 2 ( T − s ) 2 θ s d s . (55) u[1-w(\tau_0)]\theta_T+\kappa\int_0^T\psi_1(T-s)\theta_s\,\mathrm{d}s
=\psi_2(T)\theta_0+\int_0^T\sigma_\theta\psi_2(T-s)\sqrt{\theta_s}\,\mathrm{d}B_{2,s}
-\frac{\sigma_\theta^2}{2}\int_0^T\psi_2(T-s)^2\theta_s\,\mathrm{d}s. \tag*{(55)} u [ 1 − w ( τ 0 )] θ T + κ ∫ 0 T ψ 1 ( T − s ) θ s d s = ψ 2 ( T ) θ 0 + ∫ 0 T σ θ ψ 2 ( T − s ) θ s d B 2 , s − 2 σ θ 2 ∫ 0 T ψ 2 ( T − s ) 2 θ s d s . (55)
Combining these identities and exponentiating gives
exp { u VIX T , τ 0 2 } = exp { ϕ ( T ) + χ 1 ( T ) V 0 + ψ 2 ( T ) θ 0 } E T , (56) \exp\left\{u\operatorname{VIX}_{T,\tau_0}^{2}\right\}
=\exp\{\phi(T)+\chi_1(T)V_0+\psi_2(T)\theta_0\}\mathcal{E}_T, \tag*{(56)} exp { u VIX T , τ 0 2 } = exp { ϕ ( T ) + χ 1 ( T ) V 0 + ψ 2 ( T ) θ 0 } E T , (56)
where, for 0 ≤ t ≤ T 0\le t\le T 0 ≤ t ≤ T ,
E t = exp { ∫ 0 t σ V ψ 1 ( T − s ) V s d B 1 , s + ∫ 0 t σ θ ψ 2 ( T − s ) θ s d B 2 , s − 1 2 ∫ 0 t [ σ V 2 ψ 1 ( T − s ) 2 V s + σ θ 2 ψ 2 ( T − s ) 2 θ s ] d s + ∫ 0 t ψ 1 ( T − s ) d J s y − λ y ∫ 0 t [ ϱ ( ψ 1 ( T − s ) ) − 1 − μ y ψ 1 ( T − s ) ] d s } . (57) \mathcal{E}_t=\exp\left\{
\int_0^t\sigma_V\psi_1(T-s)\sqrt{V_s}\,\mathrm{d}B_{1,s}
+\int_0^t\sigma_\theta\psi_2(T-s)\sqrt{\theta_s}\,\mathrm{d}B_{2,s}
-\frac{1}{2}\int_0^t\left[\sigma_V^2\psi_1(T-s)^2V_s+\sigma_\theta^2\psi_2(T-s)^2\theta_s\right]\mathrm{d}s
+\int_0^t\psi_1(T-s)\,\mathrm{d}J_s^y
-\lambda^y\int_0^t[\varrho(\psi_1(T-s))-1-\mu_y\psi_1(T-s)]\mathrm{d}s
\right\}. \tag*{(57)} E t = exp { ∫ 0 t σ V ψ 1 ( T − s ) V s d B 1 , s + ∫ 0 t σ θ ψ 2 ( T − s ) θ s d B 2 , s − 2 1 ∫ 0 t [ σ V 2 ψ 1 ( T − s ) 2 V s + σ θ 2 ψ 2 ( T − s ) 2 θ s ] d s + ∫ 0 t ψ 1 ( T − s ) d J s y − λ y ∫ 0 t [ ϱ ( ψ 1 ( T − s )) − 1 − μ y ψ 1 ( T − s )] d s } . (57)
We first prove the martingale property for real u = η ∈ [ 0 , η ˉ T , τ 0 ) u=\eta\in[0,\bar{\eta}_{T,\tau_0}) u = η ∈ [ 0 , η ˉ T , τ 0 ) , taking ψ 1 = p η \psi_1=p_\eta ψ 1 = p η and ψ 2 = ψ 2 ( ⋅ ; η ) \psi_2=\psi_2(\mathord{\cdot};\eta) ψ 2 = ψ 2 ( ⋅ ; η ) from Lemma A.4 . Write E t = D t exp { L t − ⟨ L ⟩ t / 2 } \mathcal{E}_t=D_t\exp\{L_t-\langle L\rangle_t/2\} E t = D t exp { L t − ⟨ L ⟩ t /2 } , where L L L is the sum of the two Brownian integrals and D D D is the jump exponential in Equation B.15. The continuous and jump factors have zero quadratic covariation, so E \mathcal{E} E is a nonnegative local martingale. Since B 1 B_1 B 1 and B 2 B_2 B 2 are independent,
⟨ L ⟩ T = ∫ 0 T [ σ V 2 p η ( T − s ) 2 V s + σ θ 2 ψ 2 ( T − s ; η ) 2 θ s ] d s . \langle L\rangle_T=\int_0^T[\sigma_V^2p_\eta(T-s)^2V_s+\sigma_\theta^2\psi_2(T-s;\eta)^2\theta_s]\,\mathrm{d}s. ⟨ L ⟩ T = ∫ 0 T [ σ V 2 p η ( T − s ) 2 V s + σ θ 2 ψ 2 ( T − s ; η ) 2 θ s ] d s .
Choose a > 0 a>0 a > 0 with a σ θ 2 T ≤ 1 a\sigma_\theta^2T\le1 a σ θ 2 T ≤ 1 . The bound for E 0 Q [ e a θ s ] \mathbb{E}_0^{\mathbb{Q}}[e^{a\theta_s}] E 0 Q [ e a θ s ] in the proof of Lemma A.3 and Jensen’s inequality give E 0 Q [ exp { ( a / T ) ∫ 0 T θ s d s } ] ≤ e 2 a θ 0 \mathbb{E}_0^{\mathbb{Q}}[\exp\{(a/T)\int_0^T\theta_s\,\mathrm{d}s\}]\le e^{2a\theta_0} E 0 Q [ exp {( a / T ) ∫ 0 T θ s d s }] ≤ e 2 a θ 0 . Both coefficient functions converge uniformly to zero as η ↓ 0 \eta\downarrow0 η ↓ 0 , so reduce η ˉ T , τ 0 \bar{\eta}_{T,\tau_0} η ˉ T , τ 0 until 2 σ V 2 ∥ p η ∥ ∞ 2 < c T 2\sigma_V^2\|p_\eta\|_\infty^2<c_T 2 σ V 2 ∥ p η ∥ ∞ 2 < c T , 2 σ θ 2 ∥ ψ 2 ( ⋅ ; η ) ∥ ∞ 2 < a / T 2\sigma_\theta^2\|\psi_2(\mathord{\cdot};\eta)\|_\infty^2<a/T 2 σ θ 2 ∥ ψ 2 ( ⋅ ; η ) ∥ ∞ 2 < a / T and 2 ∥ p η ∥ ∞ < μ y − 1 2\|p_\eta\|_\infty<\mu_y^{-1} 2∥ p η ∥ ∞ < μ y − 1 throughout this interval, with norms on [ 0 , T ] [0,T] [ 0 , T ] . Conditional Cauchy–Schwarz and Lemma A.3 then yield E 0 Q [ e ⟨ L ⟩ T ] < ∞ \mathbb{E}_0^{\mathbb{Q}}[e^{\langle L\rangle_T}]<\infty E 0 Q [ e ⟨ L ⟩ T ] < ∞ almost surely.
The exponential formula for the marked Poisson process shows that D D D is a strictly positive martingale with E 0 Q [ D T ] = 1 \mathbb{E}_0^{\mathbb{Q}}[D_T]=1 E 0 Q [ D T ] = 1 and
E 0 Q [ D T 2 ] = exp { λ y ∫ 0 T [ ϱ ( 2 p η ( s ) ) − 2 ϱ ( p η ( s ) ) + 1 ] d s } < ∞ . (58) \mathbb{E}_0^{\mathbb{Q}}[D_T^2]=\exp\left\{\lambda^y\int_0^T[\varrho(2p_\eta(s))-2\varrho(p_\eta(s))+1]\,\mathrm{d}s\right\}<\infty. \tag*{(58)} E 0 Q [ D T 2 ] = exp { λ y ∫ 0 T [ ϱ ( 2 p η ( s )) − 2 ϱ ( p η ( s )) + 1 ] d s } < ∞. (58)
Define Q y \mathbb{Q}^{y} Q y on F T \mathcal{F}_{T} F T by d Q y / d Q = D T \mathrm{d}\mathbb{Q}^{y}/\mathrm{d}\mathbb{Q}=D_{T} d Q y / d Q = D T . This change of measure leaves B 1 , B 2 B_{1}, B_{2} B 1 , B 2 as a two-dimensional Brownian motion: the density martingale has zero quadratic covariation with each Brownian motion, so Girsanov’s theorem introduces no continuous drift. Moreover, Q y \mathbb{Q}^{y} Q y agrees with Q \mathbb{Q} Q on F 0 \mathcal{F}_{0} F 0 , and conditional Cauchy–Schwarz gives
E 0 Q y [ e ⟨ L ⟩ T / 2 ] = E 0 Q [ D T e ⟨ L ⟩ T / 2 ] ≤ ( E 0 Q [ D T 2 ] ) 1 / 2 ( E 0 Q [ e ⟨ L ⟩ T ] ) 1 / 2 < ∞ . (59) \mathbb{E}_{0}^{\mathbb{Q}^{y}}\left[e^{\langle L\rangle_{T}/2}\right]=\mathbb{E}_{0}^{\mathbb{Q}}\left[D_{T}e^{\langle L\rangle_{T}/2}\right]\leq\left(\mathbb{E}_{0}^{\mathbb{Q}}\left[D_{T}^{2}\right]\right)^{1/2}\left(\mathbb{E}_{0}^{\mathbb{Q}}\left[e^{\langle L\rangle_{T}}\right]\right)^{1/2}<\infty. \tag*{(59)} E 0 Q y [ e ⟨ L ⟩ T /2 ] = E 0 Q [ D T e ⟨ L ⟩ T /2 ] ≤ ( E 0 Q [ D T 2 ] ) 1/2 ( E 0 Q [ e ⟨ L ⟩ T ] ) 1/2 < ∞. (59)
To apply Novikov’s criterion, first condition Q y \mathbb{Q}^{y} Q y on the F 0 \mathcal{F}_{0} F 0 -measurable events { E 0 Q y [ e ⟨ L ⟩ T / 2 ] ≤ n } \{\mathbb{E}_{0}^{\mathbb{Q}^{y}}[e^{\langle L\rangle_{T}/2}]\leq n\} { E 0 Q y [ e ⟨ L ⟩ T /2 ] ≤ n } of positive probability. The Brownian property is preserved and the required unconditional expectation is then finite. Novikov’s criterion applies to the joint continuous exponential under each such conditional measure. These events increase to the whole space up to a null set, giving E 0 Q y [ exp { L T − ⟨ L ⟩ T / 2 } ] = 1 \mathbb{E}_{0}^{\mathbb{Q}^{y}}[\exp\{L_{T}-\langle L\rangle_{T}/2\}]=1 E 0 Q y [ exp { L T − ⟨ L ⟩ T /2 }] = 1 . Bayes’ formula therefore yields E 0 Q [ E T ] = 1 \mathbb{E}_{0}^{\mathbb{Q}}[\mathcal{E}_{T}]=1 E 0 Q [ E T ] = 1 , so E \mathcal{E} E is a uniformly integrable martingale on [ 0 , T ] [0,T] [ 0 , T ] . When λ y = 0 \lambda^{y}=0 λ y = 0 , take D ≡ 1 D\equiv1 D ≡ 1 and Q y = Q \mathbb{Q}^{y}=\mathbb{Q} Q y = Q ; the same argument applies. Taking conditional expectations in Equation B.14 therefore gives f ( η ; 0 , T ) = exp { ϕ ( T ) + χ 1 ( T ) V 0 + ψ 2 ( T ) θ 0 } < ∞ f(\eta;0,T)=\exp\{\phi(T)+\chi_{1}(T)V_{0}+\psi_{2}(T)\theta_{0}\}<\infty f ( η ; 0 , T ) = exp { ϕ ( T ) + χ 1 ( T ) V 0 + ψ 2 ( T ) θ 0 } < ∞ almost surely for every sufficiently small real η ≥ 0 \eta\geq0 η ≥ 0 , with the coefficients specified above.
It remains to extend the identity to complex arguments. By the definition of squared VIX, conditional Jensen’s inequality, and Lemma A.3, after reducing η ˉ T , τ 0 \bar{\eta}_{T,\tau_{0}} η ˉ T , τ 0 once more if necessary,
E 0 Q [ e η VIX T , τ 0 2 ] ≤ E 0 Q [ exp { η τ 0 ∫ T T + τ 0 V s d s } ] ≤ E 0 Q [ exp { η τ 0 ∫ 0 T + τ 0 V s d s } ] < ∞ . (60) \begin{aligned}
\mathbb{E}_{0}^{\mathbb{Q}}\left[e^{\eta\operatorname{VIX}_{T,\tau_{0}}^{2}}\right]\leq\mathbb{E}_{0}^{\mathbb{Q}}\left[\exp\left\{\frac{\eta}{\tau_{0}}\int_{T}^{T+\tau_{0}}V_{s}\,\mathrm{d}s\right\}\right] \\
&\leq\mathbb{E}_{0}^{\mathbb{Q}}\left[\exp\left\{\frac{\eta}{\tau_{0}}\int_{0}^{T+\tau_{0}}V_{s}\,\mathrm{d}s\right\}\right]<\infty. \tag*{(60)}
\end{aligned} E 0 Q [ e η VIX T , τ 0 2 ] ≤ E 0 Q [ exp { τ 0 η ∫ T T + τ 0 V s d s } ] ≤ E 0 Q [ exp { τ 0 η ∫ 0 T + τ 0 V s d s } ] < ∞. (60)
for every η ∈ ( 0 , η ˉ T , τ 0 ) \eta\in(0,\bar{\eta}_{T,\tau_{0}}) η ∈ ( 0 , η ˉ T , τ 0 ) . Fix u = η + i ξ u=\eta+i\xi u = η + i ξ with 0 ≤ η < η ˉ T , τ 0 0\leq\eta<\bar{\eta}_{T,\tau_{0}} 0 ≤ η < η ˉ T , τ 0 and ξ ∈ R \xi\in\mathbb{R} ξ ∈ R . Since R κ ∈ L 1 ( 0 , T ) R_{\kappa}\in L^{1}(0,T) R κ ∈ L 1 ( 0 , T ) and the nonlinearity is polynomial, contraction on short intervals gives a unique maximal continuous solution of Equation A.20. The comparison argument in Lemma A.5 applies on every compact interval of existence and gives Re ψ 1 ( s ; u ) ≤ p η ( s ) \operatorname{Re}\psi_{1}(s;u)\leq p_{\eta}(s) Re ψ 1 ( s ; u ) ≤ p η ( s ) .
To exclude explosion, write D C α g ( s ) = Γ ( 1 − α ) − 1 ∫ 0 s ( s − v ) − α g ′ ( v ) d v D_{C}^{\alpha}g(s)=\Gamma(1-\alpha)^{-1}\int_{0}^{s}(s-v)^{-\alpha}g'(v)\,\mathrm{d}v D C α g ( s ) = Γ ( 1 − α ) − 1 ∫ 0 s ( s − v ) − α g ′ ( v ) d v for the Caputo derivative of an absolutely continuous function. The resolvent identity rewrites the equation as ψ 1 = u a τ 0 + K ∗ [ − κ ψ 1 + κ u a τ 0 + ( σ V 2 / 2 ) ψ 1 2 ] \psi_{1}=ua_{\tau_{0}}+K\ast[-\kappa\psi_{1}+\kappa ua_{\tau_{0}}+(\sigma_{V}^{2}/2)\psi_{1}^{2}] ψ 1 = u a τ 0 + K ∗ [ − κ ψ 1 + κ u a τ 0 + ( σ V 2 /2 ) ψ 1 2 ] . On each compact interval of existence, boundedness first gives α \alpha α -Hölder continuity of ψ 1 \psi_{1} ψ 1 and of the integrand. Since 2 α > 1 2\alpha>1 2 α > 1 , differentiating the convolution of this integrand with K K K gives an integrable bound of the form C ( s α − 1 + s 2 α − 1 ) C(s^{\alpha-1}+s^{2\alpha-1}) C ( s α − 1 + s 2 α − 1 ) near zero. Together with absolute continuity of a τ 0 a_{\tau_{0}} a τ 0 , this proves that ψ 1 \psi_{1} ψ 1 is absolutely continuous on that interval. Its Caputo equation is therefore
D C α ψ 1 ( s ; u ) = u h ( s ) − κ ψ 1 ( s ; u ) + σ V 2 2 ψ 1 ( s ; u ) 2 , ψ 1 ( 0 ; u ) = u a τ 0 ( 0 ) , D_{C}^{\alpha}\psi_{1}(s;u)=uh(s)-\kappa\psi_{1}(s;u)+\frac{\sigma_{V}^{2}}{2}\psi_{1}(s;u)^{2},\qquad\psi_{1}(0;u)=ua_{\tau_{0}}(0), D C α ψ 1 ( s ; u ) = u h ( s ) − κ ψ 1 ( s ; u ) + 2 σ V 2 ψ 1 ( s ; u ) 2 , ψ 1 ( 0 ; u ) = u a τ 0 ( 0 ) ,
h ( s ) = − E α , 1 ( − κ ( s + τ 0 ) α ) τ 0 + 1 κ τ 0 Γ ( 1 − α ) ∫ 0 s ( s − v ) − α R κ ( v + τ 0 ) d v . (61) h(s)=-\frac{E_{\alpha,1}(-\kappa(s+\tau_{0})^{\alpha})}{\tau_{0}}+\frac{1}{\kappa\tau_{0}\Gamma(1-\alpha)}\int_{0}^{s}(s-v)^{-\alpha}R_{\kappa}(v+\tau_{0})\,\mathrm{d}v. \tag*{(61)} h ( s ) = − τ 0 E α , 1 ( − κ ( s + τ 0 ) α ) + κ τ 0 Γ ( 1 − α ) 1 ∫ 0 s ( s − v ) − α R κ ( v + τ 0 ) d v . (61)
Here h = D C α a τ 0 + κ a τ 0 h=D_{C}^{\alpha}a_{\tau_{0}}+\kappa a_{\tau_{0}} h = D C α a τ 0 + κ a τ 0 , and its displayed expression follows from Lemma A.1 . Since τ 0 > 0 \tau_{0}>0 τ 0 > 0 , the function R κ ( ⋅ + τ 0 ) R_{\kappa}(\cdot+\tau_{0}) R κ ( ⋅ + τ 0 ) is bounded on [ 0 , T ] [0,T] [ 0 , T ] , so h h h is bounded as well.
Fractional integration by parts gives D C α ∣ ψ 1 ∣ 2 ≤ 2 Re ( ψ 1 ‾ D C α ψ 1 ) D^\alpha_C|\psi_1|^2 \le2\operatorname{Re}(\overline{\psi_1}D^\alpha_C\psi_1) D C α ∣ ψ 1 ∣ 2 ≤ 2 Re ( ψ 1 D C α ψ 1 ) . Indeed, the difference between the right-hand side and the left-hand side equals
1 Γ ( 1 − α ) [ ∣ ψ 1 ( s ; u ) − ψ 1 ( 0 ; u ) ∣ 2 s α + α ∫ 0 s ∣ ψ 1 ( s ; u ) − ψ 1 ( v ; u ) ∣ 2 ( s − v ) 1 + α d v ] ≥ 0. (62) \frac{1}{\Gamma(1-\alpha)}\left[\frac{|\psi_1(s;u)-\psi_1(0;u)|^2}{s^\alpha}+\alpha\int_0^s\frac{|\psi_1(s;u)-\psi_1(v;u)|^2}{(s-v)^{1+\alpha}}\,\mathrm{d}v\right]\ge0. \tag*{(62)} Γ ( 1 − α ) 1 [ s α ∣ ψ 1 ( s ; u ) − ψ 1 ( 0 ; u ) ∣ 2 + α ∫ 0 s ( s − v ) 1 + α ∣ ψ 1 ( s ; u ) − ψ 1 ( v ; u ) ∣ 2 d v ] ≥ 0. (62)
Consequently, the real-part comparison yields
D C α ∣ ψ 1 ( s ; u ) ∣ 2 ≤ 2 Re ( ψ 1 ( s ; u ) ‾ u h ( s ) ) + [ − 2 κ + σ V 2 Re ψ 1 ( s ; u ) ] ∣ ψ 1 ( s ; u ) ∣ 2 ≤ ∣ u ∣ 2 ∥ h ∥ ∞ 2 + ( 1 + σ V 2 ∥ p η ∥ ∞ ) ∣ ψ 1 ( s ; u ) ∣ 2 . (63) \begin{aligned}
D^\alpha_C|\psi_1(s;u)|^2 &\le2\operatorname{Re}\left(\overline{\psi_1(s;u)}uh(s)\right)+\left[-2\kappa+\sigma_V^2\operatorname{Re}\psi_1(s;u)\right]|\psi_1(s;u)|^2\\
&\le|u|^2\|h\|_\infty^2+\left(1+\sigma_V^2\|p_\eta\|_\infty\right)|\psi_1(s;u)|^2. \tag*{(63)}
\end{aligned} D C α ∣ ψ 1 ( s ; u ) ∣ 2 ≤ 2 Re ( ψ 1 ( s ; u ) u h ( s ) ) + [ − 2 κ + σ V 2 Re ψ 1 ( s ; u ) ] ∣ ψ 1 ( s ; u ) ∣ 2 ≤ ∣ u ∣ 2 ∥ h ∥ ∞ 2 + ( 1 + σ V 2 ∥ p η ∥ ∞ ) ∣ ψ 1 ( s ; u ) ∣ 2 . (63)
Applying K ∗ K* K ∗ and iterating the resulting linear inequality gives
∣ ψ 1 ( s ; u ) ∣ 2 ≤ ∣ u ∣ 2 [ a τ 0 ( 0 ) 2 + ∥ h ∥ ∞ 2 T α Γ ( 1 + α ) ] E α , 1 ( [ 1 + σ V 2 ∥ p η ∥ ∞ ] T α ) (64) |\psi_1(s;u)|^2\le|u|^2\left[a_{\tau_0}(0)^2+\frac{\|h\|_\infty^2T^\alpha}{\Gamma(1+\alpha)}\right]E_{\alpha,1}\left(\left[1+\sigma_V^2\|p_\eta\|_\infty\right]T^\alpha\right) \tag*{(64)} ∣ ψ 1 ( s ; u ) ∣ 2 ≤ ∣ u ∣ 2 [ a τ 0 ( 0 ) 2 + Γ ( 1 + α ) ∥ h ∥ ∞ 2 T α ] E α , 1 ( [ 1 + σ V 2 ∥ p η ∥ ∞ ] T α ) (64)
throughout its interval of existence before T T T . This bound is independent of the endpoint of that interval. Boundedness of ψ 1 2 \psi_1^2 ψ 1 2 and translation continuity of R κ R_\kappa R κ in L 1 L^1 L 1 give a continuous limit of the convolution at any finite endpoint. The resulting continuous history term permits another short-interval contraction. Thus ψ 1 \psi_1 ψ 1 extends uniquely to [ 0 , T ] [0,T] [ 0 , T ] , with Re ψ 1 ( s ; u ) ≤ p η ( s ) < μ y − 1 \operatorname{Re}\psi_1(s;u)\le p_\eta(s)<\mu_y^{-1} Re ψ 1 ( s ; u ) ≤ p η ( s ) < μ y − 1 .
Given this solution, consider the complex Riccati equation
ψ 2 ′ ( s ; u ) = κ ψ 1 ( s ; u ) + σ θ 2 2 ψ 2 ( s ; u ) 2 , ψ 2 ( 0 ; u ) = u [ 1 − w ( τ 0 ) ] . (65) \psi_2'(s;u)=\kappa\psi_1(s;u)+\frac{\sigma_\theta^2}{2}\psi_2(s;u)^2,\qquad\psi_2(0;u)=u[1-w(\tau_0)]. \tag*{(65)} ψ 2 ′ ( s ; u ) = κ ψ 1 ( s ; u ) + 2 σ θ 2 ψ 2 ( s ; u ) 2 , ψ 2 ( 0 ; u ) = u [ 1 − w ( τ 0 )] . (65)
Its right-hand side is locally Lipschitz, so it has a unique maximal solution. Subtracting its real part from Equation A.17 gives a linear equation for ψ 2 ( s ; η ) − Re ψ 2 ( s ; u ) \psi_2(s;\eta)-\operatorname{Re}\psi_2(s;u) ψ 2 ( s ; η ) − Re ψ 2 ( s ; u ) , with zero initial value and nonnegative forcing κ [ p η ( s ) − Re ψ 1 ( s ; u ) ] + ( σ θ 2 / 2 ) [ Im ψ 2 ( s ; u ) ] 2 \kappa[p_\eta(s)-\operatorname{Re}\psi_1(s;u)]+(\sigma_\theta^2/2)[\operatorname{Im}\psi_2(s;u)]^2 κ [ p η ( s ) − Re ψ 1 ( s ; u )] + ( σ θ 2 /2 ) [ Im ψ 2 ( s ; u ) ] 2 . An integrating factor therefore gives Re ψ 2 ( s ; u ) ≤ ψ 2 ( s ; η ) \operatorname{Re}\psi_2(s;u)\le\psi_2(s;\eta) Re ψ 2 ( s ; u ) ≤ ψ 2 ( s ; η ) before any possible explosion. Hence
d d s ∣ ψ 2 ( s ; u ) ∣ 2 = 2 κ Re ( ψ 2 ( s ; u ) ‾ ψ 1 ( s ; u ) ) + σ θ 2 Re ψ 2 ( s ; u ) ∣ ψ 2 ( s ; u ) ∣ 2 ≤ κ ∣ ψ 1 ( s ; u ) ∣ 2 + [ κ + σ θ 2 ∥ ψ 2 ( ⋅ ; η ) ∥ ∞ ] ∣ ψ 2 ( s ; u ) ∣ 2 . (66) \begin{aligned}
\frac{\mathrm{d}}{\mathrm{d}s}|\psi_2(s;u)|^2&=2\kappa\operatorname{Re}\left(\overline{\psi_2(s;u)}\psi_1(s;u)\right)+\sigma_\theta^2\operatorname{Re}\psi_2(s;u)|\psi_2(s;u)|^2\\
&\le\kappa|\psi_1(s;u)|^2+\left[\kappa+\sigma_\theta^2\|\psi_2(\cdot;\eta)\|_\infty\right]|\psi_2(s;u)|^2. \tag*{(66)}
\end{aligned} d s d ∣ ψ 2 ( s ; u ) ∣ 2 = 2 κ Re ( ψ 2 ( s ; u ) ψ 1 ( s ; u ) ) + σ θ 2 Re ψ 2 ( s ; u ) ∣ ψ 2 ( s ; u ) ∣ 2 ≤ κ ∣ ψ 1 ( s ; u ) ∣ 2 + [ κ + σ θ 2 ∥ ψ 2 ( ⋅ ; η ) ∥ ∞ ] ∣ ψ 2 ( s ; u ) ∣ 2 . (66)
The preceding bound for ψ 1 \psi_1 ψ 1 and Gronwall’s inequality exclude explosion of ψ 2 \psi_2 ψ 2 before T T T . Thus both coefficient functions exist uniquely on [ 0 , T ] [0,T] [ 0 , T ] throughout the stated strip.
These bounds are uniform when u u u ranges over a compact subset of the strip. Picard iteration on successive sufficiently short intervals therefore converges locally uniformly in u u u . Since the integral equations are polynomial in their unknown functions and in u u u , their iterates are holomorphic, proving that u ↦ ψ 1 ( ⋅ ; u ) u\mapsto\psi_1(\cdot;u) u ↦ ψ 1 ( ⋅ ; u ) and u ↦ ψ 2 ( ⋅ ; u ) u\mapsto\psi_2(\cdot;u) u ↦ ψ 2 ( ⋅ ; u ) are holomorphic with values in C ( [ 0 , T ] ) C([0,T]) C ([ 0 , T ]) on 0 < Re ( u ) < η ˉ T , τ 0 0<\operatorname{Re}(u)<\bar\eta_{T,\tau_0} 0 < Re ( u ) < η ˉ T , τ 0 . The same construction gives continuous dependence at Re ( u ) = 0 \operatorname{Re}(u)=0 Re ( u ) = 0 . Define χ 1 \chi_1 χ 1 and ϕ \phi ϕ by the integrals in the real-parameter part of the proof. Their holomorphic dependence follows from that of ψ 1 \psi_1 ψ 1 and the bound Re ψ 1 < μ y − 1 \operatorname{Re}\psi_1<\mu_y^{-1} Re ψ 1 < μ y − 1 .
Finally, take the conditional moment bounds in Equation B.18 and the real-parameter transform identities simultaneously on a countable dense subset of ( 0 , η ˉ T , τ 0 ) (0,\bar{\eta}_{T,\tau_0}) ( 0 , η ˉ T , τ 0 ) . Outside a common null set, a conditional distribution of V I X T , τ 0 2 \mathrm{VIX}^{2}_{T,\tau_0} VIX T , τ 0 2 defines a holomorphic version of f ( u ; 0 , T ) f(u;0,T) f ( u ; 0 , T ) on the open strip: on each compact subset, an exponential moment at a larger real parameter dominates the integrand and its derivatives. The affine expression is holomorphic on the same strip and agrees with this version on the chosen dense set. The identity theorem therefore proves Equation (14) throughout the open strip, with the coefficients specified above. Conditional dominated convergence and continuity of the coefficient functions extend the identity to Re ( u ) = 0 \operatorname{Re}(u)=0 Re ( u ) = 0 .
Proof of Corollary 3.1
Proof. For Re ( z ) > 0 \operatorname{Re}(z)>0 Re ( z ) > 0 , the Laplace transform of the payoff x ↦ ( x − K ) + x\mapsto(\sqrt{x}-K)^{+} x ↦ ( x − K ) + is
∫ 0 ∞ e − z x ( x − K ) + d x = π 2 z 3 / 2 erfc ( K z ) . (67) \int_{0}^{\infty}e^{-zx}(\sqrt{x}-K)^{+}\,\mathrm{d}x=\frac{\sqrt{\pi}}{2z^{3/2}}\operatorname{erfc}(K\sqrt{z}). \tag*{(67)} ∫ 0 ∞ e − z x ( x − K ) + d x = 2 z 3/2 π erfc ( K z ) . (67)
Indeed, this follows by substituting x = y 2 x=y^{2} x = y 2 and integrating the resulting Gaussian terms.
Fix η ∈ ( 0 , η ˉ T , τ 0 ) \eta\in(0,\bar{\eta}_{T,\tau_0}) η ∈ ( 0 , η ˉ T , τ 0 ) . By Theorem 3.2 , the entire contour u = η + i ξ u=\eta+i\xi u = η + i ξ , ξ ∈ R \xi\in\mathbb{R} ξ ∈ R , lies in the transform domain, with f ( η ; 0 , T ) < ∞ f(\eta;0,T)<\infty f ( η ; 0 , T ) < ∞ and ∣ f ( η + i ξ ; 0 , T ) ∣ ≤ f ( η ; 0 , T ) \lvert f(\eta+i\xi;0,T)\rvert\le f(\eta;0,T) ∣ f ( η + i ξ ; 0 , T )∣ ≤ f ( η ; 0 , T ) almost surely. For z = η + i ξ z=\eta+i\xi z = η + i ξ , two integrations by parts in Equation B.25 give
π 2 z 3 / 2 erfc ( K z ) = e − z K 2 2 K z 2 − 1 4 z 2 ∫ K 2 ∞ e − z x x − 3 / 2 d x , (68) \frac{\sqrt{\pi}}{2z^{3/2}}\operatorname{erfc}(K\sqrt{z})=\frac{e^{-zK^{2}}}{2Kz^{2}}-\frac{1}{4z^{2}}\int_{K^{2}}^{\infty}e^{-zx}x^{-3/2}\,\mathrm{d}x, \tag*{(68)} 2 z 3/2 π erfc ( K z ) = 2 K z 2 e − z K 2 − 4 z 2 1 ∫ K 2 ∞ e − z x x − 3/2 d x , (68)
∣ π 2 z 3 / 2 erfc ( K z ) ∣ ≤ e − η K 2 K ( η 2 + ξ 2 ) . \left\lvert\frac{\sqrt{\pi}}{2z^{3/2}}\operatorname{erfc}(K\sqrt{z})\right\rvert\le\frac{e^{-\eta K^{2}}}{K(\eta^{2}+\xi^{2})}. 2 z 3/2 π erfc ( K z ) ≤ K ( η 2 + ξ 2 ) e − η K 2 .
The bound uses K > 0 K>0 K > 0 and ∫ K 2 ∞ x − 3 / 2 d x = 2 / K \int_{K^{2}}^{\infty}x^{-3/2}\,\mathrm{d}x=2/K ∫ K 2 ∞ x − 3/2 d x = 2/ K , and proves integrability of the payoff transform along the contour.
The function x ↦ e − η x ( x − K ) + x\mapsto e^{-\eta x}(\sqrt{x}-K)^{+} x ↦ e − η x ( x − K ) + on [ 0 , ∞ ) [0,\infty) [ 0 , ∞ ) , extended by zero to negative x x x , is continuous and integrable. Fourier inversion therefore gives its representation at every x x x , and multiplication by e η x e^{\eta x} e η x yields the Bromwich representation of the payoff. To interchange this integral with conditional expectation, restrict first to the F 0 \mathcal{F}_{0} F 0 -measurable events { f ( η ; 0 , T ) ≤ n } \{f(\eta;0,T)\le n\} { f ( η ; 0 , T ) ≤ n } . After multiplication by the corresponding indicator, the expected integral of the absolute inversion integrand is bounded by n π e − η K 2 / ( K η ) < ∞ n\pi e^{-\eta K^{2}}/(K\eta)<\infty nπ e − η K 2 / ( K η ) < ∞ . Ordinary Fubini applies on these events, which increase to the whole space up to a null set. Thus Bromwich inversion and conditional Fubini give
E 0 Q [ ( V I X T , τ 0 − K ) + ] = 1 4 π ∫ − ∞ ∞ erfc ( K η + i ξ ) ( η + i ξ ) 3 / 2 f ( η + i ξ ; 0 , T ) d ξ = 1 2 π ∫ 0 ∞ Re [ erfc ( K η + i ξ ) ( η + i ξ ) 3 / 2 f ( η + i ξ ; 0 , T ) ] d ξ . (69) \begin{aligned}
\mathbb{E}_{0}^{\mathbb{Q}}\left[(\mathrm{VIX}_{T,\tau_0}-K)^{+}\right]=\frac{1}{4\sqrt{\pi}}\int_{-\infty}^{\infty}\frac{\operatorname{erfc}(K\sqrt{\eta+i\xi})}{(\eta+i\xi)^{3/2}}f(\eta+i\xi;0,T)\,\mathrm{d}\xi\\
&=\frac{1}{2\sqrt{\pi}}\int_{0}^{\infty}\operatorname{Re}\left[\frac{\operatorname{erfc}(K\sqrt{\eta+i\xi})}{(\eta+i\xi)^{3/2}}f(\eta+i\xi;0,T)\right]\mathrm{d}\xi. \tag*{(69)}
\end{aligned} E 0 Q [ ( VIX T , τ 0 − K ) + ] = 4 π 1 ∫ − ∞ ∞ ( η + i ξ ) 3/2 erfc ( K η + i ξ ) f ( η + i ξ ; 0 , T ) d ξ = 2 π 1 ∫ 0 ∞ Re [ ( η + i ξ ) 3/2 erfc ( K η + i ξ ) f ( η + i ξ ; 0 , T ) ] d ξ . (69)
The second equality follows from conjugate symmetry. Multiplying by e − r T e^{-rT} e − r T gives Equation (16) . ◻️