A scaled version of the double-mean-reverting model for VIX derivatives

Abstract

As the Heston model is not consistent with VIX data in real market well enough, alternative stochastic volatility models including the double-mean-reverting model of Gatheral (in: Bachelier Congress, 2008) have been developed to overcome its limitation. The double-mean-reverting model is a three factor model successfully reflecting the empirical dynamics of the variance but there is no closed form solution for VIX derivatives and SPX options and thus calibration using conventional techniques may be slow. In this paper, we propose a fast mean-reverting version of the double-mean-reverting model. We obtain a closed form approximation for VIX derivatives and show how it is effective by comparing it with the Heston model and the double-mean-reverting model.

Appendix

This appendix is devoted to prove the price approximation (3) with the approximation error, stated in Sect. 3, using the asymptotic analysis of Fouque et al. [12]. Without loss of generality, the interest rate r is assumed to zero here (so, $D\left( \tau \right)$ is omitted), which makes it possible to price VIX futures and options consistently. Once pricing formulas are derived, the original prices can be obtained just by discounting the results.

First, applying the Feynman–Kac theorem (cf. [20]) to $P^{\epsilon }$, we obtain a singularly perturbed partial differential equation given by

$$\begin{aligned} \left( \frac{1}{\epsilon }\mathcal {L}_{0}+\frac{1}{\sqrt{\epsilon }} \mathcal {L}_{1}+\mathcal {L}_{2}\right) P^{\epsilon }(t,y,z)=0, \end{aligned}$$

where

$$\begin{aligned} \mathcal {L}_{0}&=\left( z-y\right) \partial _{y}+\nu ^{2}y^{2\alpha }\partial _{yy}^{2},\\ \mathcal {L}_{1}&=\sqrt{2}\rho \nu \sigma _{z}y^{\alpha }\sqrt{z}\partial _{yz}^{2},\\ \mathcal {L}_{2}&=\partial _{t}+c\left( \theta -z\right) \partial _{z}+\frac{1}{2}\sigma _{z}^{2}z\partial _{zz}^{2}. \end{aligned}$$

Then the expansion $P^{\epsilon }=P_{0}+\sqrt{\epsilon }P_{1}+\epsilon P_{2}+\cdots$ gives

$$\begin{aligned} \frac{1}{\epsilon }\mathcal {L}_{0}P_{0}&+\frac{1}{\sqrt{\epsilon }}\left( \mathcal {L}_{0}P_{1}+\mathcal {L}_{1}P_{0}\right) +\left( \mathcal {L}_{0}P_{2} +\mathcal {L}_{1}P_{1}+\mathcal {L}_{2}P_{0}\right) \nonumber \\&+\sqrt{\epsilon }\left( \mathcal {L}_{0}P_{3}+\mathcal {L}_{1}P_{2}+\mathcal {L}_{2} P_{1}\right) +\epsilon \left( \mathcal {L}_{0}P_{4}+\mathcal {L}_{1}P_{3} +\mathcal {L}_{2}P_{2}\right) +\cdots =0. \end{aligned}$$

(5)

We have the equation $\mathcal {L}_{0}P_{0}=0$ from the $\frac{1}{\epsilon }$-order term. This equation leads to the following solution (which depends on $\alpha$):

$$\begin{aligned} P_{0}=c_{1}\left( t,z\right) \int \exp \left\{ \frac{1}{\left( 2-2\alpha \right) \nu ^{2}}y^{2-2\alpha }-\frac{z}{\left( 1-2\alpha \right) \nu ^{2}}y^{1-2\alpha }\right\} dy+c_{2}\left( t,z\right) \end{aligned}$$

for $1/2<\alpha <1$,

$$\begin{aligned} P_{0}=c_{1}\left( t,z\right) \int \frac{\mathrm{e}^{y/\nu ^{2}}}{y^{z/\nu ^{2}}}dy+c_{2}\left( t,z\right) \end{aligned}$$

for $\alpha =1/2$, and

$$\begin{aligned} P_{0}=c_{1}\left( t,z\right) \int y^{\frac{1}{\nu ^{2}}}\mathrm{e}^{z/\nu ^{2}y}dy+c_{2}\left( t,z\right) \end{aligned}$$

for $\alpha =1$. Then, from the growth condition (cf. Sect. 3) assumed for the VIX derivative price (2), the above $c_{1}(t,z)$ must be zero. So, $P_{0}$ is independent of y. Similarly, the equation from the $\frac{1}{\sqrt{\epsilon }}$-order terms gives a result that $P_{1}$ is also independent of y. Then the following Poisson equation is obtained from the $\epsilon ^{0}$-order terms:

$$\begin{aligned} \mathcal {L}_{0}P_{2}+\mathcal {L}_{1}P_{1}+\mathcal {L}_{2}P_{0} =\mathcal {L}_{0}P_{2}+\mathcal {L}_{2}P_{0}=0. \end{aligned}$$

(6)

Its solvability condition (cf. [12]) is $\left\langle \mathcal {L}_{2}P_{0}\right\rangle = \left\langle \mathcal {L}_{2}\right\rangle P_{0} = 0$, where $\left\langle \,.\,\right\rangle$ denotes expectation with respect to the invariant distribution of $Y_{t}$ (cf. Remark below for the explicit formula of the invariant distribution). Since $P_{0}$ and the operator $\mathcal {L}_{2}$ do not depend on y, $\left\langle \mathcal {L}_{2}\right\rangle P_{0} =0$ leads to

$$\begin{aligned} \mathcal {L}_{2}P_{0}=\left( \partial _{t}+c(\theta -z)\partial _{z} +\frac{1}{2}\sigma _{2}^{2}z\partial _{zz}^{2}\right) P_{0}=0. \end{aligned}$$

(7)

Now, we substitute the expansion of $P^{\epsilon }$ into (2) and obtain

$$\begin{aligned} \mathrm{E}^{Q}\left[ \left. h^{\epsilon }\left( Y_{T},Z_{T}\right) \right| Y_{t} =y,Z_{t}=z\right] =P_{0}\left( t,z\right) +\sqrt{\epsilon }P_{1}\left( t,z\right) +\cdots . \end{aligned}$$

(8)

Letting $\epsilon$ go to zero, we obtain

$$\begin{aligned} P_{0}\left( t,z\right)&=\underset{\epsilon \rightarrow 0}{\lim }\mathrm{E}^{Q} \left[ \left. h^{\epsilon }\left( Y_{T},Z_{T}\right) \right| Y_{t}=y,Z_{t}=z\right] \nonumber \\&=\underset{\epsilon \rightarrow 0}{\lim }\int _{0}^{\infty }\int _{0}^{\infty } h^{\epsilon }\left( u,v\right) \phi _{y,z}\left( u|Z_{T}=v\right) du\,\psi _{z} \left( v\right) dv\nonumber \\&=\int _{0}^{\infty }\int _{0}^{\infty }\underset{\epsilon \rightarrow 0}{\lim } \left\{ h^{\epsilon }\left( u,v\right) \phi _{y,z}\left( u|Z_{T}=v\right) \right\} du\,\psi _{z}\left( v\right) dv\nonumber \\&=\int _{0}^{\infty }\int _{0}^{\infty }\tilde{h}_{0}\left( u,v\right) \phi _{v}^{*} \left( u\right) du\,\psi _{z}\left( v\right) dv\nonumber \\&=\int _{0}^{\infty }h_{0}\left( v\right) \psi _{z}\left( v\right) dv\nonumber \\&=\mathrm{E}^{Q^{*}}\left[ \left. h_{0}\left( Z_{T}\right) \right| Z_{t}=z\right] , \end{aligned}$$

(9)

where $\tilde{h}_{0}:=\underset{\epsilon \rightarrow 0}{\lim }h^{\epsilon }$, $h_{0}:=\int _{0}^{\infty }\tilde{h}_{0}\left( u,\cdot \right) \phi _{v}^{*}\left( u\right) du$, $\phi _{y,z}\left( \,\cdot \,|Z_{T}=v\right)$ is the conditional density function of $Y_{T}$ given $Z_{T}=v$, $Y_{t}=y$ and $Z_{t}=z$, $\phi _{v}^{*}$ is the invariant density function of $Y_{T}$ given $Z_{T}=v$, and $\psi _{z}$ is the density function of $Z_{T}$ given $Z_{t}=z$. Here, the bounded convergence theorem (cf. Folland [11]) has been applied to change the limit and the integral. Knowing that the operator $\mathcal {L}_{2}-\partial _{t}$ is the infinitesimal operator of the CIR process $Z_{t}$, the Feynman–Kac theorem says that Eq. (9) implies the terminal condition $P_{0}\left( T,z\right) =h_{0}\left( z\right)$. Also, one can express $P_{0}$ as a single integral of the Heston price type integrand given by

$$\begin{aligned} P_{0}(t,z)=\int _{0}^{\infty }h_{0}\left( \gamma v\right) f(v;k,\lambda )\,dv, \end{aligned}$$

where

$$\begin{aligned} \gamma =\frac{\left( 1-\mathrm{e}^{-c\tau }\right) \sigma _{z}^{2}}{4c} \end{aligned}$$

and $f\left( \nu ;k,\lambda \right)$ is a probability density function of non-central chi-square distribution with $k=\frac{4c\theta }{\sigma _{z}^{2}}$ degrees of freedom and the non-centrality parameter $\lambda =\frac{z\mathrm{e}^{-c\tau }}{\gamma }$. Note that

$$\begin{aligned} \tilde{h}_{0}\left( u,v\right) =\underset{\epsilon \rightarrow 0}{\lim } h^{\epsilon }\left( u,v\right) =\underset{\epsilon \rightarrow 0}{\lim }H\left( \sqrt{a_{1}u+a_{2}v+a_{3}\theta }\times 100\right) =H\left( \sqrt{a_{2}^{*}v+a_{3}^{*}\theta }\times 100\right) , \end{aligned}$$

or

$$\begin{aligned} h_{0}\left( Z_{T}\right) =H\left( \mathrm{VIX}_{T}^{*}\right) , \end{aligned}$$

where $\mathrm{VIX}_{T}^{*}=\sqrt{a_{2}^{*}Z_{T}+a_{3}^{*}\theta }\times 100$, $a_{2}^{*}=\frac{1}{c\varDelta T}\left( 1-\mathrm{e}^{-c\varDelta T}\right)$ and $a_{3}^{*}=1-a_{2}^{*}$. This means that $P_{0}$ is exactly the same as the VIX derivative price of the Heston model.

In what follows, we calculate the first-order correction term $P_{1}$. Because $\mathcal {L}_{0}P_{2}=0$ is obtained from (6) and (7), $P_{2}$ is independent of y also. Then, applying the solvability condition to the equation , which comes from the $\sqrt{\epsilon }$-order terms, gives

$$\begin{aligned} \mathcal {L}_{2}P_{1}=\left( \partial _{t}+c(\theta -z)\partial _{z}+\frac{1}{2}\sigma _{2}^{2}z\partial _{zz}^{2}\right) P_{1}=0. \end{aligned}$$

The terminal condition for this PDE is found by utilizing (8). If subtracting $P_{0}$ from the both sides of equation (8) and dividing it with $\sqrt{\epsilon }$, we have

$$\begin{aligned}&\frac{1}{\sqrt{\epsilon }}\left( \mathrm{E}^{Q}\left[ \left. h^{\epsilon }\left( Y_{T},Z_{T}\right) \right| Y_{t}=y,Z_{t}=z\right] -\mathrm{E}^{Q^{*}}\left[ \left. h_{0}\left( Z_{T}\right) \right| Z_{t}=z\right] \right) \\&\quad =P_{1}\left( t,z\right) +\sqrt{\epsilon }P_{2}\left( t,z\right) +\cdots \end{aligned}$$

Letting $\epsilon$ go to zero, we obtain

$$\begin{aligned} P_{1}\left( t,z\right)&=\underset{\epsilon \rightarrow 0}{\lim } \frac{1}{\sqrt{\epsilon }}\left[ \mathrm{E}^{Q}\left[ \left. h^{\epsilon } \left( Y_{T},Z_{T}\right) \right| Y_{t}=y,Z_{t}=z\right] -\mathrm{E}^{Q^{*}} \left[ \left. h_{0}\left( Z_{T}\right) \right| Z_{t}=z\right] \right] \nonumber \\&=\underset{\epsilon \rightarrow 0}{\lim }\frac{1}{\sqrt{\epsilon }} \left[ \int _{0}^{\infty }\int _{0}^{\infty }h^{\epsilon }\left( u,v\right) \phi _{y,z}\left( u|Z_{T}=v\right) du\,\psi _{z}\left( v\right) dv-\int _{0}^{\infty } h_{0}\left( v\right) \psi _{z}\left( v\right) dv\right] \nonumber \\&=\underset{\epsilon \rightarrow 0}{\lim }\int _{0}^{\infty }\int _{0}^{\infty } \left( \frac{h^{\epsilon }\left( u,v\right) -\tilde{h}_{0}\left( u,v\right) }{\sqrt{\epsilon }} \right) \phi _{y,z}\left( u|Z_{T}=v\right) du\,\psi _{z}\left( v\right) dv\nonumber \\&=\int _{0}^{\infty }\int _{0}^{\infty }\underset{\epsilon \rightarrow 0}{\lim } \left\{ \left( \frac{h^{\epsilon }\left( u,v\right) -\tilde{h}_{0}\left( u,v\right) }{\sqrt{\epsilon }}\right) \phi _{y,z}\left( u|Z_{T}=v\right) \right\} du\,\psi _{z} \left( v\right) dv\nonumber \\&=\int _{0}^{\infty }\int _{0}^{\infty }\tilde{h}_{1}\left( u,v\right) \phi _{v}^{*} \left( u\right) du\,\psi _{z}\left( v\right) dv\nonumber \\&=\int _{0}^{\infty }h_{1}\left( v\right) \psi _{z}\left( v\right) dv\nonumber \\&=\mathrm{E}^{Q^{*}}\left[ \left. h_{1}\left( Z_{T}\right) \right| Z_{t}=z\right] , \end{aligned}$$

(10)

where $\tilde{h}_{1}:=\underset{\epsilon \rightarrow 0}{\lim }\frac{h^{\epsilon }\left( u,v\right) -\tilde{h}_{0}\left( u,v\right) }{\sqrt{\epsilon }}$ and $h_{1}:=\int _{0}^{\infty }\tilde{h}_{1}\left( u,\cdot \right) \phi _{v}^{*}\left( u\right) du$. Then we have

$$\begin{aligned} P_{1}(t,z)=\int _{0}^{\infty }h_{1}\left( \gamma v\right) f(v;k,\lambda )\,dv. \end{aligned}$$

In fact, one can extend the results for $P_0$ and $P_1$ to all $P_{i}$ for $i\ge 2$ by induction and obtain

$$\begin{aligned} P_{i}(t,z)=\int _{0}^{\infty }h_{i}\left( \gamma v\right) f(v;k,\lambda )\,dv \quad (i\ge 2) \end{aligned}$$

and they all are independent of y because $\mathcal {L}_{0}P_{i}=0$ holds. At first glance, it might look strange because VIX derivatives are certainly affected by the process $Y_t$. However, the perturbation theory used in our model framework may not reflect the current value y of the fast mean-reverting variance directly. Its mean level z (the intermediate-term variance) is a major factor of the price. Recall $\mathrm{VIX}_{T}^{2}=a_{1}Y_{T}+a_{2}Z_{T}+a_{3}\theta$. For the calibration purpose, this is desirable in the sense that it is difficult to observe the short-term value of volatility whereas it may not be the case with its mean level. Even if there is no explicit dependence of each term $P_i$ on y, we note that the dependence on the parameters $\alpha$ and $\nu$ controlling the process $Y_t$ can be identified through the invariant distribution $\phi ^*_z$ which is given by Remark at the end of this appendix.

Now, we explain how to obtain each $h_{i}$. We first note that $\tilde{h}_{i}$ is defined in such a way that $h^{\epsilon }=\tilde{h}_{0}+\sqrt{\epsilon }\tilde{h}_{1}+\epsilon \tilde{h}_{2}+\cdots$ holds, which leads to the idea that $\tilde{h}_{i}$ could be easily derived by deploying the Talyor expansion associated with $h^{\epsilon }$. By the definitions of $a_{1}$ and $a_{2}$ in Sect, 2, $a_{1}u+a_{2}v+a_{3}\theta =a_{2}^{*}v+a_{3}^{*}\theta +\left( ca_{2}^{*}\left( v-\theta \right) +\frac{1}{\varDelta T}\left( u-v\right) \right) \epsilon +O\left( \epsilon ^{2}\right)$ for some variables u and v. So, for a given (regularized) payoff H, $h^{\epsilon }\left( u,v\right)$ can be expanded as follows:

$$\begin{aligned} h^{\epsilon }\left( u,v\right)&=H\left( 100\sqrt{a_{1}u+a_{2}v+a_{3}\theta }\right) \\&=\tilde{H}\left( a_{1}u+a_{2}v+a_{3}\theta \right) \\&=\tilde{H}\left( a_{2}^{*}v+a_{3}^{*}\theta \right) +\tilde{H}'\left( a_{2}^{*}v+a_{3}^{*}\theta \right) \left( ca_{2}^{*}\left( v-\theta \right) +\frac{1}{\varDelta T}\left( u-v\right) \right) \epsilon +\cdots , \end{aligned}$$

where $\tilde{H}\left( x\right) :=H\left( 100\sqrt{x}\right)$ is a single variable function. This implies

$$\begin{aligned} \tilde{h}_{0}\left( u,v\right)= & {} \tilde{H}\left( a_{2}^{*}v+a_{3}^{*}\theta \right) ,\\ \tilde{h}_{2}\left( u,v\right)= & {} \tilde{H}'\left( a_{2}^{*}v+a_{3}^{*} \theta \right) \left( ca_{2}^{*}\left( v-\theta \right) +\frac{1}{\varDelta T}\left( u-v\right) \right) , \end{aligned}$$

and $\tilde{h}_{i}\left( u,v\right) =0$ for every odd number i. So, we obtain $h_{0}$ and $h_{2}$ given by

$$\begin{aligned} h_{0}\left( v\right)&=\int _{0}^{\infty }\tilde{h}_{0}\left( u,v\right) \phi _{v}^{*}\left( u\right) du\\&=\int _{0}^{\infty }\tilde{H}\left( a_{2}^{*}v+a_{3}^{*}\theta \right) \phi _{v}^{*}\left( u\right) du\\&=\tilde{H}\left( a_{2}^{*}v+a_{3}^{*}\theta \right) , \end{aligned}$$

$$\begin{aligned} h_{2}\left( v\right)&=\int _{0}^{\infty }\tilde{H}'\left( a_{2}^{*} v+a_{3}^{*}\theta \right) \left( ca_{2}^{*}\left( v-\theta \right) +\frac{1}{\varDelta T} \left( u-v\right) \right) \phi _{v}^{*}\left( u\right) du\\&=ca_{2}^{*}\left( v-\theta \right) \tilde{H}'\left( a_{2}^{*}v+a_{3}^{*}\theta \right) , \end{aligned}$$

respectively and $h_{i}\left( v\right) =0$ for every odd number i. Here, $\int _{0}^{\infty }u\phi _{v}^{*}\left( u\right) du=v$ has been used for $h_{2}$. It is notable that $P_{i}=0$ for every odd number i. Particularly, $P_{1}=0$ and $P_{3}=0$ improve the accuracy of our approximation. Specific forms of $h_{i}$ for each VIX derivative are obtained by applying this method to the corresponding payoff H given by

• VIX futures: $H\left( x\right) =x$

• VIX call option: $H\left( x\right) =\left( x-K\right) \varvec{1}_{\left\{ x\ge K\right\} }$

• VIX put option: $H\left( x\right) =\left( K-x\right) \varvec{1}_{\left\{ x\le K\right\} }$

• VIX binary call option: $H\left( x\right) =\varvec{1}_{\left\{ x\ge K\right\} }$

• VIX binary put option: $H\left( x\right) =\varvec{1}_{\left\{ x\le K\right\} }$

Here, K is the strike of an option.

Next, we discuss the accuracy of the price approximation $\tilde{P}^{\epsilon }\left( =P_{0}+\epsilon P_{2}\right)$ for the futures. One can do it similarly for the options. For convenience, let $\hat{P}^{\epsilon }:=\tilde{P}^{\epsilon }+\epsilon ^{2}P_{4}+\epsilon ^{2}\sqrt{\epsilon }P_{5}$, $R^{\epsilon }:=P^{\epsilon }-\hat{P}^{\epsilon }$ and $\mathcal {L}^{\epsilon }:=\frac{1}{\epsilon }\mathcal {L}_{0}+\frac{1}{\sqrt{\epsilon }}\mathcal {L}_{1}+\mathcal {L}_{2}$. Then

$$\begin{aligned}&\mathcal {L}^{\epsilon }R^{\epsilon }+\mathcal {L}^{\epsilon }\left( P_{0}+\epsilon P_{2}+\epsilon ^{2}P_{4}+\epsilon ^{2}\sqrt{\epsilon }P_{5}\right) \\&\quad =\mathcal {L}^{\epsilon }R^{\epsilon }+\frac{1}{\epsilon }\mathcal {L}_{0} P_{0}+\frac{1}{\sqrt{\epsilon }}\mathcal {L}_{1}P_{0}+\left( \mathcal {L}_{0}P_{2} +\mathcal {L}_{2}P_{0}\right) +\sqrt{\epsilon }\mathcal {L}_{1}P_{2}+\epsilon \left( \mathcal {L}_{0}P_{4}+\mathcal {L}_{2}P_{2}\right) \\&\qquad +\epsilon \sqrt{\epsilon }\left( \mathcal {L}_{0}P_{5}+\mathcal {L}_{1} P_{4}\right) +\epsilon ^{2}\left( \mathcal {L}_{1}P_{5}+\mathcal {L}_{2}P_{4}\right) +\epsilon ^{2}\sqrt{\epsilon }\mathcal {L}_{2}P_{5}=0 \end{aligned}$$

and so from (5) we have

$$\begin{aligned} \mathcal {L}^{\epsilon }R^{\epsilon }+\epsilon ^{2}\left( \mathcal {L}_{1}P_{5} +\mathcal {L}_{2}P_{4}\right) +\epsilon ^{2}\sqrt{\epsilon }\mathcal {L}_{2}P_{5}=0, \end{aligned}$$

where $R^{\epsilon }$ has the terminal condition

$$\begin{aligned} R^{\epsilon }\left( T,y,z\right)&=P^{\epsilon }\left( T,y,z\right) -\left( P_{0}\left( T,z\right) +\epsilon P_{2}\left( T,z\right) +\epsilon ^{2} P_{4}\left( T,z\right) +\epsilon ^{2}\sqrt{\epsilon }P_{5}\left( T,z\right) \right) . \end{aligned}$$

Then, from the Feynman-Kac theorem (probabilistic representation), we have

$$\begin{aligned} R^{\epsilon }\left( t,y,z\right)&=\mathrm{E}^{Q}\left[ \left. P^{\epsilon } \left( T,Y_{T},Z_{T}\right) -P_{0}\left( T,Z_{T}\right) -\epsilon P_{2} \left( T,Z_{T}\right) \right| Y_{t}=y,Z_{t}=z\right] \\&+\epsilon ^{2}\mathrm{E}^{Q}\left[ -P_{4}\left( T,Y_{T},Z_{T}\right) +\int _{t}^{T}\left\{ \mathcal {L}_{1}P_{5}\left( s,Y_{s},Z_{s}\right) +\mathcal {L}_{2}P_{4}\left( s,Y_{s},Z_{s}\right) \right\} \right. \\&\qquad \left. \left. ds\right| Y_{t}=y,Z_{t} =z\right] \\&+\epsilon ^{2}\sqrt{\epsilon }\mathrm{E}^{Q}\left[ \left. -P_{5} \left( T,Y_{T},Z_{T}\right) +\int _{t}^{T}\mathcal {L}_{2}P_{5}\left( s,Y_{s},Z_{s} \right) ds\right| Y_{t}=y,Z_{t}=z\right] . \end{aligned}$$

Here, the first term of the right hand side is convergent in the sense that

$$\begin{aligned}&\underset{\epsilon \rightarrow 0}{\lim }\frac{1}{\epsilon ^{2}}\mathrm{E}^{Q} \left[ \left. P^{\epsilon }\left( T,Y_{T},Z_{T}\right) -P_{0}\left( T,Z_{T}\right) -\epsilon P_{2}\left( T,Z_{T}\right) \right| Y_{t}=y,Z_{t}=z\right] \\&\quad =E^{Q^{*}}\left[ \left. \frac{1}{8h_{0}(Z_{T})^{3}}\left\{ b_{1}^{*} \theta +b_{2}^{*}\left( Z_{T}-\theta \right) \right\} \left( Z_{T}-\theta \right) \right| Z_{t}=z\right] \\&\qquad -\underset{\epsilon \rightarrow 0}{\lim }E^{Q}\left[ \left. \frac{1}{8\tau _{0}^{2}h_{0}(Z_{T})^{3}}\left( Y_{T}^{2}-Z_{T}^{2}\right) \right| Y_{t} =y,Z_{t}=z\right] , \end{aligned}$$

where $b_{1}^{*}$ and $b_{2}^{*}$ are given by

$$\begin{aligned} b_{1}^{*}&=4c^{2}a_{2}^{*}-\frac{4c}{\varDelta T},\\ b_{2}^{*}&=\left( b_{1}^{*}-c^{2}a_{2}^{*}\right) a_{2}^{*}, \end{aligned}$$

respectively, and

$$\begin{aligned}&\underset{\epsilon \rightarrow 0}{\lim }E^{Q}\left[ \left. \frac{1}{8\tau _{0}^{2} h_{0}(Z_{T})^{3}}\left( Y_{T}^{2}-Z_{T}^{2}\right) \right| Y_{t}=y,Z_{t}=z\right] \\&\quad =\underset{\epsilon \rightarrow 0}{\lim }\int _{0}^{\infty }\int _{0}^{\infty } \left( \frac{1}{8\tau _{0}^{2}h_{0}(v)^{3}}\left( u^{2}-v^{2}\right) \right) \phi _{y,z}\left( u|Z_{T}=v\right) du\,\psi _{z}\left( v\right) dv\\&\quad =\int _{0}^{\infty }\int _{0}^{\infty }\left( u^{2}-v^{2}\right) \phi _{z}^{*}\left( u\right) du\,\frac{1}{8\tau _{0}^{2}h_{0}(v)^{3}} \psi _{z}\left( v\right) dv \; < \; C_{1}\left( z\right) \end{aligned}$$

for some bounded function $C_{1}$. Here, all the moments of $Y_{T}$ with respect to invariant distribution of $Y_{t}$ are assumed to be well-defined. Hence, there is a bounded function $C_{2}\left( t,y,z\right)$ independent of $\epsilon$ such that

$$\begin{aligned} \mathrm{E}^{Q}\left[ \left. P^{\epsilon }\left( T,Y_{T},Z_{T}\right) -P_{0} \left( T,Z_{T}\right) -\epsilon P_{2}\left( T,Z_{T}\right) \right| Y_{t}=y,Z_{t} =z\right] \le C_{2}\left( t,y,z\right) \epsilon ^{2}. \end{aligned}$$

Consequently, $\left| R^{\epsilon }\left( t,y,z\right) \right| \le C_{3}\left( t,y,z\right) \epsilon ^{2}$ for some bounded function $C_{3}\left( t,y,z\right)$ independent of $\epsilon$ and thus

$$\begin{aligned} \left| P^{\epsilon }\left( t,y,z\right) -\tilde{P}^{\epsilon }\left( t,z\right) \right| \le \left| R^{\epsilon }\left( t,y,z\right) \right| +\left| \hat{P}^{\epsilon }\left( t,y,z\right) -\tilde{P}^{\epsilon }\left( t,z\right) \right| \le C\left( t,y,z\right) \epsilon ^{2} \end{aligned}$$

for some bounded function $C\left( t,y,z\right)$ independent of $\epsilon$.

Remark

The invariant distributions $\phi ^*_z$ of the mean-reverting process $Y_t$ given $Z_{t}=z$ for various $\alpha$ are as follows.

• $\alpha =1$:

  $$\begin{aligned} \phi ^*_z\left( y\right) =\frac{\zeta ^{\eta }}{\Gamma \left( \eta \right) } \frac{e^{-\zeta /y}}{y^{1+\eta }}, \end{aligned}$$

  which is a density of the inverse gamma distribution with shape parameter $\eta =1+1/\nu ^{2}$ and scale parameter $\zeta =z/\nu ^{2}$.

• $\alpha =\frac{1}{2}$:

  $$\begin{aligned} \phi ^*_z\left( y\right) =\frac{1}{\nu ^{2\theta }\Gamma \left( \theta \right) } y^{\theta -1}e^{-y/\nu ^{2}}, \end{aligned}$$

  which is a density of the gamma distribution with shape parameter $\theta =z/\nu ^{2}$ and scale parameter $\nu ^{2}$.

• $\frac{1}{2}<\alpha <1$:

  $$\begin{aligned} \phi ^*_z\left( y\right) =\frac{1}{Cy^{2\alpha }}\mathrm{exp}\left( \frac{z}{\nu ^{2}\left( 1-2\alpha \right) }y^{1-2\alpha }-\frac{1}{\nu ^{2}\left( 2-2\alpha \right) }y^{2-2\alpha }\right) , \end{aligned}$$

  where C is a constant such that $\int _{0}^{\infty }\phi ^*_z=1$ holds.