3.1. Fractional stochastic volatility model

Let $(\Omega ,\mathcal {F},\mathbb {P})$ be a probability space. We investigate an option pricing problem when the dynamics of the underlying are driven by stochastic volatility by the following fractional SDE:

(3.1) $$ \begin{align} \left \{ \begin{aligned}\, dS_t&=\mu{S_t}\,dt+{|v_t|}S_t\,dB^S_t,\\ dv_t&=\beta(\alpha-v_t)\,dt+{\gamma_1}d{B}^H_t+\gamma_2\,dB^v_t, \end{aligned} \right. \end{align} $$

where $S_t$ is risky asset price process. Here, $\mu $ is the drift rate of the risk asset price process and $\beta $ is the average recurrent rate of the volatility process. Since the volatility process is a mean-reverting process, $v_t$ tends towards a long-term value $\alpha $ with rate $\beta $ . Here, $\gamma _1$ and $\gamma _2$ are two constants. We suppose that $B^S_t$ and $B^v_t$ are two sBms with correlation coefficient $\rho $ . $B^H_t$ is an fBm with Hurst parameter $H>1/2$ , and it is independent of $B^S_t$ and $B^v_t$ . In order to derive the PDE, we need the following lemma.

Lemma 3.1. The equation of volatility has a unique solution of the form

(3.2) $$ \begin{align} v_t=e^{-\beta t}v_0+\beta\alpha\int^t_0e^{\beta(s-t)}\,ds+\int^t_0\gamma_2 e^{\beta(s-t)}\,dB^v_s+\gamma_1B^H_t-\beta\int^t_0\gamma_1 e^{\beta(s-t)}B^H_s\,ds. \end{align} $$

The Malliavin derivatives of the items with fBm of equation (3.2) are equal to

$$ \begin{align*} \mathbb{D}^H_uv_t=\bigg(\gamma_1Z_H(t,u)-\beta\int^t_0\gamma_1 e^{\beta(s-t)}Z_H(s,u)\,ds\bigg)\mathbf{1}_{\{u<t\}}, \end{align*} $$

where $Z_H(\cdot ,\cdot )$ is defined in (2.1).

We now consider a contingent claim $U(t,S_t,v_t)$ under model (3.1). To obtain the PDE of U, we try to set up a portfolio to replicate it. On the other hand, we let the market price of volatility risk be given by the function $\lambda (t,S_t,v_t)$ . We assume that $\lambda (t,S_t,v_t)$ is independent of both the underlying $S_t$ and the time t [Reference Breeden5, Reference Cox, Ingersoll and Ross9]. The consumption process is supposed to be a Cox–Ingersoll–Ross process [Reference Cox, Ingersoll and Ross9] and the investors have log-utility. Then, by Breeden’s consumption-based model [Reference Breeden5], $\lambda (t,S_t,v_t)$ is proportional to the volatility, that is, $\lambda (t,S_t,v_t)=\lambda |v|$ , where $\lambda $ is a positive constant.

Theorem 3.2. The contingent claim $U(t,S,v)$ satisfies the PDE

$$ \begin{align*} \begin{split} 0&=\partial_tU+rS\partial_SU+(\beta(\alpha-v)-\lambda{|v|})\partial_vU-rU+\tfrac{1}{2}v^2S^2\partial^2_SU\\ &\quad+\rho\gamma_2|v|S\partial^2_{Sv}U+\big(\tfrac{1}{2}\gamma_2^2+\gamma_1\Phi\big)\partial^2_vU, \end{split} \end{align*} $$

where r denotes risk-free interest rate and

$$ \begin{align*} \Phi=\Phi(t,u)=\bigg(\gamma_1Z_H(t,u)-\beta\int^t_0\gamma_1 e^{\beta(s-t)}Z_H(s,u)\,ds\bigg)\mathbf{1}_{\{u<t\}}. \end{align*} $$

Proof. Here, we establish a replicating portfolio. It consists of the underlying asset and a security with price function $R(t,S_t,v_t)$ . Let $V_t=U(t,S_t,v_t)$ . Then the portfolio is given by

$$ \begin{align*} V_t=\Delta_1\,S_t+\Delta_2\,R_t. \end{align*} $$

The SDE of $V_t$ can be written as

$$ \begin{align*} \begin{split} dV_t&=\,\Delta_1\,dS_t+\Delta_2\,dR_t\\ &=\,\Delta_1\,dS_t+\Delta_2\,dR_t+r(V_t-\Delta_2\,R_t)\,dt-r\Delta_1\,S_t\,dt, \end{split} \end{align*} $$

where $\Delta _1$ and $\Delta _2$ are the number of shares of the underlying asset and the security. By It $\hat {\textrm {o}}$ ’s formula, Theorem 2.1 and Lemma 3.1, we have

$$ \begin{align*} \begin{split} dU=\,&(\partial_tU+{\mu}S_t\partial_SU+\beta(\alpha-v_t)\partial_vU+\tfrac{1}{2}v^2_tS^2_t\partial^2_SU+{\gamma_1}\mathbb{D}^H_tv_t\partial^2_vU\\ &+\tfrac{1}{2}\gamma_2^2\partial^2_vU+\rho\gamma_2|v_t|S_t\partial^2_{Sv}U)\,dt+{|v_t|}S_t\partial_SUdB^S_t+\gamma_1\partial_vUdB^H_t+\gamma_2\partial_vUdB^v_t\\ =\,&(\partial_tU+{\mu}S_t\partial_SU+\beta(\alpha-v_t)\partial_vU+\tfrac{1}{2}v^2_tS^2_t\partial^2_SU+\tfrac{1}{2}\gamma_2^2\partial^2_vU\\ &+\gamma_1\Phi\partial^2_vU+\rho\gamma_2|v_t|S_t\partial^2_{Sv}U)\,dt+{|v_t|}S_t\partial_SUdB^S_t+\gamma_1\partial_vUdB^H_t+\gamma_2\partial_vUdB^v_t\\ =\,&U^*dt+{|v_t|}S_t\partial_SUdB^S_t+\gamma_1\partial_vUdB^H_t+\gamma_2\partial_vUdB^v_t, \end{split} \end{align*} $$

where

$$ \begin{align*} U^*=\partial_tU+{\mu}S_t\partial_SU+\beta(\alpha-v_t)\partial_vU+\tfrac{1}{2}v^2_tS^2_t\partial^2_SU+(\tfrac{1}{2}\gamma_2^2+\gamma_1\Phi)\partial^2_vU+\rho\gamma_2|v_t|S_t\partial^2_{Sv}U. \end{align*} $$

On the other hand, we have $R=R(t,S_t,v_t)$ and

$$ \begin{align*} dR_t=R^*dt+{|v_t|}S_t\partial_SRdB^S_t+\gamma_1\partial_vRdB^H_t+\gamma_2\partial_vRdB^v_t, \end{align*} $$

where $R^*$ has the same form as $U^*$ .

By Itô’s formula,

$$ \begin{align*} \begin{split} dV_t&=(\Delta_1\,{\mu}S_t+r(V_t-\Delta_2\,R_t)-r\Delta_1\,S_t+\Delta_2\,R^*)\,dt+(\Delta_1\,{|v_t|}S_t+\Delta_2\,{|v_t|}S_t\partial_SR)\,dB^S_t\\ &\quad +\Delta_2\,\gamma_1\partial_vRd{B}^H_t+\Delta_2\,\gamma_2\partial_vRd{B}^v_t. \end{split} \end{align*} $$

We now compare $dU_t$ and $dV_t$ . Let the coefficients of $dB^S_t$ and $d{B}^H_t$ be equal. We have

$$ \begin{align*} \begin{split} \Delta_1{|v_t|}S_t+\Delta_2{|v_t|}S_t\partial_SR&=\ {|v_t|}S_t\partial_SU,\\ \Delta_2\gamma_1\partial_vR&=\ \gamma_1\partial_vU. \end{split} \end{align*} $$

Then we get

$$ \begin{align*} \left \{ \begin{aligned} &\Delta_1=\partial_SU-\frac{\partial_vU\partial_SR}{\partial_vR},\\ &\Delta_2=\frac{\partial_vU}{\partial_vR}. \end{aligned} \right. \end{align*} $$

Substituting $\Delta _1$ and $\Delta _2$ into the coefficients of $dt$ yields

$$ \begin{align*} U^*=({\mu}S_t-rS_t)\bigg(\partial_SU-\frac{\partial_vU\partial_SR}{\partial_vR}\bigg)+r\bigg(V_t-R_t\frac{{\partial_vU}}{{\partial_vR}}\bigg)+R^*\frac{{\partial_vU}}{{\partial_vR}}. \end{align*} $$

Since $V_t=U_t$ , there is

(3.3) $$ \begin{align} \frac{1}{\partial_vU}(U^*-{\mu}S_t\partial_SU+rS_t\partial_SU-rU_t)=\frac{1}{\partial_vR}(R^*-{\mu}S_t\partial_SR+rS\partial_SR-rR_t). \end{align} $$

Note that the right-hand side depends only on R, while the left-hand side depends only on U. Here R is arbitrary. So both sides must be $\lambda (t,S_t,v_t)$ , which is called the market price of volatility risk. By the previous discussion, we write $\lambda $ as $\lambda (t,S_t,v_t)={\lambda }|v_t|. $ Let $\lambda (t,S_t,v_t)$ be the right-hand side of equation (3.3). Then the proof is complete.

3.2. The expression for the option price

In this section we consider the price of the European call option. Suppose that $\mathfrak {X}$ is a European call option with strike price K and maturity T. The option price can be written as

$$ \begin{align*} U(t,S_t,v_t)=e^{-r(T-t)}E[\Psi(S_T)|\mathcal{F}_t], \quad\text{where}\ \Psi(S_T)=(S_T-K)^+. \end{align*} $$

By Theorem 3.2, we can get the PDE of European call options. To solve the PDE, we first simplify it through a series of transformations.

Let $\tau =T-t$ . The PDE is given by

(3.4) $$ \begin{align} \left \{ \begin{aligned} &\partial_{\tau}U-rS\partial_SU-\tilde{\kappa}(\tilde{\theta}-v)\partial_vU+rU-\tfrac{1}{2}v^2 S^2\partial^2_SU \\ &\quad-\big(\tfrac{1}{2}\gamma_2^2+\gamma_1\Phi\big)\partial^2_vU-\rho\gamma_2|v|S\partial^2_{Sv}U=0, \\ &U(0,S,v)=(S-K)^+, \end{aligned} \right. \end{align} $$

where $\tilde {\kappa }=\beta +\lambda I$ , $\tilde {\theta }=\beta \alpha /(\beta +\lambda I)$ and

$$ \begin{align*} I=\left \{ \begin{aligned} 1\quad &v>0,\\ -1\quad &v\leq 0. \end{aligned} \right. \end{align*} $$

Let $Q=e^{r\tau }S$ and $\hat {U}(\tau ,Q,v)=U(\tau ,S,v)e^{r\tau }.$ Then we have $U(\tau ,S,v)=e^{-r\tau }\hat {U}(\tau ,Q,v)$ .

Through the above transformation, we rewrite equations (3.4) as follows:

(3.5) $$ \begin{align} \left \{ \begin{aligned} &\partial_{\tau}\hat{U}-\tilde{\kappa}(\tilde{\theta}-v)\partial_v\hat{U}-\tfrac{1}{2}v^2Q^2\partial^2_Q\hat{U}-\big(\tfrac{1}{2}\gamma_2^2+\gamma_1\Phi\big)\partial^2_v\hat{U}-\rho\gamma_2Q|v|\partial_{Qv}\hat{U}=0,\\ &\hat{U}(0,Q,v)=(Q-K)^+. \end{aligned} \right. \end{align} $$

Let

(3.6) $$ \begin{align} X=\ln(Q/K),\quad \tilde{U}(\tau,X,v)=\frac{e^{-X/2}}{K}\hat{U}(\tau,Q,v). \end{align} $$

Substituting equations (3.6) into equations (3.5), we obtain

(3.7) $$ \begin{align} \left \{ \begin{aligned} &\partial_{\tau}\tilde{U}-\tilde{\kappa}(\tilde{\theta}-v)\partial_v\tilde{U}-\tfrac{1}{2}v^2\partial^2_X\tilde{U}-\big(\tfrac{1}{2}\gamma_2^2+\gamma_1\Phi\big)\partial^2_v\tilde{U}-\rho\gamma_2|v|\partial_{Xv}\tilde{U}+\tfrac{1}{8}v^2\tilde{U}=0,\\ &\tilde{U}(0,X,v)=(e^{X/2}-e^{-X/2})^+. \end{aligned} \right. \end{align} $$

Next, we give the solution of the PDE (3.7).

Proposition 3.3. Let $a=2\gamma _2^2+4\gamma _1\Phi $ , $b=-2\kappa +2ik\rho \gamma _2,$

$$ \begin{align*} &\quad\xi(k)=-\dfrac{1}{2}k^2-\dfrac{1}{8}, \quad \psi(\tau,k)=\sqrt{\frac{b^2-4a\xi(k)}{4a^2}},\\ &\tilde{\Phi}=\displaystyle\int^\tau_02\gamma_2^2+4\gamma_1\Phi(s,u)\,ds, h(\tau,k)=e^{2\psi(\tau,k)\tilde{\Phi}}\frac{{b}/{2a}-\psi(\tau,k)}{{b}/{2a}+\psi(\tau,k)}. \end{align*} $$

Then the solution of equations (3.7) is

(3.8) $$ \begin{align} \tilde{U}=\frac{1}{2\pi}\int^{+\infty}_{-\infty}\int^{+\infty}_{-\infty}\exp\bigg{(}ik(X-X')+A(\tau,k)+B(\tau,k)|v|+C(\tau,k)v^2\bigg{)}\varphi(X')\,dk\,dX', \end{align} $$

where i is the imaginary unit, $\varphi (X')=(e^{X'/2}-e^{-X'/2})^+$ and

(3.9) $$ \begin{align} \left \{ \begin{aligned} A(\tau,k)&=\int^\tau_0\bigg{[}(\gamma_2^2+2\gamma_1\Phi) C(s,k) +2\kappa^2\theta^2\int^s_0 \exp\bigg{(}(\kappa-ik\rho\gamma_2)(u-\tau)\\ &\quad +2\gamma_2^2\int^\tau_u C(\mu,k)d\mu+4\gamma_1\int^\tau_s\Phi C(\mu,k)d\mu\bigg{)}C(u,k)\,du \\ &\quad +2\gamma_2^2\kappa^2\theta^2\bigg{(}\int^s_{0}\exp\bigg{(}\kappa-ik\rho\gamma_2)(u-\tau)+2\gamma_2^2\int^\tau_u C(\nu,k)\,d\nu\\ &\quad+4\gamma_1\int^\tau_s\Phi C(\nu,k)\,d\nu\bigg{)} C(u,k)\,du\bigg{)}^2 +4\gamma_1\kappa^2\theta^2\bigg{(}\int^s_0\Phi \exp\bigg{(}(\kappa-ik\rho\gamma_2)(u-\tau)\\ &\quad+2\gamma_2^2\int^\tau_u C(\nu,k)\,d\nu C(u,k)\,du+4\gamma_1\int^\tau_s\Phi C(\nu,k)\,d\nu\bigg{)} C(u,k)\,du\bigg{)}^2\bigg{]}\,ds,\\ B(\tau,k)&=2\kappa\theta\int^\tau_0\exp\bigg{(}(\kappa-ik\rho\gamma_2)(s-\tau)+2\gamma_2^2\int^\tau_s C(u,k)\,du \\ &\quad +4\gamma_1\int^\tau_s\Phi C(u,k)\,du\bigg{)} C(s,k)\,ds,\\ C(\tau,k)&=\psi(k)\frac{1+h(\tau,k)}{1-h(\tau,k)}-\frac{b}{2a}.\\\\[-24pt] \end{aligned} \right. \end{align} $$

Proof. Define Green’s function as

$$ \begin{align*} p(t,X,v,t',X',v')=P(X_{t'}=X',|v_{t'}|=v'|X_t=X,|v_t|=v). \end{align*} $$

We observe that $p(t,X,v,t',X',v')$ is the transition probability of $(X_t,|v_t|)$ . Then the pricing formula $\tilde {U}$ can be written as

(3.10) $$ \begin{align} \tilde{U}=\int^\infty_0\int^{+\infty}_{-\infty}p(t,X,v,t',X',v')\varphi(X')\,dX'\,dv', \end{align} $$

where $\varphi (X')=(e^{X'/2}-e^{-X'/2})^+$ . Let Green’s function $p(t,X,v,t',X',v')=p(\tau ,Y,v,v')$ where $\tau =t'-t,{Y}=X'-X$ . Because of the independence of the coefficients and the killing rate being $v^2/8$ with respect to $t,X$ , we can concentrate on the differences $\tau $ and Y. Notice that, actually, v stands for $|v|$ . Substituting equation (3.10) into equations (3.7) yields

(3.11) $$ \begin{align} \left \{ \begin{aligned} &\partial_{\tau}p-\kappa(\theta-v)\partial_vp-\tfrac{1}{2}v^2\partial^2_Yp-\big(\tfrac{1}{2}\gamma_2^2+\gamma_1\Phi\big)\partial^2_vp-\rho\gamma_2vI\partial_{Yv}p+\tfrac{1}{8}v^2p=0,\\ &p(\tau\rightarrow0,Y,v,v')=\delta(Y)\delta(v'-v), \end{aligned} \right. \end{align} $$

where $\kappa =(\beta +\lambda )I,~\theta =\beta \alpha /(\beta +\lambda )$ and $\delta (\cdot )$ is the Dirac delta function with the properties

$$ \begin{align*} \left \{ \begin{aligned} &\delta(x)=0\qquad\qquad\qquad\qquad\ \ \ \text{for all } x\neq 0,\\ &\int^\varepsilon_\varepsilon\delta(x)\,dx=1\qquad\qquad\quad\ \ \,\,\text{for all } \varepsilon>0,\\ &\int^T_0\delta(x-x')f(x)\,dx=f(x')\quad \text{for all } x'\in[0,T]. \end{aligned} \right. \end{align*} $$

Let

$$ \begin{align*} \left \{ \begin{aligned} q(\tau,Y,v)&=\int^\infty_0p(\tau,Y,v,v')\,dv',\\ q(0,Y,v)&=\int^\infty_0\delta(Y)\delta(v'-v)\,dv'\\ &=\delta(Y). \end{aligned} \right. \end{align*} $$

If we get the expression for q, by equation (3.10), we will obtain $\widetilde {U}$ . Suppose that

(3.12) $$ \begin{align} q(\tau,Y,v)=\frac{1}{2\pi}\int^{+\infty}_{-\infty}\exp{(}ikY+A(\tau,k)+B(\tau,k)v+C(\tau,k)v^2{)}\,dk. \end{align} $$

Substitute equation (3.12) into equations (3.11). Note that the exponential function in the above integral still satisfies equations (3.11). We get the following derivatives by letting $g(k)=\exp (ikY+A(\tau ,k)+B(\tau ,k)v+C(\tau ,k)v^2)$ :

$$ \begin{align*} \begin{split} \partial_\tau{g}&=\exp{(}ikY+A(\tau,k)+B(\tau,k)v+C(\tau,k)v^2{)}\cdot(\partial_\tau{A}+\partial_\tau{B}v+\partial_\tau{C}v^2),\\ \partial_vg&=\exp{(}ikY+A(\tau,k)+B(\tau,k)v+C(\tau,k)v^2{)}\cdot(B+2Cv),\\ \partial^2_vg&=\exp{(}ikY+A(\tau,k)+B(\tau,k)v+C(\tau,k)v^2{)}\cdot(2C+(B+2Cv)^2),\\ \partial^2_Yg&=-k^2\exp{(}ikY+A(\tau,k)+B(\tau,k)v+C(\tau,k)v^2{)},\\ \partial_v g&=\exp{(}ikY+A(\tau,k)+B(\tau,k)v+C(\tau,k)v^2{)}\cdot(ik(B+2Cv)). \end{split} \end{align*} $$

Then we rewrite equations (3.11) as

$$ \begin{align*} \begin{split} \partial_\tau{A}&+\partial_\tau{B}v+\partial_\tau{C}v^2-\kappa(\theta-v)(B+2Cv)\\&+\tfrac{1}{2}v^2 k^2 \big(\tfrac{1}{2}\gamma_2^2+\gamma_1\Phi\big)(2C+(B+2Cv)^2)-ik\rho\gamma_2vI(B+2Cv)+\tfrac{1}{8}v^2=0. \end{split} \end{align*} $$

Notice that A, B and C are independent of v. We have

(3.13) $$ \begin{align} \left \{ \begin{aligned} &\partial_\tau{A}-\kappa\theta{B}-\gamma_2^2C-\tfrac{1}{2}\gamma_2^2B^2-2C\gamma_1\Phi-B^2\gamma_1\Phi=0,\\ &\partial_\tau{B}+\kappa{B}-2\kappa\theta{C}-2\gamma_2^2BC-4BC\gamma_1\Phi-ik\rho\gamma_2IB=0,\\ &\partial_\tau{C}+2\kappa{C}+\tfrac{1}{2}k^2-2\gamma_2^2C^2-4C^2\gamma_1\Phi-2ikC\rho\gamma_2I+\tfrac{1}{8}=0. \end{aligned} \right. \end{align} $$

On the other hand, we have

$$ \begin{align*} \begin{split} q(0,Y,v)&=\frac{1}{2\pi}\int^{+\infty}_{-\infty}\exp{(}ikY+A(0,k)+B(0,k)v+C(0,k)v^2{)}\,dk =\delta(Y). \end{split} \end{align*} $$

So, this yields

(3.14) $$ \begin{align} A(0,k)=B(0,k)=C(0,k)=0. \end{align} $$

We solve the equations (3.13) with the boundary conditions (3.14) and get the expressions for A, B and C. This completes the proof of the proposition.

Next, by the transformations, we have

(3.15) $$ \begin{align} U(t,S,v)&=\,U(\tau,S,v) =e^{-r\tau}\hat{U}(\tau,Q,v) =Ke^{-r\tau+X/2}\tilde{U}(\tau,X,v)\nonumber\\ &=\,Ke^{-r\tau+X/2}\frac{1}{2\pi}\int^{+\infty}_{-\infty}\int^{+\infty}_{-\infty}\exp{(}ik(X-X')+A(\tau,k)+B(\tau,k)|v|\nonumber\\ &\quad+C(\tau,k)v^2{)}\varphi(X')\,dk\,dX'\nonumber\\ &=\,\frac{\sqrt{SK}}{2\pi}e^{-r(T-t)/2}\int^{+\infty}_{-\infty}\int^{+\infty}_{-\infty}\exp\bigg{(}ik\,\bigg(\ln\frac{S}{K}+r(T-t)-X'\bigg)+A(T-t,k)\nonumber\\ &\quad+B(T-t,k)|v|+C(T-t,k)v^2\bigg{)}\,\varphi(X')\,dk\,dX'. \end{align} $$

Here, $A(\tau ,k)$ , $B(\tau ,k)$ and $C(\tau ,k)$ can be expressed as in equations (3.9).

4.1. Model comparison

When it comes to Heston stochastic volatility model, we can see that he models variance [Reference Heston21], while the Schöbel–Zhu model [Reference Schöbel and Zhu33] and our model give the different SDEs of volatility. On the other hand, in Heston’s model, the volatility follows an OU process [Reference Stein and Stein32] with zero mean-reversion parameter. Indeed, the volatility process of Heston’s model can be expressed as follows [Reference Heston21]:

(4.1) $$ \begin{align} dv_t=-\widetilde{\beta} v_t\, dt +\widetilde{\gamma}\, dB_t^v. \end{align} $$

Then, by Itô’s formula, the variance process $\Gamma _t=v_t^2$ is

(4.2) $$ \begin{align} d\Gamma_t=\widetilde{\kappa}(\widetilde{\theta}-\Gamma_t)\,dt+\gamma_3\sqrt{\Gamma_t}\,dB_t^v, \end{align} $$

where

(4.3) $$ \begin{align} \widetilde{\beta}=\frac{\widetilde{\kappa}}{2},\quad\widetilde{\gamma}=\frac{\gamma_3}{2},\quad \widetilde{\theta}=\frac{\widetilde{\gamma}^2}{\widetilde{\kappa}}. \end{align} $$

We know that the mean-reversion parameter gives the level of volatility in the long run. So, from equation (4.1), $v_t$ does not seem very reasonable, because the mean-reversion level equals zero. However, according to the discussion in [Reference Schöbel and Zhu33], we know that the two processes (4.1) and (4.2) are not always consistent. For some values of $\widetilde {\kappa },~\widetilde {\theta },~\gamma _3$ , the volatility process cannot be derived from (4.2), because the values of the parameters do not satisfy equations (4.3). When we set the parameters as

(4.4) $$ \begin{align} \beta=\frac{\widetilde{\kappa}}{2},\quad\gamma_1=0,\quad\alpha=0,\quad\gamma_2=\frac{\gamma_3}{2},\quad\widetilde{\theta}=\frac{\gamma_2^2}{\widetilde{\kappa}}, \end{align} $$

our model degenerates to Heston’s model.

In our model and the Schöbel–Zhu model, both volatility and squared volatility perform mean-reversion, which is as an excellent property of our model (the derivation of Schöbel–Zhu model is shown in [Reference Schöbel and Zhu33]). For our model, by Itô’s formula with respect to fBm, we can obtain the SDE of the variance process as follows:

(4.5) $$ \begin{align} \begin{array}{rcl} d\Gamma_t=(2\beta\alpha\sqrt{\Gamma_t}-2\beta\Gamma_t+\gamma^2_2+4H\gamma_1^2t^{2H-1})\,dt+2\gamma_1\sqrt{\Gamma_t}d B_t^H +2\gamma_2\sqrt{\Gamma_t}d B_t^v. \end{array} \end{align} $$

By equation (4.5), the variance process is a fractional mean-reverting double square-root process with the additional drift term $2\beta \alpha \sqrt {\Gamma _t}+4H\gamma _1^2t^{2H-1}$ .

When compared with the Schöbel–Zhu model, we add fBm to the volatility. The stochastic volatility model with fBm is more realistic than the model with sBm, because a large number of empirical studies have shown that the volatility process has a long-term or short-term correlation and the correlation function of the volatility process decays to the fractional power of time deviation [Reference Comte and Renault7, Reference Comte, Coutin and Renault8, Reference Garnier and Sølna17].

4.2. Analysis of our model by numerical calculation

We now analyse our model through numerical calculation. For comparison, we need to choose suitable benchmarks. The first benchmark we choose is the BS option price with the expected average variance as input. The reason for the use of expected average variance is that the option price can be approximated very accurately by applying the variance in the BS formula (see [Reference Ball and Roma1, Reference Hull and White23]). Here, we suppose that the volatility process is a fractional OU process as shown in equation (3.1) and denote the BS prices calculated in this way by I-BS. Without going through the lengthy derivation procedures, the expected average variance can be written as

$$ \begin{align*} E\bigg[\frac{1}{T-t}&\int_t^Tv^2(u)\,du\bigg]\nonumber\\ \quad&=\,\alpha^2+\frac{\gamma_2^2}{2\beta}+\frac{2\alpha}{T-t}(v_0-\alpha)e^{-\beta(T-t)} +\frac{(\alpha-v_0)^2-\gamma_2^2/2\beta}{T-t}e^{-2\beta(T-t)}\nonumber\\ &\quad +\frac{1}{T-t}\int_t^T\gamma_1^2e^{-2\beta s}\int_0^s\int_0^s \{H(2H-1)e^{\beta u}e^{\beta v}|u-v|^{2H-2}\}du\,dv\,ds\nonumber\\ &\quad +\int_t^T\gamma_1^2e^{-2\beta s}\int_0^s\int_0^s E[\mathbb{D}_u^H e^{\beta v}\mathbb{D}_v^H e^{\beta u}]\,du \,dv \,d s. \end{align*} $$

The second benchmark is the BS prices with the target volatility as shown in [Reference Stein and Stein32]. By using the target volatility as constant volatility in the BS model, we can obtain the option prices referred to the case as II-BS [Reference Stein and Stein32]. Next, we choose the option prices under the Schöbel–Zhu model (III-SZ) as the third benchmark. Finally, we write our model as IV-OM.

By using Proposition 3.3, we obtain the closed-form solution of equations (3.7). Through the derivations of equation (3.15), we can get the European call option price. Next, we will give the numerical evaluation of the prices. Note that, in all the numerical experiments, we set $H=3/4$ .

First, we will show the impact of model correlation $\rho $ on the option prices (Table 1 and Figure 1). To compare with other models, we choose the parameter values suggested in [Reference Stein and Stein32, Reference Schöbel and Zhu33].

Table 1 Option prices with different $\rho $ .

Figure 1 Comparison chart of option prices with different correlation coefficients.

In Table 1 we let the correlation range from $-1$ to $1$ . We observe that options with different moneyness can be influenced in different ways. While the values of at-the-money options ( $K=100$ ) and in-the-money (ITM) options ( $K=90,95$ ) decrease with increasing correlation $\rho $ , out-of-the-money (OTM) options ( $K=105,110,115,120$ ) increase in value. This is consistent with the results of Schöbel and Zhu [Reference Schöbel and Zhu33], and the comparison is shown in Figure 1. On the other hand, from panels (a) and (b) of Table 1, we find that whether $\gamma _1$ and $\gamma _2$ are the same or not, ITM (OTM) option prices in our model are greater (less) than I-BS option prices with a negative correlation $\rho $ . This can be used to explain volatility skews. Indeed, the implied volatility with moneyness is monotonically downward sloping, displayed as a skew. This means that the ITM (OTM) option prices are undervalued (overvalued) by the BS formula. In addition, through a simple calculation, we obtain the relative changes of option prices in panel (a). With different $\rho $ , the relative change rates of option prices are $4.09\%$ , $4.39\%$ , $4.65\%$ , $8.67\%$ , $14.77\%$ , $21.14\%$ , $44.78\%$ when $K=90,95,100,105,110,115,120$ , respectively. These results show that the option prices with different moneyness have different sensitivity to $\rho $ . The sensitivity of ITM option prices to $\rho $ is less than that of OTM option prices. Finally, we observe that the price differences between I-BS and our model values with $\rho =0$ are the smallest in most cases. Therefore, I-BS is not a suitable estimation method for the case of nonzero correlation. In panel (a), the prices of I-BS and our model values ( $\rho =0$ ) are all slightly higher than that of II-BS.

For Figure 1, we compare the option prices of our model and that of the Schöbel–Zhu model with the same parameters besides $\gamma _1=0$ . At this point, our model degenerates to the Schöbel–Zhu model. To replicate the results of the Schöbel–Zhu model, the model parameters are set to $\beta =4, ~\alpha =0.2, ~\gamma _1=0,~\gamma _2=0.1, T-t=0.5, ~S=100, ~r=0.0953, ~v=0.2$ , as they should result in the same options as those presented in Table 1 of [Reference Schöbel and Zhu33]. From Figure 1, we can see that the two models match well in this case and this also verifies our results.

Secondly, we analyse the option prices with different time to maturity ( $T-t$ ), and the results are shown in Table 2. In both panels (a) and (b), as time to maturity increases, the option prices with the same strike price increase whatever the moneyness is. Through simple calculation, we find that the sensitivity of option prices with smaller $T-t$ to strike price K is larger than that of option prices with larger $T-t$ .

Table 2 Option prices with different $T-t$ .

Thirdly, we examine how the option prices behave with different mean-reversion level $\alpha $ shown in Table 3 and Figure 2. Notice that, in Figure 2, we choose the same parameters as in Table 3. Through observation, we find that the option prices increase as the mean-reversion level $\alpha $ goes up regardless of the correlation coefficient. From Table 3 and Figure 2, we see that the option prices are sensitive to $\alpha $ . As shown above, the parameter $\alpha $ indicates the long-run level of volatility. So, the volatility values can influence option prices in the long run. Through Figure 2, we see that different $\alpha $ cannot affect the differences of option prices, when we change the strike price K. So, we infer that there is no obvious correlation between the mean-reversion level and the strike price.

Table 3 Option prices with different $\alpha $ .

Figure 2 Comparison chart of option prices with different mean-reversion levels.

4.3. Implied volatility

In this section we analyse the implied volatility through numerical calculation. We know that the limitations of the BS model can be shown by implied volatility skew or smile. We use S $\&$ P 500 stock index data to simulate the option prices and get the values of implied volatilities. The two typical features of the implied volatilities from stock index options are that they are higher than the historical volatilities of the index (see [Reference Fouque, Papanicolaou, Sircar and Sølna16]), and the skews are steeper for shorter maturities as shown in Figure 3.

Figure 3 S $\&$ P 500 implied volatilities.

In Figure 3 we illustrate several sets of implied volatilities with different maturities. Each strand corresponds to a given maturity as a function of moneyness $K/S$ from S $\&$ P 500 index options on 1 June 2007 [Reference Fouque, Papanicolaou, Sircar and Sølna16]. As discussed above, the implied volatility changes more obviously for shorter maturities when we change moneyness. Thus, next we analyse implied volatilities skews for short maturities.

We compare implied volatilities corresponding to our model and Heston’s stochastic volatility model [Reference Heston21] with realistic implied volatilities in Figure 4. The results of implied volatilities with times to maturity 1 month, 3 months and 6 months are shown in panels (a), (b) and (c) of Figure 4, respectively. Through long derivation and calibration, we found the appropriate parameters to simulate the real implied volatilities. For the case of the time to maturity equal to 1 month, the parameters are set to be $v=0.2512,~\alpha =0.2431,~\beta =7.5682,~\gamma _1=\gamma _2=0.2784,~T-t=1/12 \,\text {(year) and} \rho =-0.1125$ . Moreover, for Heston’s model, we set parameters according to equations (4.2)–(4.4). So, with the same $T-t$ and $\rho $ , other parameters should be $v=0.2521,~\widetilde {\kappa }=15.1256,~\gamma _3=0.5568$ and $\widetilde {\theta }=0.0051$ . Through a similar process, we obtain the values of parameters in the other two cases. When $T-t=0.25$ (year), the parameters are $v=0.1523,~\alpha =0.1892,~\beta =3.7268,~\gamma _1=\gamma _2=0.1124,~T-t=0.25$ (year), $~\rho =-0.1328,~\widetilde {\kappa }=7.4536,~\gamma _3=0.2248$ and $\widetilde {\theta }=0.001\,695$ . For the case of 6 months, we set the parameters to be $v=0.1328,~\alpha =0.1428,~\beta =2.4653,~\gamma _1=\gamma _2=0.0928,~T-t=0.5$ (year), $~\rho =-0.1243,~\widetilde {\kappa }=4.9306,~\gamma _3=0.1856$ and $\widetilde {\theta }=0.001\,747$ . We use the sum of squares of the difference between the simulated implied volatilities and the realistic implied volatilities to measure the results of the two models. Let the implied volatility corresponding to our model and Heston’s model be $I^{(1)}$ and $I^{(2)}$ , respectively. The realistic implied volatility is $I_R$ . Therefore, the sum of squared errors can be expressed as

$$ \begin{align*} E_i=\sum\limits_{j=1}^{n}(I_j^{(i)}-I_R)^2,\quad i=1,2, \end{align*} $$

where n represents the number of simulated implied volatilities; $E_1$ and $E_2$ are the sum of the squared errors corresponding to the two models.

Figure 4 Implied volatilities with different maturities.

Table 4 The sum of squared errors under two models.

The simulation error results are shown in Table 4. On the whole, the errors of our model are smaller than that of Heston’s model. Therefore, from the perspective of implied volatilities, this can be seen as an advantage of our model. On the other hand, with increasing the time to maturity, both $E_1$ and $E_2$ decrease. So, the simulation result in the case of longer time to maturity is better than that for the shorter one.