Option Pricing under Double Heston Model with Approximative Fractional Stochastic Volatility

Abstract

We establish double Heston model with approximative fractional stochastic volatility in this article. Since approximative fractional Brownian motion is a better choice compared with Brownian motion in financial studies, we introduce it to double Heston model by modeling the dynamics of the stock price and one factor of the variance with approximative fractional process and it is our contribution to the article. We use the technique of Radon–Nikodym derivative to obtain the semianalytical pricing formula for the call options and derive the characteristic functions. We do the calibration to estimate the parameters. The calibration demonstrates that the model provides the best performance among the three models. The numerical result demonstrates that the model has better performance than the double Heston model in fitting with the market implied volatilities for different maturities. The model has a better fit to the market implied volatilities for long-term options than for short-term options. We also examine the impact of the positive approximation factor and the long-memory parameter on the call option prices.

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

We propose a double Heston model with approximative fractional stochastic volatility in this article. Though the Black and Scholes [1] model has been a benchmark in European option pricing studies, the existence of the volatility smile made the stochastic volatility model more proper for application. The Hull and White [2] model, the Scott [3] model, the Wiggins [4] model, the Stein and Stein [5] model, and the Lewis [6] model are developed, but most of these stochastic volatility models are not proper for application. The Heston [7] model which is determined by a Cox–Ingersoll–Ross (CIR) [8] process displays the property of mean reversion effect and the character of stationarity. It is the pioneering work in the history of stochastic volatility model studies since Heston technically derived the closed-form formula for the option prices by using characteristic functions. Bakshi et al. [9] empirically showed that stochastic volatility model has better performance than the Black–Scholes [1] model. However, the one-factor stochastic volatility model has some restrictions on exploring the correlation between the volatility and the smile. For instance, one-factor stochastic volatility model cannot explain some specific phenomena that the figure of the smile is unconnected with the volatility. Christoffersen et al. [10] developed a double Heston model for exploring the correlation between the volatility and the smile. Their study showed that a two-factor model can explain the correlation between the volatility and the figure of the smile and has more flexibleness to present the volatility term structure. One factor having greater mean-reverting level explains the correlation between short-period returns and the square of volatility. The other one having smaller mean-reverting level explains the correlation between long-period returns and the square of volatility.

The financial market owns the properties of self-similarity and long-range dependence. One should consider it in financial studies. Fractional Brownian motion is centered Gaussian process and it also owns the two properties. Therefore, compared with Brownian motion, it is a better financial tool. Hu and Øksendal [11] developed Wick products to analyze it because it is not a Markov process, or semimartingale, and therefore, the traditional Itô theory is not appropriate for the application of it. Its biggest disadvantage is that it cannot be economically explained very well [12]. If some specific conditions hold, Brownian motion can be expressed equivalently by using the mixed fractional Brownian motion [13], and therefore, one can use the latter for financial application [14–16]. Approximative fractional Brownian motion can also be used to replace fractional Brownian motion [17]. Two terms multiplied by a function can be used to express fractional Brownian motion; one can only concern the second term since it has long-range dependence and the first term has continuous trajectories [17]. Approximative fractional Brownian motion is the approximation of the second term and it is a semimartingale [17]. Thao [17] has proved that when the dynamics of the stock price is driven by approximative fractional Brownian motion, it has no arbitrage opportunity. Several papers have used it as an efficient tool in financial modeling [18, 19]. Several papers presented empirical studies on option pricing models with approximative fractional Brownian motion. Mrázek et al. [20] developed an approximative fractional stochastic volatility model for option pricing application. They applied the Lewis formula [6] for obtaining the semi-closed-form pricing formula and the fundamental transform and used the optimization approach for calibration. They considered two errors and five weights for calibration. Pospíšil and Sobotka [21] also developed an approximative fractional stochastic volatility model. They used the Heston approach [7] for obtaining the semi-closed-form pricing formula and the characteristic functions and used the optimization approach for calibration. They considered three errors and three weights for calibration. Several authors also developed an option pricing model with approximative fractional Brownian motion from a different point of view. Kang et al. [22] developed an approximative fractional stochastic volatility model to study foreign exchange (FX) options. They derived the pricing formula and used some numerical approach to examine the impact of the long-memory parameter and the positive approximation factor on the option prices.

There are two typical approaches to derive the formulae for the characteristic functions. The first approach is developed by Duffie et al. [23]. They expended the Heston [7] model on the multifactor models which are called the affine jump-diffusion (AJD) models. Using the Feynman–Kac theorem, the closed form of characteristic functions can be derived. It is generally used to derive the characteristic functions for obtaining the pricing formula for the options. The advantage of the approach is that the solutions for the characteristic functions can be obtained when some specific conditions are met. The second one is introduced by Schöbel and Zhu [24] and developed by Zhu [25]. Schöbel and Zhu [24] used the expectation approach to derive the characteristic functions, and Zhu [25] extended it on the multifactor models and used the modular approach to derive the characteristic functions. Specifically, the formulae for the characteristic functions can be expressed as a structure consisting of several modules, and the formulae for the modules are derived by the expectation approach. The advantage of this approach is that when the option pricing models are determined by several stochastic factors, one can obtain the characteristic functions by integrating the corresponding parts of the formulae for the characteristic functions to be a unified structure. Ahlip et al. [26] proposed a double Heston hybrid model with jumps. They used expectation approach and martingale approach to obtain the pricing formula for the options. Meanwhile, another approach can also be used in several papers. For instance, He and Zhu [27] developed an option pricing model; they expressed the pricing formula for the options with the forward characteristic function. To obtain the formula for the options and corresponding characteristic function, they changed the dynamics of their model to the forward measure.

The loss function approach is generally used for the estimation of the parameters and there are two classes of the loss functions [28]. The first class is minimizing the error between the quoted and model prices. There are two approaches generally used in the first class. The first one is to minimize the absolute value of the mean-squared error (MSE) between the quoted and model prices. Since this approach puts too much weight on long-term and in-the-money options, some authors prefer the second approach for estimating parameters [29]. The second one is to minimize the relative value of the mean-squared error (RVMSE) between the quoted and model prices. The second class is minimizing the error between the quoted and model implied volatilities. If one considers the denominator to be the weight attached to the loss function, there are several choices one can decide to be the weight. Christoffersen et al. [10] developed a loss function that they used the Black–Scholes sensitivity as the weight. The difference between the market prices of the bid and ask also can be taken as the weight [21]. Though there is no concern that which approach is better, the approach used for estimation and evaluation should be the same [29].

We develop double Heston model with approximative fractional stochastic volatility in this article. It is an extensive study based on the study of Ahlip et al. [26] and Kang et al. [22]. Our contribution is that we introduce approximative fractional Brownian motion to double Heston model by modeling the dynamics of the stock price and one factor of the variance with approximative fractional process and estimate the parameters under this model. We do some pioneering work to study this innovative model. We derive the semianalytical pricing formula. We do calibration to estimate the parameters and examine the impact of the long-memory parameter and the positive approximation factor on the call option prices.

We form the article as follows. We develop the model and then derive the semianalytical pricing formula in Section 2. We derive the characteristic functions in Section 3. We estimate the parameters and give some numerical study in Section 4. The conclusion is presented in Section 5.

2. The Double Heston Model with Approximative Fractional Stochastic Volatility

Let $\mathrm{\Omega },\mathrm{ℱ},\mathrm{\mathbb{P}}$ be a probability space; $\mathbb{F}={{\mathrm{ℱ}}_t}_{t\in 0,T}$ is the filtration and $\mathrm{\mathbb{P}}$ is the risk-neutral probability measure. The stock price process $S={S_t}_{t\in 0,T}$ and the variance processes $v={v_t}_{t\in 0,T}$ and $\widehat{v}={{\widehat{v}}_t}_{t\in 0,T}$ are given by the following dynamic system:$$\begin{array}{}\text{(1)}& \begin{array}{c}\frac{\text{d}{S}_t}{S_t}=r\text{d}t+{\sigma }_s\text{d}{W}_t^S+\xi \sqrt{v_t}\text{d}{\tilde{B}}_t^{\epsilon ,H}+\sqrt{{\widehat{v}}_t}\text{d}{\widehat{W}}_t^S,\\ \text{d}{v}_t=\theta -\kappa v_t\text{d}t+{\sigma }_v\sqrt{v_t}\text{d}{\tilde{B}}_t^{\epsilon ,H},\\ \text{d}{v}_t=\widehat{\theta }-\widehat{\kappa }{\widehat{v}}_t\text{d}t+{\sigma }_{\widehat{v}}\sqrt{{\widehat{v}}_t}\text{d}{\widehat{W}}_t^v.\end{array}\end{array}$$

First, we give some preliminaries about fractional Brownian motion ${B_t^H}_{t\ge 0}$ and the approximation of it. It is a Gaussian process including Hurst index $H\in \mathrm{0,1}$. It has zero mean and the covariance equation $\mathbb{E}{B}_t^H{B}_s^H=1/2t^{2H}+s^{2H}-{t-s}^{2H}$. For $H=1/2$, it is a standard Brownian motion. We usually assume $H\in 1/2,1$ for financial application since it displays long-range dependence in that case [19]. It can be decomposed as$$\begin{array}{}\text{(2)}& B_t^H=\frac{1}{\Gamma H+1/2}{Z}_t+{\displaystyle {\int }_0^t{t-s}^{H-1/2}\text{d}{W}_s},\end{array}$$where $Z_t={\displaystyle {\int }_{-\infty }^0{t-s}^{H-1/2}-{-s}^{H-1/2}\text{d}{W}_s}$, $W_t$ denotes a standard Brownian motion, and $\mathrm{\Gamma }$ denotes the gamma function. Since $Z_t$ has continuous trajectories, we can just consider the second term in the bracket which is ${\tilde{B}}_t={\displaystyle {\int }_0^t{t-s}^{H-1/2}\text{d}{W}_s}$ and it has long-range dependence property [17, 19–21]. ${\tilde{B}}_t$ can be approximated by approximative fractional Brownian motion ${\tilde{B}}_t^{\epsilon ,H}$ which is given by [17]$$\begin{array}{}\text{(3)}& {\tilde{B}}_t^{\epsilon ,H}={\displaystyle {\int }_0^t{t-s+\epsilon }^{H-1/2}\text{d}{W}_s}.\end{array}$$

${\tilde{B}}_t^{\epsilon ,H}\stackrel{L^2\mathrm{\Omega }}{⟶}{\tilde{B}}_t$ as $\epsilon ⟶0^{+}$. ${\tilde{B}}_t^{\epsilon ,H}$ is semimartingale [17]. The dynamics of it can be given by [17]$$\begin{array}{}\text{(4)}& \text{d}{\tilde{B}}_t^{\epsilon ,H}=H-\frac{1}{2}{\psi }_t\text{d}t+{\epsilon }^{H-1/2}\text{d}{W}_t^v,\end{array}$$where $H$ represents a long-memory parameter, we consider $H\in 1/2,1$ since we only concern the long-range dependence property [20, 21], $\epsilon$ is an approximation factor, and ${\psi }_t$ is a stochastic process given by$$\begin{array}{}\text{(5)}& {\psi }_t={\displaystyle {\int }_0^t{t-s+\epsilon }^{H-3/2}\text{d}{W}_s^{\psi }}.\end{array}$$

Some assumptions are given by the following statements:

(1) ${\widehat{W}}^S={{\widehat{W}}_t^S}_{t\in 0,T}$ and ${\widehat{W}}^v={{\widehat{W}}_t^v}_{t\in 0,T}$ are correlated Brownian motions that $d{{\widehat{W}}^S,{\widehat{W}}^v}_t=\widehat{\rho }\text{d}t$, $\widehat{\rho }\in -\mathrm{1,1}$

(2) $W^S={W_t^S}_{t\in 0,T}$, $W^v={W_t^v}_{t\in 0,T}$, and $W^{\psi }={W_t^{\psi }}_{t\in 0,T}$ are mutually independent Brownian motions

The call option price formula at time $t\in 0,T$ with strike $K$ under $\mathrm{\mathbb{P}}$ is given by$$\begin{array}{c}\text{(6)}& C_{t}T,K={\mathbb{E}}_t^{\mathrm{\mathbb{P}}}{e}^{-rT-t}{C}_{T}T,K={\mathbb{E}}_t^{\mathrm{\mathbb{P}}}{e}^{-rT-t}{S_T-K}^{+}=e^{-rT-t}{\mathbb{E}}_t^{\mathrm{\mathbb{P}}}{S}_T{1}_{S_T\succ K}-K{e}^{-rT-t}{\mathbb{E}}_t^{\mathrm{\mathbb{P}}}{1}_{S_T\succ K}.\end{array}$$

The forward price for the stock $Ft,T$ at time $t$ can be defined by $Ft,T:=S_t/e^{-rT-t}$, and then, we can have that $S_T=FT,T$, and $C_{t}T,K$ can be expressed as$$\begin{array}{}\text{(7)}& C_{t}T,K=e^{-rT-t}{\mathbb{E}}_t^{\mathrm{\mathbb{P}}}FT,T{1}_{S_T\succ K}-K{e}^{-rT-t}{\mathbb{E}}_t^{\mathrm{\mathbb{P}}}{1}_{S_T\succ K}.\end{array}$$

Lemma 1.

The dynamics of $Ft,T$ under $\mathrm{\mathbb{P}}$ can be given by the following equation:$$\begin{array}{}\text{(8)}& FT,T=Ft,T\exp {\displaystyle {\int }_t^T{\tilde{\sigma }}_{F}u,T\text{d}{\tilde{W}}_u-\frac{1}{2}{\displaystyle {\int }_t^T{{\tilde{\sigma }}_{F}u,T}^2\text{d}u}},\end{array}$$where ${{\tilde{\sigma }}_{F}u,T}_{t\in 0,T}$ and $\tilde{W}={{\tilde{W}}_t}_{t\in 0,T}$ are ${ℝ}^3$-valued processes given by$$\begin{array}{c}\text{(9)}& {\tilde{\sigma }}_{F}u,T={\sigma }_S,\xi {\epsilon }^{H-1/2}\sqrt{v_t},\sqrt{{\widehat{v}}_t},\text{d}{\tilde{W}}_t={\text{d}{W}_t^S+\frac{1}{{\sigma }_S}\xi H-\frac{1}{2}{\psi }_t\sqrt{v_t}\text{d}t,\text{d}{W}_t^v,\text{d}{\widehat{W}}_t^S}^{\ast }.\end{array}$$

Proof.

We can obtain the following equation using Itô’s lemma:$$\begin{array}{}\text{(10)}& \text{d}Ft,T=Ft,T{\sigma }_S\text{d}{W}_t^S+\xi H-\frac{1}{2}{\psi }_t\sqrt{v_t}\text{d}t+\xi {\epsilon }^{H-1/2}\sqrt{v_t}\text{d}{W}_t^v+\sqrt{{\widehat{v}}_t}\text{d}{\widehat{W}}_t^S.\end{array}$$

We can obtain the expression for the dynamics of $Ft,T$ by applying Itô’s lemma, and the proof is complete.

Therefore, we can obtain the following equation:$$\begin{array}{}\text{(11)}& e^{-rT-t}{\mathbb{E}}_t^{\mathrm{\mathbb{P}}}FT,T{1}_{s_T\succ K}=S_t{\mathbb{E}}_t^{\mathrm{\mathbb{P}}}\exp {\displaystyle {\int }_t^T{\tilde{\sigma }}_{F}u,T\text{d}{\tilde{W}}_u-\frac{1}{2}{\displaystyle {\int }_t^T{{\tilde{\sigma }}_{F}u,T}^2\text{d}u}}{1}_{s_T\succ K}.\end{array}$$

Definition 1.

We define the forward measure ${\widehat{\mathrm{\mathbb{P}}}}_T$ which is equivalent to $\mathrm{\mathbb{P}}$ on $\mathrm{\Omega },{\mathrm{ℱ}}_T$:$$\begin{array}{}\text{(12)}& {\widehat{\eta }}_t={\frac{d{\widehat{\mathrm{\mathbb{P}}}}_T}{d\mathrm{\mathbb{P}}}|}_{{\mathrm{ℱ}}_t}=\exp {\displaystyle {\int }_0^t{\tilde{\sigma }}_{F}u,T,\text{d}{\tilde{W}}_u-\frac{1}{2}{\displaystyle {\int }_0^t{{\tilde{\sigma }}_{F}u,T}^2\text{d}u}}.\end{array}$$

Therefore, the formula for the call option price can be expressed as follows:$$\begin{array}{c}\text{(13)}& C_{t}T,K=S_t{\mathbb{E}}_t^{{\widehat{\mathrm{\mathbb{P}}}}_T}{1}_{S_T\succ K}-K{e}^{-rT-t}{\mathbb{E}}_t^{\mathrm{\mathbb{P}}}{1}_{S_T\succ K}=S_t{\widehat{\mathrm{\mathbb{P}}}}_T{S}_T\succ K|S_t,v_t,{\widehat{v}}_t-K{e}^{-rT-t}\mathrm{\mathbb{P}}{S}_T\succ K|S_t,v_t,{\widehat{v}}_t\\ =S_t{\widehat{\mathrm{\mathbb{P}}}}_T{x}_T\succ \ln K|S_t,v_t,{\widehat{v}}_t-K{e}^{-rT-t}\mathrm{\mathbb{P}}{x}_T\succ \ln K|S_t,v_t,{\widehat{v}}_t\\ =S_t{P}_{1}t,S_t,v_t,{\widehat{v}}_t,K-K{e}^{-rT-t}{P}_{2}t,S_t,v_t,{\widehat{v}}_t,K,\end{array}$$where $x_t=\ln Ft,T$.$$\begin{array}{c}\text{(14)}& P_{1}t,S_t,v_t,{\widehat{v}}_t,K={\widehat{\mathrm{\mathbb{P}}}}_T{x}_T\succ \ln K|S_t,v_t,{\widehat{v}}_t=\frac{1}{2}+\frac{1}{\pi }{\displaystyle {\int }_0^{\infty }\mathrm{Re}\frac{e^{-i\varphi \ln K}{f}_1\varphi }{i\varphi }}\text{d}\varphi ,\\ P_{2}t,S_t,v_t,{\widehat{v}}_t,K=\mathrm{\mathbb{P}}{x}_T\succ \ln K|S_t,v_t,{\widehat{v}}_t\\ =\frac{1}{2}+\frac{1}{\pi }{\displaystyle {\int }_0^{\infty }\mathrm{Re}\frac{e^{-i\varphi \ln K}{f}_2\varphi }{i\varphi }}\text{d}\varphi ,\end{array}$$where $f_1\varphi :f_1\varphi ,t,S_t,v_t,{\widehat{v}}_t={\mathbb{E}}_t^{{\widehat{\mathrm{\mathbb{P}}}}_T}\exp i\varphi x_T$ and $f_2\varphi :f_2\varphi ,t,S_t,v_t,{\widehat{v}}_t={\mathbb{E}}_t^{\mathrm{\mathbb{P}}}\exp i\varphi x_T$, and $f_1\varphi$ and $f_1\varphi$ are the characteristic functions.

Therefore, we only need to derive the formulae for the two characteristic functions for obtaining the formula for the call option price.

3. The Derivation of the Characteristic Functions

We derive the formulae for the characteristic functions in this section.

Theorem 1.

Suppose that the stock price has the dynamics given by the model, then the formulae for the characteristic functions $f_1\varphi$ and $f_2\varphi$ can be expressed as follows:$$\begin{array}{c}\text{(15)}& f_1\varphi =c_t\exp \frac{\alpha \beta {\sigma }_S^2\tau }{2}-\frac{\alpha \xi }{{\sigma }_v}{v}_t+\theta \tau -\frac{\alpha \widehat{\rho }}{{\sigma }_{\widehat{v}}}{\widehat{v}}_t+\widehat{\theta }\tau \times {\mathbb{E}}_t^{\mathrm{\mathbb{P}}}\exp -s_1{v}_T-s_2{\displaystyle {\int }_t^T{v}_u\text{d}u}-s_3{\widehat{v}}_T-s_4{\displaystyle {\int }_t^T{\widehat{v}}_u\text{d}u},\\ \text{(16)}& f_2\varphi =c_t\exp -\frac{\beta \gamma {\sigma }_S^2\tau }{2}-\frac{\beta \xi }{{\sigma }_v}{v}_t+\theta \tau -\frac{\beta \widehat{\rho }}{{\sigma }_{\widehat{v}}}{\widehat{v}}_t+\widehat{\theta }\tau \times {\mathbb{E}}_t^{\mathrm{\mathbb{P}}}\exp -q_1{v}_T-q_2{\displaystyle {\int }_t^T{v}_u\text{d}u}-q_3{\widehat{v}}_T-q_4{\displaystyle {\int }_t^T{\widehat{v}}_u\text{d}u},\end{array}$$where $\tau =T-t$, $c_t=\exp i\varphi x_t$, $\alpha =1+i\varphi$, $\beta =i\varphi$, and $\gamma =1-i\varphi$. $s_1$, $s_2$, $s_3$, and $s_4$ are given by $s_1=-\alpha \xi /{\sigma }_v$, $s_2=-\alpha \xi \kappa /{\sigma }_v-\xi {\epsilon }^{2H-1}/2$, $s_3=-\alpha \widehat{\rho }/{\sigma }_{\widehat{v}}$, and $s_4=-{\alpha }^{2}1-{\widehat{\rho }}^2/2+\alpha \widehat{\kappa }\widehat{\rho }/{\sigma }_{\widehat{v}}-\alpha /2$. $q_1$, $q_2$, $q_3$, and $q_4$ are given by $q_1=-\beta \xi /{\sigma }_v$, $q_2=-\beta \xi \kappa /{\sigma }_v-\xi {\epsilon }^{2H-1}/2$, $q_3=-\beta \widehat{\rho }/{\sigma }_{\widehat{v}}$, and $q_4=-{\beta }^{2}1-{\widehat{\rho }}^2/2+\beta \widehat{\kappa }\widehat{\rho }/{\sigma }_{\widehat{v}}-\beta /2$.

Proof.

We can obtain the following equation:$$\begin{array}{c}\text{(17)}& f_1\varphi ={\mathbb{E}}_t^{{\widehat{\mathrm{\mathbb{P}}}}_T}\exp i\varphi x_T={\mathbb{E}}_t^{\mathrm{\mathbb{P}}}\exp i\varphi x_T\exp {\displaystyle {\int }_t^T{\sigma }_S\text{d}{W}_u^S}+{\displaystyle {\int }_t^T\xi H-\frac{1}{2}{\psi }_u\sqrt{v_u}\text{d}u}+{\displaystyle {\int }_t^T\xi {\epsilon }^{H-1/2}\sqrt{v_u}\text{d}{W}_u^v}+{\displaystyle {\int }_t^T\sqrt{{\widehat{v}}_u}\text{d}{\widehat{W}}_u^S}-\frac{1}{2}{\displaystyle {\int }_t^T{\sigma }_S^2+{\xi }^2{\epsilon }^{2H-1}{v}_u+{\widehat{v}}_u\text{d}u}.\end{array}$$

Therefore, we can obtain that$$\begin{array}{c}\text{(18)}& f_1\varphi ={\mathbb{E}}_t^{\mathrm{\mathbb{P}}}\exp i\varphi x_t+{\displaystyle {\int }_t^T{\sigma }_S\text{d}{W}_u^S}+{\displaystyle {\int }_t^T\xi H-\frac{1}{2}{\psi }_u\sqrt{v_u}\text{d}u}+{\displaystyle {\int }_t^T\xi {\epsilon }^{H-1/2}\sqrt{v_u}\text{d}{W}_u^v}+{\displaystyle {\int }_t^T\sqrt{{\widehat{v}}_u}\text{d}{\widehat{W}}_u^S}\times \exp -\frac{i\varphi }{2}{\displaystyle {\int }_t^T{\sigma }_S^2+{\xi }^2{\epsilon }^{2H-1}{v}_u+{\widehat{v}}_u\text{d}u}\\ \times \exp {\displaystyle {\int }_t^T{\sigma }_S\text{d}{W}_u^S}+{\displaystyle {\int }_t^T\xi H-\frac{1}{2}{\psi }_u\sqrt{v_u}\text{d}u}+{\displaystyle {\int }_t^T\xi {\epsilon }^{H-1/2}\sqrt{v_u}\text{d}{W}_u^v}+{\displaystyle {\int }_t^T\sqrt{{\widehat{v}}_u}\text{d}{\widehat{W}}_u^S}\\ \times \exp -\frac{1}{2}{\displaystyle {\int }_t^T{\sigma }_S^2+{\xi }^2{\epsilon }^{2H-1}{v}_u+{\widehat{v}}_u\text{d}u}.\end{array}$$

We insert $c_t=\exp i\varphi x_t$, $\alpha =1+i\varphi$, and $\beta =i\varphi$ into the above equation, and then we can obtain that$$\begin{array}{c}\text{(19)}& f_1\varphi =c_t{\mathbb{E}}_t^{\mathrm{\mathbb{P}}}\exp \alpha {\displaystyle {\int }_t^T{\sigma }_S\text{d}{W}_u^S}+{\displaystyle {\int }_t^T\xi H-\frac{1}{2}{\psi }_u\sqrt{v_u}\text{d}u}+{\displaystyle {\int }_t^T\xi {\epsilon }^{H-1/2}\sqrt{v_u}\text{d}{W}_u^v}+{\displaystyle {\int }_t^T\sqrt{{\widehat{v}}_u}\text{d}{\widehat{W}}_u^S}\times \exp -\frac{\alpha }{2}{\displaystyle {\int }_t^T{\sigma }_S^2+{\xi }^2{\epsilon }^{2H-1}{v}_u+{\widehat{v}}_u\text{d}u}.\end{array}$$

We decompose ${\widehat{W}}^S$ into two parts:$$\begin{array}{}\text{(20)}& {\widehat{W}}_t^S=\widehat{\rho }{\widehat{W}}_t^v+\sqrt{1-{\widehat{\rho }}^2}{\widehat{W}}_t,\end{array}$$where $\widehat{W}={{\widehat{W}}_t}_{t\in 0,T}$ is a Brownian motion under $\mathrm{\mathbb{P}}$ and not correlated with $W^S$, ${\widehat{W}}^S$, $W^v$, ${\widehat{W}}^v$, and $W^{\psi }$. Consequently, we can obtain that$$\begin{array}{c}\text{(21)}& f_1\varphi =c_t{\mathbb{E}}_t^{\mathrm{\mathbb{P}}}\exp \alpha {\displaystyle {\int }_t^T{\sigma }_S\text{d}{W}_u^S}+{\displaystyle {\int }_t^T\xi H-\frac{1}{2}{\psi }_u\sqrt{v_u}\text{d}u}+{\displaystyle {\int }_t^T\xi {\epsilon }^{H-1/2}\sqrt{v_u}\text{d}{W}_u^v}+\widehat{\rho }{\displaystyle {\int }_t^T\sqrt{{\widehat{v}}_u}\text{d}{\widehat{W}}_u^v}+\sqrt{1-{\widehat{\rho }}^2}{\displaystyle {\int }_t^T\sqrt{{\widehat{v}}_u}\text{d}{\widehat{W}}_u}\times \exp -\frac{\alpha }{2}{\displaystyle {\int }_t^T{\sigma }_S^2+{\xi }^2{\epsilon }^{2H-1}{v}_u+{\widehat{v}}_u\text{d}u}.\end{array}$$

We denote $\tau =T-t$ and we can obtain the following equations according to the dynamics of the model:$$\begin{array}{c}\text{(22)}& {\displaystyle {\int }_t^T\sqrt{v_u}\text{d}{W}_u^v=\frac{{\epsilon }^{1/2-H}}{{\sigma }_v}{v}_T-v_t-\theta \tau +\kappa {\displaystyle {\int }_t^T{v}_u\text{d}u}-{\displaystyle {\int }_t^{T}H-\frac{1}{2}}{\sigma }_v{\psi }_u\sqrt{v_u}\text{d}u,}{\displaystyle {\int }_t^T\sqrt{{\widehat{v}}_u}\text{d}{\widehat{W}}_u^v}=\frac{1}{{\sigma }_{\widehat{v}}}{\widehat{v}}_T-{\widehat{v}}_t-\widehat{\theta }\tau +\widehat{\kappa }{\displaystyle {\int }_t^T{\widehat{v}}_u\text{d}u}.\end{array}$$

We insert the above equations into the above expression for $f_1\varphi$ and rearrange it, and then we can obtain that$$\begin{array}{c}\text{(23)}& f_1\varphi =c_t{\mathbb{E}}_t^{\mathrm{\mathbb{P}}}\exp \frac{\alpha \beta {\sigma }_S^2\tau }{2}-\frac{\alpha \xi }{{\sigma }_v}{v}_t+\theta \tau -\frac{\alpha \widehat{\rho }}{{\sigma }_{\widehat{v}}}{\widehat{v}}_t+\widehat{\theta }\tau \times \exp \frac{\alpha \xi }{{\sigma }_v}{v}_T+\alpha \xi \frac{\kappa }{{\sigma }_v}-\frac{\xi {\epsilon }^{2H-1}}{2}{\displaystyle {\int }_t^T{v}_u\text{d}u}+\frac{\alpha \widehat{\rho }}{{\sigma }_{\widehat{v}}}{\widehat{v}}_T+\frac{{\alpha }^{2}1-{\widehat{\rho }}^2}{2}+\frac{\alpha \widehat{\kappa }\widehat{\rho }}{{\sigma }_{\widehat{v}}}-\frac{\alpha }{2}{\displaystyle {\int }_t^T{\widehat{v}}_u\text{d}u}.\end{array}$$

We can rearrange it to obtain that$$\begin{array}{c}\text{(24)}& f_1\varphi =c_t\exp \frac{\alpha \beta {\sigma }_S^2\tau }{2}-\frac{\alpha \xi }{{\sigma }_v}{v}_t+\theta \tau -\frac{\alpha \widehat{\rho }}{{\sigma }_{\widehat{v}}}{\widehat{v}}_t+\widehat{\theta }\tau \times {\mathbb{E}}_t^{\mathrm{\mathbb{P}}}\exp \frac{\alpha \xi }{{\sigma }_v}{v}_T+\alpha \xi \frac{\kappa }{{\sigma }_v}-\frac{\xi {\epsilon }^{2H-1}}{2}{\displaystyle {\int }_t^T{v}_u\text{d}u}+\frac{\alpha \widehat{\rho }}{{\sigma }_{\widehat{v}}}{\widehat{v}}_T+\frac{{\alpha }^{2}1-{\widehat{\rho }}^2}{2}+\frac{\alpha \widehat{\kappa }\widehat{\rho }}{{\sigma }_{\widehat{v}}}-\frac{\alpha }{2}{\displaystyle {\int }_t^T{\widehat{v}}_u\text{d}u}.\end{array}$$

Therefore, we can obtain the following equation:$$\begin{array}{c}\text{(25)}& f_1\varphi =c_t\exp \frac{\alpha \beta {\sigma }_S^2\tau }{2}-\frac{\alpha \xi }{{\sigma }_v}{v}_t+\theta \tau -\frac{\alpha \widehat{\rho }}{{\sigma }_{\widehat{v}}}{\widehat{v}}_t+\widehat{\theta }\tau \times {\mathbb{E}}_t^{\mathrm{\mathbb{P}}}\exp -s_1{v}_T-s_2{\displaystyle {\int }_t^T{v}_u\text{d}u}-s_3{\widehat{v}}_T-s_4{\displaystyle {\int }_t^T{\widehat{v}}_u\text{d}u},\end{array}$$where $s_1$, $s_2$, $s_3$, and $s_4$ are given by equation (15).

We can obtain the formula for $f_2\varphi$ by applying the same method used for deriving the formula for $f_1\varphi$:$$\begin{array}{c}\text{(26)}& f_2\varphi ={\mathbb{E}}_t^{\mathrm{\mathbb{P}}}\exp i\varphi x_T={\mathbb{E}}_t^{\mathrm{\mathbb{P}}}\exp i\varphi x_t+{\displaystyle {\int }_t^T{\sigma }_S\text{d}{W}_u^S}+{\displaystyle {\int }_t^T\xi H-\frac{1}{2}{\psi }_u\sqrt{v_u}\text{d}u}+{\displaystyle {\int }_t^T\xi {\epsilon }^{H-1/2}\sqrt{v_u}\text{d}{W}_u^v}+{\displaystyle {\int }_t^T\sqrt{{\widehat{v}}_u}\text{d}{\widehat{W}}_u^S}\\ \times \exp -\frac{i\varphi }{2}{\displaystyle {\int }_t^T{\sigma }_S^2+{\xi }^2{\epsilon }^{2H-1}{v}_u+{\widehat{v}}_u\text{d}u}.\end{array}$$

Inserting $\beta =i\varphi$ and decomposition of ${\widehat{W}}^S$ into the above equation, then we can obtain that$$\begin{array}{c}\text{(27)}& f_2\varphi =c_t{\mathbb{E}}_t^{\mathrm{\mathbb{P}}}\exp \beta {\displaystyle {\int }_t^T{\sigma }_S\text{d}{W}_u^S}+{\displaystyle {\int }_t^T\xi H-\frac{1}{2}{\psi }_u\sqrt{v_u}\text{d}u}+{\displaystyle {\int }_t^T\xi {\epsilon }^{H-1/2}\sqrt{v_u}\text{d}{W}_u^v}+\widehat{\rho }{\displaystyle {\int }_t^T\sqrt{{\widehat{v}}_u}\text{d}{\widehat{W}}_u^v}+\sqrt{1-{\widehat{\rho }}^2}{\displaystyle {\int }_t^T\sqrt{{\widehat{v}}_u}\text{d}{\widehat{W}}_u}\times \exp -\frac{\beta }{2}{\displaystyle {\int }_t^T{\sigma }_S^2+{\xi }^2{\epsilon }^{2H-1}{v}_u+{\widehat{v}}_u\text{d}u}.\end{array}$$

We denote $\gamma =1-i\varphi$ and rearrange it to obtain that$$\begin{array}{c}\text{(28)}& f_2\varphi =c_t{\mathbb{E}}_t^{\mathrm{\mathbb{P}}}\exp -\frac{\beta \gamma {\sigma }_S^2\tau }{2}-\frac{\beta \xi }{{\sigma }_v}{v}_t+\theta \tau -\frac{\beta \widehat{\rho }}{{\sigma }_{\widehat{v}}}{\widehat{v}}_t+\widehat{\theta }\tau \times \exp \frac{\beta \xi }{{\sigma }_v}{v}_T+\beta \xi \frac{\kappa }{{\sigma }_v}-\frac{\xi {\epsilon }^{2H-1}}{2}{\displaystyle {\int }_t^T{v}_u\text{d}u}+\frac{\beta \widehat{\rho }}{{\sigma }_{\widehat{v}}}{\widehat{v}}_T+\frac{{\beta }^{2}1-{\widehat{\rho }}^2}{2}+\frac{\beta \widehat{\kappa }\widehat{\rho }}{{\sigma }_{\widehat{v}}}-\frac{\beta }{2}{\displaystyle {\int }_t^T{\widehat{v}}_u\text{d}u}.\end{array}$$

Moreover,$$\begin{array}{c}\text{(29)}& f_2\varphi =c_t\exp -\frac{\beta \gamma {\sigma }_S^2\tau }{2}-\frac{\beta \xi }{{\sigma }_v}{v}_t+\theta \tau -\frac{\beta \widehat{\rho }}{{\sigma }_{\widehat{v}}}{\widehat{v}}_t+\widehat{\theta }\tau \times {\mathbb{E}}_t^{\mathrm{\mathbb{P}}}\frac{\beta \xi }{{\sigma }_v}{v}_T+\beta \xi \frac{\kappa }{{\sigma }_v}-\frac{\xi {\epsilon }^{2H-1}}{2}{\displaystyle {\int }_t^T{v}_u\text{d}u+\frac{\beta \widehat{\rho }}{{\sigma }_{\widehat{v}}}{\widehat{v}}_T}+\frac{{\beta }^{2}1-{\widehat{\rho }}^2}{2}+\frac{\beta \widehat{\kappa }\widehat{\rho }}{{\sigma }_{\widehat{v}}}-\frac{\beta }{2}{\displaystyle {\int }_t^T{\widehat{v}}_u\text{d}u}.\end{array}$$

Therefore, we can obtain the following equation:$$\begin{array}{c}\text{(30)}& f_2\varphi =c_t\exp -\frac{\beta \gamma {\sigma }_S^2\tau }{2}-\frac{\beta \xi }{{\sigma }_v}{v}_t+\theta \tau -\frac{\beta \widehat{\rho }}{{\sigma }_{\widehat{v}}}{\widehat{v}}_t+\widehat{\theta }\tau \times {\mathbb{E}}_t^{\mathrm{\mathbb{P}}}\exp -q_1{v}_T-q_2{\displaystyle {\int }_t^T{v}_u\text{d}u}-q_3{\widehat{v}}_T-q_4{\displaystyle {\int }_t^T{\widehat{v}}_u\text{d}u}.\end{array}$$where $q_1$, $q_2$, $q_3$, and $q_4$ are given by equation (16).

To complete the formulae for $f_1\varphi$ and $f_2\varphi$, we present the calculation of the expectations in the formulae for $f_1\varphi$ and $f_2\varphi$.

Lemma 2.

If the processes $v$ and $\widehat{v}$ follow the dynamics in the model with independent Brownian motions $W^v$ and ${\widehat{W}}^v$, if$$\begin{array}{}\text{(31)}& F\tau ,v_t,{\widehat{v}}_t={\mathbb{E}}_t^{\mathrm{\mathbb{P}}}\exp -{\lambda }_1{v}_T-{\mu }_1{\displaystyle {\int }_t^T{v}_u\text{d}u}-{\lambda }_2{\widehat{v}}_T-{\mu }_2{\displaystyle {\int }_t^T{\widehat{v}}_u\text{d}u},\end{array}$$then the formula for it is given by$$\begin{array}{}\text{(32)}& F\tau ,v_t,{\widehat{v}}_t=\exp -G_1\tau ,{\lambda }_1,{\mu }_1{v}_t-G_2\tau ,{\lambda }_2,{\mu }_2{\widehat{v}}_t-\theta H_1\tau ,{\lambda }_1,{\mu }_1-\widehat{\theta }{H}_2\tau ,{\lambda }_2,{\mu }_2,\end{array}$$where$$\begin{array}{c}\text{(33)}& G_1\tau ,{\lambda }_1,{\mu }_1=\frac{{\lambda }_1{\gamma }_1+\kappa +e^{{\gamma }_1\tau }{\gamma }_1-\kappa +2{\mu }_1{e}^{{\gamma }_1\tau }-1}{{\epsilon }^{2H-1}{\sigma }_v^2{\lambda }_1{e}^{{\gamma }_1\tau }-1+{\gamma }_1-\kappa +e^{{\gamma }_1\tau }{\gamma }_1+\kappa },G_2\tau ,{\lambda }_2,{\mu }_2=\frac{{\lambda }_2{\gamma }_2+\widehat{\kappa }+e^{{\gamma }_2\tau }{\gamma }_2-\widehat{\kappa }+2{\mu }_2{e}^{{\gamma }_2\tau }-1}{{\sigma }_{\widehat{v}}^2{\lambda }_2{e}^{{\gamma }_2\tau }-1+{\gamma }_2-\widehat{\kappa }+e^{{\gamma }_2\tau }{\gamma }_2+\widehat{\kappa }},\\ H_1\tau ,{\lambda }_1,{\mu }_1=-\frac{2}{{\epsilon }^{2H-1}{\sigma }_v^2}\ln \frac{2{\gamma }_1{e}^{\kappa +{\gamma }_1\tau /2}}{{\epsilon }^{2H-1}{\sigma }_v^2{\lambda }_1{e}^{{\gamma }_1\tau }-1+{\gamma }_1-\kappa +e^{{\gamma }_1\tau }{\gamma }_1+\kappa },\\ H_2\tau ,{\lambda }_2,{\mu }_2=-\frac{2}{{\sigma }_{\widehat{v}}^2}\ln \frac{2{\gamma }_2{e}^{\widehat{\kappa }+{\gamma }_2\tau /2}}{{\sigma }_{\widehat{v}}^2{\lambda }_2{e}^{{\gamma }_2\tau }-1+{\gamma }_2-\widehat{\kappa }+e^{{\gamma }_2\tau }{\gamma }_2+\widehat{\kappa }},\\ {\gamma }_1=\sqrt{{\kappa }^2+2{\epsilon }^{2H-1}{\sigma }_v^2{\mu }_1},\\ {\gamma }_2=\sqrt{{\widehat{\kappa }}^2+2{\sigma }_{\widehat{v}}^2{\mu }_2}.\end{array}$$

Proof.

For $t\in 0,T$, we denote$$\begin{array}{}\text{(34)}& M_t=F\tau ,v_t,{\widehat{v}}_t\exp -{\mu }_1{\displaystyle {\int }_0^t{v}_u\text{d}u}-{\mu }_2{\displaystyle {\int }_0^t{\widehat{v}}_u\text{d}u}.\end{array}$$

It can also be expressed as$$\begin{array}{}\text{(35)}& M_t={\mathbb{E}}_t^{\mathrm{\mathbb{P}}}\exp -{\lambda }_1{v}_T-{\mu }_1{\displaystyle {\int }_0^T{v}_u\text{d}u}-{\lambda }_2{\widehat{v}}_t-{\mu }_2{\displaystyle {\int }_0^T{\widehat{v}}_u\text{d}u}.\end{array}$$

$M_t$ is an $\mathbb{F}$-martingale under $\mathrm{\mathbb{P}}$ [26, 30]. $F\tau ,v,\widehat{v}$ satisfies the following PIDE by applying Feynman–Kac theorem:$$\begin{array}{}\text{(36)}& -\frac{\partial F}{\partial \tau }+\theta -\kappa v\frac{\partial F}{\partial v}+\widehat{\theta }-\widehat{\kappa }\widehat{v}\frac{\partial F}{\partial \widehat{v}}+\frac{1}{2}{\epsilon }^{2H-1}{\sigma }_v^{2}v\frac{{\partial }^{2}F}{\partial v^2}+\frac{1}{2}{\sigma }_{\widehat{v}}^2\widehat{v}\frac{{\partial }^{2}F}{\partial {\widehat{v}}^2}-{\mu }_{1}v+{\mu }_2\widehat{v}F=0,\end{array}$$with the boundary condition $F0,v,\widehat{v}=\exp -{\lambda }_{1}v-{\lambda }_2\widehat{v}$; we recall that ${\psi }_t$ is a martingale and ${\psi }_0={\mathbb{E}}^{\mathrm{\mathbb{P}}}{\psi }_t=0$ is used in the above equation.

We conjecture that the solution is$$\begin{array}{}\text{(37)}& F\tau ,v,\widehat{v}=\exp -G_1\tau ,{\lambda }_1,{\mu }_{1}v-G_2\tau ,{\lambda }_2,{\mu }_2\widehat{v}-\theta H_1\tau ,{\lambda }_1,{\mu }_1-\widehat{\theta }{H}_2\tau ,{\lambda }_2,{\mu }_2,\end{array}$$with the boundary condition $G_{1}0,{\lambda }_1,{\mu }_1={\lambda }_1$, $G_{2}0,{\lambda }_2,{\mu }_2={\lambda }_2,\text{\hspace{0.17em}and\hspace{0.17em}}{H}_{1}0,{\lambda }_1,{\mu }_1=H_{2}0,{\lambda }_2,{\mu }_2=0$.

We rearrange equation (36) by instituting the above expression for conjecture into it$$\begin{array}{c}\text{(38)}& 0=\frac{\partial G_1}{\partial \tau }v+\frac{\partial G_2}{\partial \tau }\widehat{v}+\theta \frac{\partial H_1}{\partial \tau }+\widehat{\theta }\frac{\partial H_2}{\partial \tau }-\theta -\kappa v{G}_1-\widehat{\theta }-\widehat{\kappa }\widehat{v}{G}_2+\frac{1}{2}{\epsilon }^{2H-1}{\sigma }_v^{2}v{G}_1^2+\frac{1}{2}{\sigma }_{\widehat{v}}^2\widehat{v}{G}_2^2-{\mu }_{1}v+{\mu }_2\widehat{v}.\end{array}$$

We can obtain the following ODEs:$$\begin{array}{}\text{(39)}& \begin{array}{c}\frac{\partial G_1}{\partial \tau }=-\frac{1}{2}{\epsilon }^{2H-1}{\sigma }_v^2{G}_1^2-\kappa G_1+{\mu }_1,\\ \frac{\partial H_1}{\partial \tau }=G_1,\\ \frac{\partial G_2}{\partial \tau }=-\frac{1}{2}{\sigma }_{\widehat{v}}^2{G}_2^2-\widehat{\kappa }{G}_2+{\mu }_2,\\ \frac{\partial H_2}{\partial \tau }=G_2.\end{array}\end{array}$$

We can obtain the formulae for $G_1\tau ,{\lambda }_1,{\mu }_1$, $G_2\tau ,{\lambda }_2,{\mu }_2$, $H_1\tau ,{\lambda }_1,{\mu }_1$, and $H_2\tau ,{\lambda }_2,{\mu }_2$ by solving the above ODEs. Therefore, the proof is complete.

4. Calibration

In this section, we use the loss function to estimate the parameters in the Heston model, double Heston model, and the model we developed and compare the performance among the three models. The market data offered on May 10, 2012, we use are the implied volatilities of the put options on Dow Jones Industrial Average ETF [28]. Since there are four maturities, we consider 0.001 as the risk-free interest rate. We set $\epsilon =$1e − 5. We need to estimate the following parameters:$$\begin{array}{}\text{(40)}& \Theta =\kappa ,\theta ,{\sigma }_v,v_0,\widehat{\kappa },\widehat{\theta },{\sigma }_{\widehat{v}},{\widehat{v}}_0,\widehat{\rho },\epsilon ,H,{\sigma }_s.\end{array}$$

As we mentioned before, there are two approaches generally used for calibration, to minimize the absolute value of the mean-squared error (MSE) and the relative value of the mean-squared error (RVMSE) between the quoted and model prices. We consider the first approach for calibration. To be specific, first, we choose to use the bisection algorithm [28] to obtain the market prices since we only have the quoted implied volatilities data. Second, we use the market prices and the first approach to estimate the parameters. Suppose that the total number of the quotes is $N$, for a set number of $N_T$ maturities ${\tau }_t$ ($t=1,\dots ,N_T$), there are a set number of $N_K$ strikes $K_k$$k=1,\dots ,N_K$. For each ${\tau }_t,K_k$, the corresponding market price is $C_{t,k}$ and the model price is $C_{t,k}^{\mathrm{\Theta }}$. We consider the following loss function:$$\begin{array}{}\text{(41)}& \text{MSE}=\frac{1}{N}{\displaystyle \sum _{t,k}{C_{t,k}-C_{t,k}^{\Theta }}^2}.\end{array}$$

Table 1 shows the calibration results. It demonstrates that the model provides the smallest MSE and IVMSE. The estimation of the parameters also accords with the study based on Christoffersen et al. [10]. For both double Heston model and the model, one factor of the variance offers a higher mean-reverting level and the other factor of the variance offers a lower mean-reverting level. Since double Heston model offers better performance on MSE and IVMSE compared with the Heston model under the loss function we use, we plot the implied volatilities based on the former and the model using the estimated parameters, and the market implied volatilities for different maturities. Figure 1 shows the result. It demonstrates the model gives better performance than double Heston model on fitting with the market implied volatilities for different maturities. Compared with short-term options, the model has a better fit to the market implied volatilities on long-term options.

Table 1

Calibration of the parameters.

Heston	${\kappa }_H$	${\theta }_H$	${\sigma }_H$	$v_{H0}$	${\rho }_H$	MSE	IVMSE
	8.9814	0.0409	0.2970	0.0244	−0.9621	0.02395	5.266e − 5
Double Heston	${\kappa }_{D1}$	${\theta }_{D1}$	${\sigma }_{D1}$	$v_{D10}$	${\rho }_{D1}$
	10.7526	0.0330	0.3613	0.0252	−0.8916	0.01743	4.928e − 5
	${\kappa }_{D2}$	${\theta }_{D2}$	${\sigma }_{D2}$	$v_{D20}$	${\rho }_{D2}$
	0.9491	0.0257	0.0517	0.0003	0.7009
The model	$\kappa$	$\theta$	${\sigma }_v$	$v_0$	${\sigma }_s$
	11.6490	1.0898	0.5346	0.5344	0.0037	0.00277	1.316e − 5
	$\widehat{\kappa }$	$\widehat{\theta }$	${\sigma }_{\widehat{v}}$	${\widehat{v}}_0$	$\widehat{\rho }$
	0.1673	0.1068	0.9906	0.0316	−0.4422
	$\xi$	$H$
	0.4316	0.7976

[figures omitted; refer to PDF]

We also give some numerical results to examine the impact of the approximation factor $\epsilon$ and the long-memory parameter $H$ on the call option prices. We set the values of parameters as follows: S = 100, K = 100, $r$ = 0.005, and $\xi =0.5$. The values of the other parameters are set in Table 2.

Table 2

The values of the parameters.

Parameters	Values	Parameters	Values	Parameters	Values	Parameters	Values
$\theta$	2.2	$\kappa$	1.5	${\sigma }_v$	0.5	$v_0$	0.5
$\widehat{\theta }$	1.2	$\widehat{\kappa }$	0.5	${\sigma }_{\widehat{v}}$	0.9	${\widehat{v}}_0$	0.25
$\widehat{\rho }$	−0.5	${\sigma }_s$	0.15

Table 3 shows the result. It demonstrates the impact of the parameters $\epsilon$ and $H$ on the prices of the call options. When the value of $\epsilon$ is fixed, the call option price decreases as $H$ increases. When the value of $H$ is fixed, the call option price decreases as $\epsilon$ decreases.

Table 3

The prices of the call options.

$T$		1M	2M	3M	6M	9M	1Y
$\epsilon$	$H$
1e − 5	0.5	7.7304	11.5134	14.7379	23.0400	30.2008	36.5632
	0.6	6.5633	9.8007	12.5793	19.8259	26.1925	31.9473
	0.7	6.4345	9.6118	12.3414	19.4736	25.7555	31.4459
	0.8	6.4215	9.5927	12.3174	19.4380	25.7113	31.3953
1e − 6	0.5	7.7304	11.5134	14.7379	23.0400	30.2008	36.5632
	0.6	6.5108	9.7237	12.4823	19.6823	26.0143	31.7427
	0.7	6.4258	9.5990	12.3254	19.4498	25.7259	31.4120
	0.8	6.4204	9.5911	12.3154	19.4350	25.7076	31.3910
1e − 7	0.5	7.7304	11.5134	14.7379	23.0400	30.2008	36.5632
	0.6	6.4775	9.6748	12.4208	19.5911	25.9011	31.6129
	0.7	6.4223	9.5939	12.3189	19.4403	25.7142	31.3985
	0.8	6.4201	9.5907	12.3149	19.4343	25.7067	31.3900

5. Conclusion

We developed a double Heston model with approximative fractional stochastic volatility in this article. It is an extensive study based on the study of Ahlip et al. [26] and Kang et al. [22]. Since approximative fractional Brownian motion can replace fractional Brownian motion, we introduced it to double Heston model by modeling the dynamics of the stock price and one factor of the variance with approximative fractional process and it is our contribution to this article. We did some pioneering work to study this innovative model. We derived the pricing formula for the options and used the expectation approach to derive the characteristic functions. We did the calibration to estimate the parameters using the market data. The calibration demonstrates that the model we developed provides the smallest MSE and IVMSE. In both double Heston model and the model, one factor of the variance has a higher mean-reverting level and the other factor of the variance has a lower mean-reverting level. The result accords with the study of Christoffersen et al. [10]. Since double Heston model offers better performance on MSE and IVMSE compared with the Heston model under the loss function we use, we plotted the implied volatilities based on the former and the model using the estimated parameters, and the market implied volatilities for different maturities. It demonstrates the model gives better performance than double Heston model on fitting with the market implied volatilities for different maturities. Compared with short-term options, the model has a better fit to the market implied volatilities on long-term options. We also examined the impact of the positive approximation factor and the long-memory parameter on the call option prices. The result demonstrates that when the value of the positive approximation factor is fixed, the call option price decreases as the long-memory parameter increases. When the value of the long-memory parameter is fixed, the call option price decreases as the positive approximation factor decreases.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (grant no. 11601420) and the Natural Science Foundation of Shaanxi Province, China (grant no. 2020JM-577).

References

[1] F. Black, M. Scholes, "The pricing of options and corporate liabilities," Journal of Political Economy, vol. 81 no. 3, pp. 637-654, DOI: 10.1086/260062, 1973.

[2] J. Hull, A. White, "The pricing of options on assets with stochastic volatilities," The Journal of Finance, vol. 42 no. 2, pp. 281-300, DOI: 10.1111/j.1540-6261.1987.tb02568.x, 1987.

[3] L. O. Scott, "Option pricing when the variance changes randomly: theory, estimation, and an application," The Journal of Financial and Quantitative Analysis, vol. 22 no. 4, pp. 419-438, DOI: 10.2307/2330793, 1987.

[4] J. B. Wiggins, "Option values under stochastic volatility: theory and empirical estimates," Journal of Financial Economics, vol. 19 no. 2, pp. 351-372, DOI: 10.1016/0304-405x(87)90009-2, 1987.

[5] E. M. Stein, J. C. Stein, "Stock price distributions with stochastic volatility: an analytic approach," Review of Financial Studies, vol. 4 no. 4, pp. 727-752, DOI: 10.1093/rfs/4.4.727, 1991.

[6] A. L. Lewis, Option Valuation under Stochastic Volatility with Mathematica Code, 2000.

[7] S. L. Heston, "A closed-form solution for options with stochastic volatility with applications to bond and currency options," Review of Financial Studies, vol. 6 no. 2, pp. 327-343, DOI: 10.1093/rfs/6.2.327, 1993.

[8] J. C. Cox, J. E. Ingersoll, S. A. Ross, "A theory of the term structure of interest rates," Econometrica, vol. 53 no. 2, pp. 385-407, DOI: 10.2307/1911242, 1985.

[9] G. Bakshi, C. Cao, Z. Chen, "Empirical performance of alternative option pricing models," The Journal of Finance, vol. 52 no. 5, pp. 2003-2049, DOI: 10.1111/j.1540-6261.1997.tb02749.x, 1997.

[10] P. Christoffersen, S. Heston, K. Jacobs, "The shape and term structure of the index option smirk: why multifactor stochastic volatility models work so well," Management Science, vol. 55 no. 12, pp. 1914-1932, DOI: 10.1287/mnsc.1090.1065, 2009.

[11] Y. Hu, B. Øksendal, "Fractional white noise calculus and applications to finance," Infinite Dimensional Analysis, Quantum Probability and Related Topics, vol. 06 no. 01,DOI: 10.1142/s0219025703001110, 2003.

[12] T. Björk, H. Hult, "A note on Wick products and the fractional Black-Scholes model," Finance and Stochastics, vol. 9, pp. 197-209, 2005.

[13] P. Cheridito, "Mixed fractional Brownian motion," Bernoulli, vol. 7 no. 6, pp. 913-934, DOI: 10.2307/3318626, 2001.

[14] L. Sun, "Pricing currency options in the mixed fractional Brownian motion," Physica A: Statistical Mechanics and Its Applications, vol. 392 no. 16, pp. 3441-3458, DOI: 10.1016/j.physa.2013.03.055, 2013.

[15] Z. Guo, H. Yuan, "Pricing European option under the time-changed mixed Brownian-fractional Brownian model," Physica A: Statistical Mechanics and Its Applications, vol. 406, pp. 73-79, DOI: 10.1016/j.physa.2014.03.032, 2014.

[16] K.-H. Kim, S. Yun, N.-U. Kim, J.-H. Ri, "Pricing formula for European currency option and exchange option in a generalized jump mixed fractional Brownian motion with time-varying coefficients," Physica A: Statistical Mechanics and Its Applications, vol. 522 no. C, pp. 215-231, DOI: 10.1016/j.physa.2019.01.145, 2019.

[17] T. H. Thao, "An approximate approach to fractional analysis for finance," Nonlinear Analysis: Real World Applications, vol. 7 no. 1, pp. 124-132, DOI: 10.1016/j.nonrwa.2004.08.012, 2006.

[18] A. Intarasit, P. Sattayatham, "A geometric Brownian motion model with compound Poisson process and fractional stochastic volatility," Advances and Application in Statistics, vol. 16 no. 1, pp. 25-47, 2010.

[19] P. Sattayatham, A. Intarasit, "An approximate formula of European option for fractional stochastic volatility jump-diffusion model," Journal of Mathematics and Statistics, vol. 7 no. 3, pp. 230-238, DOI: 10.3844/jmssp.2011.230.238, 2011.

[20] M. Mrázek, J. Pospíšil, T. Sobotka, "On calibration of stochastic and fractional stochastic volatility models," European Journal of Operational Research, vol. 254 no. 3, pp. 1036-1046, 2016.

[21] J. Pospíšil, T. Sobotka, "Market calibration under a long memory stochastic volatility model," Applied Mathematical Finance, vol. 23 no. 5, pp. 323-343, 2016.

[22] J. Kang, B. Yang, N. Huang, "Pricing of FX options in the MPT/CIR jump-diffusion model with approximative fractional stochastic volatility," Physica A: Statistical Mechanics and Its Applications, 2019.

[23] D. Duffie, J. Pan, K. Singleton, "Transform analysis and asset pricing for affine jump-diffusions," Econometrica, vol. 68 no. 6, pp. 1343-1376, DOI: 10.1111/1468-0262.00164, 2000.

[24] R. Schöbel, J. Zhu, "Stochastic volatility with an Ornstein-Uhlenbeck process: an extension," Review of Finance, vol. 3 no. 1, pp. 23-46, DOI: 10.1023/a:1009803506170, 1999.

[25] J. Zhu, Applications of Fourier Transform to Smile Modeling: Theory and Implementation, 2010.

[26] R. Ahlip, L. A. F. Park, A. Prodan, "Semi-Analytical option pricing under double Heston jump-diffusion hybrid model," Journal of Mathematical Sciences and Modelling, vol. 1 no. 3, pp. 138-152, DOI: 10.33187/jmsm.432019, 2018.

[27] X.-J. He, S.-P. Zhu, "A closed-form pricing formula for European options under the Heston model with stochastic interest rate," Journal of Computational and Applied Mathematics, vol. 335, pp. 323-333, DOI: 10.1016/j.cam.2017.12.011, 2018.

[28] F. D. Rouah, The Heston Model and its Extensions in Matlab and C#, 2013.

[29] P. Christoffersen, K. Jacobs, "The importance of the loss function in option valuation," Journal of Financial Economics, vol. 72 no. 2, pp. 291-318, DOI: 10.1016/j.jfineco.2003.02.001, 2004.

[30] R. J. Elliott, P. E. Kopp, Mathematics of Financial Markets, 2005.