Zhen Yuan¹, Haomin Zhang^1,2*, Songyu Hong^3
1School of Mathematics and Statistics, Guilin University of Technology, Guilin, China.
²Guangxi Colleges and Universities Key Laboratory of Applied Statistics, Guilin University of Technology, Guilin, China.
³School of Economics and Statistics, Guangzhou University, Guangzhou, China.
DOI: 10.4236/ajibm.2026.164023   PDF    HTML   XML   60 Downloads   289 Views   Citations

Abstract

This paper proposes a novel stochastic volatility model for pricing European call options. Based on the double Heston model, the model introduces a stochastic long-term average process and additional volatility terms for each volatility component, and assumes that the long-term mean itself has dynamic evolution characteristics. Moreover, the model regulates some key parameters through a Markov state transition mechanism. This study uses the characteristic function method to derive a closed-form pricing formula for European call options. The numerical accuracy of the formula is verified through Monte Carlo simulation, and further numerical experiments are conducted to average the process. Finally, based on a rigorously designed empirical analysis, it is shown that the proposed model outperforms the two comparison models in terms of option pricing accuracy.

Keywords

Regime-Switching, Double Heston Model, Nonlinear, Closed-Form Solution, European Call Option

Share and Cite:

1. Introduction

European options represent one of the most significant derivative instruments in the global financial markets. To accurately price these instruments, numerous scholars have developed various pricing models. However, as most of these models rely on numerical simulation methods—which demand substantial computational resources—the development of a closed-form pricing solution for such securities remains of paramount importance for practitioners in the financial industry. Since the Black and Scholes (B-S) model (Black & Scholes, 1973) was proposed, it has still been widely adopted by financial practitioners. In the B-S model, the underlying asset price follows a lognormal distribution, and there is a simple formula to calculate the price of European options. However, despite the popularity of the B-S model, some of the oversimplified assumptions it makes to achieve analytical tractability lead to potential mispricing problems. The existence of the “volatility smile” (Dumas et al., 1998) demonstrates the unrealistic nature of the model’s assumption of constant volatility. Non-constant volatility models have been significantly developed to improve the model further. Non-constant volatility models can be divided into two main categories: local volatility models and stochastic volatility models. The former was proposed by Dupire (Dupire, 1994); in this model, volatility is defined as a deterministic function of the underlying price and time. However, numerous empirical studies have demonstrated that the local volatility model cannot capture the “volatility smile” (Hagan et al., 2002). This limitation has led to the growing popularity of stochastic volatility models.

The stochastic volatility model complicates the discovery of a closed-form pricing formula for European options because it introduces volatility as a new random variable. As a result, most existing models must rely on numerical methods. Scott (Scott, 1987) uses the Monte Carlo simulation to calculate the option price. Wiggins (Wiggins, 1987) uses the finite difference method to calculate the option price. However, these numerical methods often need a long time to obtain the results of option pricing, making them inapplicable in actual financial markets due to the extensive time required for model calibration and option valuation. Therefore, stochastic volatility models with closed-form solutions have become a better and more relevant research direction for the needs of the financial market. Hull and White (Hull & White, 1987) derived a series solution under their model in which one volatility follows another geometric Brownian motion. However, the model is still unsatisfactory. On one hand, the assumption of independence between the underlying price and volatilities violates the leverage effect, as empirical studies confirm a negative correlation between the underlying price and volatility (Bakshi et al., 1997; Jacquier et al., 2004). On the other hand, the volatility process of the model lacks a mean-reverting property, contradicting the fact that the volatility process is mean-reverting (Beckers, 1983).

Heston (Heston, 1993) proposed a more effective and widely utilized model among financial practitioners, which enables the derivation of a closed-form pricing formula for European call options through the Cox-Ingersoll-Ross (CIR) process. The model not only satisfies the assumption of arbitrary correlation between the underlying price and volatility but also upholds the non-negative and mean-reverting properties of the volatility process. More importantly, due to the existence of an analytical solution, the Heston model has overcome the shortcomings associated with numerical methods, which require considerable time and energy for model calibration and option valuation when applied to real financial markets. Despite its numerous advantages, the Heston model is not without flaws, as it exhibits significant nonlinear mean-reversion phenomena in asset volatility (Bakshi et al., 2006). To mitigate this issue as much as possible, He and Chen (He & Chen, 2021) replace the constant term of the long-term mean volatility with a variable that follows a normal distribution. Notably, this model still retains the essential advantage of the Heston model: analytical tractability, allowing for the derivation of closed pricing formulas for European options. To model the volatility structure more flexibly and better fit empirical data to European option prices, Christoffersen et al. (Christoffersen et al., 2009) constructed the double Heston model by adding another stochastic volatility factor based on the CIR process to the Heston model. Mehrdoust et al. (Mehrdoust et al., 2023) and Zhang and Feng (Zhang & Feng, 2019) support this with their research on American options. More profoundly, two mutually independent CIR mean-reversion processes are simultaneously included in the double Heston model, with computational tools very similar to those of the standard Heston model; however, the results obtained from option pricing under this model can be highly desirable.

Regime-switching models have been widely used in the finance field to simulate the prices of financial derivatives that are impacted by economic cycles (Hamilton, 1990; Eraker, 2004). Recently, much literature has applied the regime-switching mechanism to a stochastic volatility model with regime-switching. Elliott and Lian (Elliott & Lian, 2013) provide closed-form exact solutions for pricing discrete sampling variance swaps and volatility swaps based on the Heston stochastic volatility model featuring regime-switching. Lin and He (2021) and He and Lin (2023) constructed two nonlinear stochastic volatility models based on the He-Chen model (He & Chen, 2021) by modelling the long-term mean of volatility and introduced a regime-switching model in each, deriving closed solutions easily. Mehrdoust (Mehrdoust et al., 2022) examined the pricing of American options under a new model developed by incorporating regime-switching at the interest rate level and mean reversion into the double Heston model.

Building on previous research, this paper establishes a new model based on the double Heston model by replacing the mean reversion level of each volatility process with a stochastic long-term mean process and substituting the constant parameter of the stochastic long-term mean process with a regime-switching term governed by a two-state Markov process. If the regime-switching term degenerates into a constant parameter, we obtain a new nonlinear double without regime-switching. We refer to the new models with and without regime-switching as the Markov Regime-Switching Nonlinear Double Heston Model (MRSNDH) and the Nonlinear Double Heston Model (NDH). Both MRSNDH and NDH models retain the fundamental advantages of the He-Chen model and the double Heston model, namely, analytical tractability. In addition to the derivation details, we provide verification of the accuracy of the newly derived formulas. On this basis, the MRSNDH model is compared with the NDH model to explore the effect of introducing regime-switching on European call option prices, while the NDH model is compared with the double Heston model to investigate the impact of incorporating two stochastic long-term mean processes in the double Heston model on European call option prices.

This paper involves the estimation of numerous parameters, and selecting an appropriate estimation method to achieve rapid and accurate parameter estimation under limited computing resources poses a significant challenge. The Particle Swarm Optimisation (PSO) algorithm, proposed by electrical engineer Russell Eberhart and American social psychologist James Kennedy in 1995, is an evolutionary computing method based on swarm intelligence (Kennedy & Eberhart, 1995). In 1998, Shi and Eberhart (Shi & Eberhart, 1998) introduced the inertia weight into the original PSO, using its value to represent the contribution of historical velocity to current velocity, which was later called the standard PSO. He (He, 2017) proposed a particle swarm optimisation algorithm based on correcting the global optimal position to estimate the parameters of the B-S model. There is still room for improvement in the traditional PSO, so Gong and Zhang (Gong & Zhang, 2016) introduced two hybrid optimisation techniques, the hybrid PSO algorithm and the hybrid Differential Evolution (DE) algorithm, into the parameter calibration scheme to enhance the calibration quality of the new model. Ratnaweera et al. (Ratnaweera et al., 2004) proposed an adaptive adjustment strategy based on the most up-to-date information from each individual. This paper will adopt the Adaptive Particle Swarm Optimisation (APSO) to estimate the parameters of MRSNDH.

The remainder of this paper is organized as follows. In Section 2, we introduce the newly proposed model, followed by the closed pricing formula for the European call option based on this model. Section 3 investigates various properties of the new formula through numerical experiments. The final section presents the conclusion. In Section 4, we utilize actual data on European call options related to the S&P 500 index to analyze the performance of different models in practical market applications and demonstrate the superiority of our MRSNDH model over others.

2. The Closed Solution of the MRSNDH Model

In this section, we propose a new model, namely the MRSNDH model, for modelling the price of the underlying asset and option pricing based on the He-Chen model (He & Chen, 2021) and the double Heston model (Christoffersen et al., 2009). This model draws on the design idea of the He-Chen model in introducing a random long-term mean in the volatility process and makes improvements on the basis of the double Heston model. Meanwhile, it further introduces a volatility term in the volatility dynamics to enhance the model’s ability to describe market complexity. First, we introduce the He-Chen model that

$\{\begin{array}{l}\frac{\text{d}{S}_t}{S_t}=r\text{d}t+\sqrt{v_t}\text{d}{W}_t^1,\\ \text{d}{v}_t=\kappa \left({\theta }_t-v_t\right)\text{d}t+{\sigma }_1\sqrt{v_t}\text{d}{W}_t^2,\\ \text{d}{\theta }_t=\lambda \text{d}t+\eta \text{d}{B}_t.\end{array}$ (1)

Here, $t\ge 0$ , $S_t$ and $v_t$ denote the underlying asset price and volatility, respectively. $W_t^1$ and $W_t^2$ are two standard Brownian motions with a correlation coefficient of $\rho$ , $r$ denotes the risk-free interest rate, and $\sigma$ is the so-called volatility of volatility. $\kappa$ represents the speed of mean-reversion. The long-term average of volatility is composed of a stochastic component ${\theta }_t$ . To capture the volatility structure with greater flexibility and price European options more accurately, we propose a nonlinear double-Heston model augmented with a Markov regime-switching mechanism—hereinafter referred to as the MRSNDH model.

$\{\begin{array}{l}\frac{\text{d}{S}_t}{S_t}=r\text{d}t+\sqrt{v_1}\text{d}{W}_t^1+\sqrt{v_2}\text{d}{W}_t^2,\\ \text{d}{v}_1={\kappa }_1\left({\overline{v}}_1+{\theta }_1-v_1\right)\text{d}t+{\sigma }_1\sqrt{v_1}\text{d}{W}_t^3+{\epsilon }_{X_t}\text{d}{Z}_t^1,\\ \text{d}{v}_2={\kappa }_2\left({\overline{v}}_2+{\theta }_2-v_2\right)\text{d}t+{\sigma }_2\sqrt{v_2}\text{d}{W}_t^4+{\delta }_{X_t}\text{d}{Z}_t^2,\\ \text{d}{\theta }_1={\lambda }_{X_t}\text{d}t+{\eta }_{X_t}\text{d}{Z}_t^3,\\ \text{d}{\theta }_2={\alpha }_{X_t}\text{d}t+{\beta }_{X_t}\text{d}{Z}_t^4,\end{array}$ (2)

where, $W_t^1$ and $W_t^3$ have correlation ${\rho }_1$ , $W_t^2$ and $W_t^4$ have correlation ${\rho }_2$ . $Z_t^1$ through $Z_t^4$ are independent standard Brownian motions. The Markov chain $X_t$ is defined as:

$X_t=\{\begin{array}{ll}{\left(1,0\right)}^{\top }\hfill & \text{state}=1,\hfill \\ {\left(0,1\right)}^{\top }\hfill & \text{state}=2.\hfill \end{array}$

Here, states 1 and 2 are employed to characterise the prevailing economic conditions. In this study, the classification of these states is determined by the level of implicit volatility observed in the actual data. The transition between states follows a Poisson process:

$P\left(t_{ij}\ge t\right)={\text{e}}^{-{\lambda }_{ij}t},i,j=1,2,i\ne j,$

where, ${\lambda }_{ij}$ represents the transition rate of the random variable $X_t$ from state $i$ to state $j$ , while $t_{ij}$ denotes the dwell time in state $j$ prior to transitioning to state $i$ . ${\epsilon }_{X_t}$ , ${\delta }_{X_t}$ , ${\lambda }_{X_t}$ , ${\eta }_{X_t}$ , ${\alpha }_{X_t}$ and ${\beta }_{X_t}$ are parameters that depend on the Markov chain, expressed as follows:

${\epsilon }_{X_t}=\langle \widehat{\epsilon },X_t\rangle ,$

${\delta }_{X_t}=\langle \widehat{\delta },X_t\rangle ,$

${\lambda }_{X_t}=\langle \widehat{\lambda },X_t\rangle ,$

${\eta }_{X_t}=\langle \widehat{\eta },X_t\rangle ,$

${\alpha }_{X_t}=\langle \widehat{\alpha },X_t\rangle ,$

${\beta }_{X_t}=\langle \widehat{\beta },X_t\rangle .$

Here, $\widehat{\epsilon }={\left({\epsilon }_1,{\epsilon }_2\right)}^{\text{T}}$, $\widehat{\delta }={\left({\delta }_1,{\delta }_2\right)}^{\text{T}}$, $\widehat{\lambda }={\left({\lambda }_1,{\lambda }_2\right)}^{\text{T}}$, $\widehat{\eta }={\left({\eta }_1,{\eta }_2\right)}^{\text{T}}$, $\widehat{\alpha }={\left({\alpha }_1,{\alpha }_2\right)}^{\text{T}}$, $\widehat{\beta }={\left({\beta }_1,{\beta }_2\right)}^{\text{T}}$, Note: ${\widehat{u}}^{\text{T}}$ denotes the transpose of vector $\widehat{u}$ , and the symbol $\langle \cdot ,\cdot \rangle$ indicates the inner product of two vectors. When the conversion rate is 0, i.e., ${\lambda }_{12}={\lambda }_{21}=0$ and $\widehat{\epsilon }$ , $\widehat{\delta }$ , $\widehat{\lambda }$ , $\widehat{\eta }$ , $\widehat{\alpha }$ and $\widehat{\beta }$ are all set to zero vectors, the MRSNDH model reduces to a special case corresponding to a nonlinear two-dimensional Heston model with no state transitions, which we refer to as the NDH model.

$\{\begin{array}{l}\frac{\text{d}S}{S}=r\text{d}t+\sqrt{v_1}\text{d}{W}_t^1+\sqrt{v_2}\text{d}{W}_t^2,\\ \text{d}{v}_1={\kappa }_1\left({\overline{v}}_1+{\theta }_1-v_1\right)\text{d}t+{\sigma }_1\sqrt{v_1}\text{d}{W}_t^3+\epsilon \text{d}{Z}_t^1,\\ \text{d}{v}_2={\kappa }_2\left({\overline{v}}_2+{\theta }_2-v_2\right)\text{d}t+{\sigma }_2\sqrt{v_2}\text{d}{W}_t^4+\delta \text{d}{Z}_t^2,\\ \text{d}{\theta }_1=\lambda \text{d}t+\eta \text{d}{Z}_t^3,\\ \text{d}{\theta }_2=\alpha \text{d}t+\beta \text{d}{Z}_t^4.\end{array}$ (3)

Theorem 1: Let $U\left(S,v_1,v_2,{\theta }_1,{\theta }_2,X_t,t\right)$ be the European call option price satisfying the MRSNDH model. Then

$U\left(S,v,\theta ,X_t,t\right)=S_t{P}_1-K{\text{e}}^{-r\left(T-t\right)}{P}_2,$ (4)

where

$P_1=\frac{1}{2}+\frac{1}{\pi }{\displaystyle {\int }_0^{\infty }\text{Re}}\left[\frac{\text{exp}\left(-\text{i}\varphi \ln K\right)\overline{f}\left(\varphi ;\tau ,y,v_1,v_2,{\theta }_1,{\theta }_2,X_t\right)}{\text{i}\varphi }\right]\text{d}\varphi ,$

$P_2=\frac{1}{2}+\frac{1}{\pi }{\displaystyle {\int }_0^{\infty }\text{Re}}\left[\frac{\text{exp}\left(-\text{i}\varphi \ln K\right)f\left(\varphi ;\tau ,y,v_1,v_2,{\theta }_1,{\theta }_2,X_t\right)}{\text{i}\varphi }\right]\text{d}\varphi ,$

$\begin{array}{c}f\left(\varphi ;\tau ,y,v_1,v_2,{\theta }_1,{\theta }_2,X_t\right)=\text{exp}(\overline{C}\left(\tau ;\varphi \right)+D_1\left(\tau ;\varphi \right)v_1+D_2\left(\tau ;\varphi \right)v_2\\ +E_1\left(\tau ;\varphi \right){\theta }_1+E_2\left(\tau ;\varphi \right){\theta }_2+\text{i}\varphi y)\langle {\text{e}}^M{X}_t,I\rangle ,\end{array}$

$D_1\left(\tau ;\varphi \right)=\frac{{\kappa }_1-{\sigma }_1{\rho }_1\text{i}\varphi -d_1}{{\sigma }_1^2}\cdot \frac{1-{\text{e}}^{-d_1\tau }}{1-c_1{\text{e}}^{-d_1\tau }},$

$D_2\left(\tau ;\varphi \right)=\frac{{\kappa }_2-{\sigma }_2{\rho }_2\text{i}\varphi -d_2}{{\sigma }_2^2}\cdot \frac{1-{\text{e}}^{-d_2\tau }}{1-c_2{\text{e}}^{-d_2\tau }},$

$E_1\left(\tau ;\varphi \right)=\frac{{\kappa }_1}{{\sigma }_1^2}\left[\left({\kappa }_1-{\sigma }_1{\rho }_1\text{i}\varphi -d_1\right)\tau -2\ln \left(\frac{1-c_1{\text{e}}^{-d_1\tau }}{1-c_1}\right)\right],$

$E_2\left(\tau ;\varphi \right)=\frac{{\kappa }_2}{{\sigma }_2^2}\left[\left({\kappa }_2-{\sigma }_2{\rho }_2\text{i}\varphi -d_2\right)\tau -2\ln \left(\frac{1-c_2{\text{e}}^{-d_2\tau }}{1-c_2}\right)\right],$

$\overline{C}\left(\tau ;\varphi \right)=r\text{i}\varphi \tau +{\overline{v}}_1{E}_1\left(\tau ;\varphi \right)+{\overline{v}}_2{E}_2\left(\tau ;\varphi \right),$

$d_1=\sqrt{{\left({\sigma }_1{\rho }_1\text{i}\varphi -{\kappa }_1\right)}^2+{\sigma }_1^2\left(\text{i}\varphi +{\varphi }^2\right)},$

$d_2=\sqrt{{\left({\sigma }_2{\rho }_2\text{i}\varphi -{\kappa }_2\right)}^2+{\sigma }_2^2\left(\text{i}\varphi +{\varphi }^2\right)},$

$c_1=\frac{1}{g_1},c_2=\frac{1}{g_2},g_1=\frac{{\kappa }_1-{\sigma }_1{\rho }_1\text{i}\varphi +d_1}{{\kappa }_1-{\sigma }_1{\rho }_1\text{i}\varphi -d_1},g_2=\frac{{\kappa }_2-{\sigma }_2{\rho }_2\text{i}\varphi +d_2}{{\kappa }_2-{\sigma }_2{\rho }_2\text{i}\varphi -d_2},$

$M=A^{\top }\tau +B,$

$A=\left(\begin{array}{cc}-{\lambda }_{12}& {\lambda }_{12}\\ {\lambda }_{21}& -{\lambda }_{21}\end{array}\right),B=\left(\begin{array}{cc}{p}_1& 0\\ 0& p_2\end{array}\right),$

$\begin{array}{c}{p}_1={\displaystyle {\int }_0^{\tau }({\lambda }_1{E}_1\left(s;\varphi \right)+{\alpha }_1{E}_2\left(s;\varphi \right)+\frac{1}{2}{\epsilon }_1^2{D}_1^2\left(s;\varphi \right)+\frac{1}{2}{\delta }_1^2{D}_2^2\left(s;\varphi \right)}\\ +\frac{1}{2}{\eta }_1^2{E}_1^2\left(s;\varphi \right)+\frac{1}{2}{\beta }_1^2{E}_2^2\left(s;\varphi \right))\text{d}s,\end{array}$

$\begin{array}{c}{p}_2={\displaystyle {\int }_0^{\tau }({\lambda }_2{E}_1\left(s;\varphi \right)+{\alpha }_2{E}_2\left(s;\varphi \right)+\frac{1}{2}{\epsilon }_2^2{D}_1^2\left(s;\varphi \right)+\frac{1}{2}{\delta }_2^2{D}_2^2\left(s;\varphi \right)}\\ +\frac{1}{2}{\eta }_2^2{E}_1^2\left(s;\varphi \right)+\frac{1}{2}{\beta }_2^2{E}_2^2\left(s;\varphi \right))\text{d}s,\end{array}$

$I={\left(1,1\right)}^{\top },\tau =T-t,\text{i}=\sqrt{-1}.$

Proof. Let $K$ represent the strike price and $y$ denote the logarithm of the asset price ($\text{ln}\left(S\right)$ ). The process of deriving European call options $U\left(S,v_1,v_2,{\theta }_1,{\theta }_2,X_t,t\right)$ is as follows:

$\begin{array}{c}U\left(S,v_1,v_2,{\theta }_1,{\theta }_2,X_t,t\right)={\text{e}}^{-r\left(T-t\right)}{\mathbb{E}}^Q\left[\text{max}\left(S_T-K,0\right)|S_t,v_{1t},v_{2t},{\theta }_{1t},{\theta }_{2t},X_t\right]\\ ={\text{e}}^{-r\left(T-t\right)}{\mathbb{E}}^Q\left[\left(S_T-K\right)1_{S_T>K}|y_t,v_{1t},v_{2t},{\theta }_{1t},{\theta }_{2t},X_t\right]\\ ={\text{e}}^{-r\left(T-t\right)}{\mathbb{E}}^Q\left[S_T{1}_{S_T>K}\right]-K{\text{e}}^{-r\left(T-t\right)}{\mathbb{E}}^Q\left[1_{S_T>K}\right].\end{array}$ (5)

Here, $\mathbb{E}$ represents the formula for calculating the expected value. $1_{s_t>k}$ . $p\left(y\right)$ represents the density function of $y_t$ . If this probability density function can be directly calculated, determining the option price becomes straightforward; however, in practice, deriving the density function of the underlying asset’s price is challenging. Consequently, based on the Gil-Peláez theorem, further derivations can be carried out.

$P_2\triangleq {\displaystyle {\int }_{\text{ln}\left(K\right)}^{+\infty }p\left(y\right)\text{d}y}=\frac{1}{2}+\frac{1}{\pi }{\displaystyle {\int }_0^{\infty }\text{Re}}\left[\frac{{\text{e}}^{-i\varphi \ln \left(K\right)}f\left(\varphi ;\tau ,y,v_1,v_2,{\theta }_1,{\theta }_2,X_t\right)}{i\varphi }\right]\text{d}\varphi ,$

where, $\text{Re}\left[a+bi\right]$ denotes real part $a$ of complex number $a+bi$ . $f\left(\varphi ;\tau ,y,v_1,v_2,{\theta }_1,{\theta }_2,X_t\right)$ denotes the characteristic function of $y_T$ . It is important to note that we define a new measure, $Q^S$ :

$\frac{\text{d}{Q}^S}{\text{d}Q}=\frac{S_T/S_t}{B_T/B_t}=\frac{{\text{e}}^{y_T}}{{\mathbb{E}}^Q\left[{\text{e}}^{y_T}\right]},$ (6)

where, $B_t={\text{e}}^{{\displaystyle {\int }_0^{t}r\text{d}u}}={\text{e}}^{rt}$ . Therefore, we can derive

${\mathbb{E}}^Q\left[S_T{1}_{S_T>K}\right]=S_t{\text{e}}^{r\left(T-t\right)}=\mathbb{E}\left[S_T\right].$

Under the risk-neutral ($Q$ ) measure, asset prices evolve deterministically at the

risk-free rate; consequently, ${\text{e}}^{-r\left(T-t\right)}{\mathbb{E}}^Q\left[S_T{1}_{S_T>K}\right]$ simplifies to

$\begin{array}{c}{\text{e}}^{-r\left(T-t\right)}{\mathbb{E}}^Q\left[S_T{1}_{S_T>K}\right]=S_t{\mathbb{E}}^Q\left[\frac{S_T/S_t}{B_T/B_t}{1}_{S_T>K}\right]=S_t{\mathbb{E}}^{Q^S}\left[\frac{S_T/S_t}{B_T/B_t}{1}_{S_T>K}\frac{\text{d}Q}{\text{d}{Q}^S}\right]\\ =S_t{\mathbb{E}}^{Q^S}\left[1_{S_T>K}\right].\end{array}$

According to Formula (6), the new density function $p^S\left(y\right)$ is defined as $p\left(y\right)$ via the Radon-Nikodym derivative with respect to

$p^S\left(y\right)\text{d}y=\frac{{\text{e}}^y}{{\mathbb{E}}^Q\left[{\text{e}}^{y_T}\right]}p\left(y\right)\text{d}y.$

Therefore, the characteristic function of $p^S\left(y\right)$ is given by:

${\mathbb{E}}^{{\text{Q}}^{\text{S}}}\left[{\text{e}}^{i\varphi y}\right]={\displaystyle {\int }_{-\infty }^{\infty }{\text{e}}^{i\varphi y}p\left(y\right)\text{d}y}=\frac{1}{{\mathbb{E}}^Q\left[{\text{e}}^{y_T}\right]}{\displaystyle {\int }_{-\infty }^{\infty }{\text{e}}^{i\varphi y}{\text{e}}^{y}p\left(y\right)\text{d}y},$

where ${\mathbb{E}}^Q\left[{\text{e}}^{y_T}\right]$ is a constant; moreover, given that the characteristic function of $y_T$ is $f\left(\varphi ;\tau ,y,v_1,v_2,{\theta }_1,{\theta }_2,X_t\right)={\mathbb{E}}^Q\left[{\text{e}}^{i\varphi y_T}\right]$ , it follows that

${\mathbb{E}}^Q\left[{\text{e}}^{y_T}\right]=f\left(-i;\tau ,y,v_1,v_2,{\theta }_1,{\theta }_2,X_t\right)=S{\text{e}}^{r\left(T-t\right)},$

and

${\displaystyle {\int }_{-\infty }^{\infty }{\text{e}}^{i\varphi y}{\text{e}}^{y}p\left(y\right)\text{d}y}={\mathbb{E}}^Q\left[{\text{e}}^{i\left(\varphi -i\right)y_T}\right],$

Therefore, the characteristic function of the density function $p^S\left(y\right)$ is given by

${\mathbb{E}}^{Q^S}\left[{\text{e}}^{i\varphi y_T}\right]=\frac{f\left(\varphi -i;\tau ,y,v_1,v_2,{\theta }_1,{\theta }_2,X_t\right)}{f\left(-i;\tau ,y,v_1,v_2,{\theta }_1,{\theta }_2,X_t\right)}\triangleq \overline{f}\left(\varphi ;\tau ,y,v_1,v_2,{\theta }_1,{\theta }_2,X_t\right).$

For $\overline{f}\left(\varphi ;\tau ,y,v_1,v_2,{\theta }_1,{\theta }_2,X_t\right)$ , the Gil-Peláez theorem yields:

$P_1\triangleq {\mathbb{E}}^{Q^S}\left[1_{S_T>K}\right]=\frac{1}{2}+\frac{1}{\pi }{\displaystyle {\int }_0^{\infty }\text{Re}}\left[\frac{{\text{e}}^{-i\varphi \ln \left(K\right)}\overline{f}\left(\varphi -i;\tau ,y,v_1,v_2,{\theta }_1,{\theta }_2,X_t\right)}{i\varphi }\right]\text{d}\varphi$

Our next objective is to define the characteristic function $f\left(\varphi ;\tau ,y,v_1,v_2,{\theta }_1,{\theta }_2,X_t\right)$ , which will allow us to determine the final price of the European call option. To determine the analytical solution of the characteristic function, we begin with its definition.

$f\left(\varphi ;\tau ,y,v_1,v_2,{\theta }_1,{\theta }_2,X_t\right)=\mathbb{E}\left[{\text{e}}^{i\varphi y_T}|y_t,v_{1t},v_{2t},{\theta }_{1t},{\theta }_{2t},X_t\right],$

based on the expected value formula, this characteristic function can also be represented as

$f\left(\varphi ;\tau ,y,v_1,v_2,{\theta }_1,{\theta }_2,X_t\right)=\mathbb{E}\left[\mathbb{E}\left({\text{e}}^{i\varphi y_T}|y_t,v_{1t},v_{2t},{\theta }_{1t},{\theta }_{2t},X_T\right)|y_t,v_{1t},v_{2t},{\theta }_{1t},{\theta }_{2t},X_t\right],$

this expression suggests that we should initially compute the inner expectation, assuming the Markov chain is known prior to the maturity date. The inner expectation can be defined as

$h\left(\varphi ;y_t,v_{1t},v_{2t},{\theta }_{1t},{\theta }_{2t},X_t|X_T\right)=\mathbb{E}\left[{\text{e}}^{i\varphi y_T}|y_t,v_{1t},v_{2t},{\theta }_{1t},{\theta }_{2t},X_T\right].$

Based on the data from the Markov chain, parameters ${\kappa }_1$ , ${\kappa }_2$ , ${\sigma }_1$ , ${\sigma }_2$ , ${\rho }_1$ and ${\rho }_2$ are deterministic functions dependent solely on time. Consequently, these parameters can be denoted as ${\kappa }_1\left(t\right)$ , ${\kappa }_2\left(t\right)$ , ${\sigma }_1\left(t\right)$ , ${\sigma }_2\left(t\right)$ , ${\rho }_1\left(t\right)$ and ${\rho }_2\left(t\right)$ , respectively. According to the Feynman-Kac theorem, $f$ must satisfy the following Partial Differential Equation (PDE):

$\begin{array}{c}\frac{\partial h}{\partial \tau }=\left[r-\frac{1}{2}\left(v_1+v_2\right)\right]\frac{\partial h}{\partial x}+{\kappa }_1\left({\overline{v}}_1+{\theta }_1-v_1\right)\frac{\partial h}{\partial v_1}+{\kappa }_2\left({\overline{v}}_2+{\theta }_2-v_2\right)\frac{\partial h}{\partial v_2}\\ +{\lambda }_{X_t}\frac{\partial h}{\partial {\theta }_1}+{\alpha }_{X_t}\frac{\partial h}{\partial {\theta }_2}+\frac{1}{2}\left(v_1+v_2\right)\frac{{\partial }^{2}h}{\partial x^2}+\frac{1}{2}\left({\sigma }_1^2{v}_1+{\epsilon }_{X_t}^2\right)\frac{{\partial }^{2}h}{\partial v_1^2}\\ +\frac{1}{2}\left({\sigma }_2^2{v}_2+{\delta }_{X_t}^2\right)\frac{{\partial }^{2}h}{\partial v_2^2}+\frac{1}{2}{\eta }_{X_t}^2\frac{{\partial }^{2}h}{\partial {\theta }_1^2}+\frac{1}{2}{\beta }_{X_t}^2\frac{{\partial }^{2}h}{\partial {\theta }_2^2}\\ +{\sigma }_1{v}_1{\rho }_1\frac{{\partial }^{2}h}{\partial x\partial v_1}+{\sigma }_2{v}_2{\rho }_2\frac{{\partial }^{2}h}{\partial x\partial v_2}.\end{array}$ (7)

Given the previously established conditions, we obtain

$h\left(\varphi ;0,y,v_1,v_2,{\theta }_1,{\theta }_2,X_t|X_T\right)={\text{e}}^{i\varphi y_T}.$ (8)

Based on the results of previous researchers (He & Chen, 2021; Christoffersen et al., 2009; Ratnaweera et al., 2004), the internal function can be structured as follows.

$h\left(\varphi ;\tau ,y,v_1,v_2,{\theta }_1,{\theta }_2,X_t|X_T\right)={\text{e}}^{C\left(\tau ;\varphi \right)+D_1\left(\tau ;\varphi \right)v_1+D_2\left(\tau ;\varphi \right)v_2+E_1\left(\tau ;\varphi \right){\theta }_1+E_2\left(\tau ;\varphi \right){\theta }_2+i\varphi y},$ (9)

when substituted into the PDF equation above, we obtain the following five Ordinary Differential Equations (ODEs).

$\begin{array}{l}\frac{\partial D_1}{\partial \tau }=\frac{1}{2}{\sigma }_1^2{D}_1{\left(\tau ;\varphi \right)}^2-\left({\kappa }_1-i\varphi {\rho }_1{\sigma }_1\right)D_1\left(\tau ;\varphi \right)-\frac{1}{2}\left(i\varphi +{\varphi }^2\right),\\ \frac{\partial D_2}{\partial \tau }=\frac{1}{2}{\sigma }_2^2{D}_2{\left(\tau ;\varphi \right)}^2-\left({\kappa }_2-i\varphi {\rho }_2{\sigma }_2\right)D_2\left(\tau ;\varphi \right)-\frac{1}{2}\left(i\varphi +{\varphi }^2\right),\\ \frac{\partial E_1}{\partial \tau }={\kappa }_1{D}_1\left(\tau ;\varphi \right),\\ \frac{\partial E_2}{\partial \tau }={\kappa }_2{D}_2\left(\tau ;\varphi \right),\\ \frac{\partial C}{\partial \tau }=ri\varphi +{\kappa }_1{\overline{v}}_1{D}_1\left(\tau ;\varphi \right)+{\kappa }_2{\overline{v}}_2{D}_2\left(\tau ;\varphi \right)+{\lambda }_t{E}_1\left(\tau ;\varphi \right)+{\alpha }_t{E}_2\left(\tau ;\varphi \right)\\ +\frac{1}{2}{\epsilon }_t^2{D}_1{\left(\tau ;\varphi \right)}^2+\frac{1}{2}{\delta }_t^2{D}_2{\left(\tau ;\varphi \right)}^2+\frac{1}{2}{\eta }_t^2{E}_1{\left(\tau ;\varphi \right)}^2+\frac{1}{2}{\beta }_t^2{E}_2{\left(\tau ;\varphi \right)}^2.\end{array}$ (10)

The initial conditions are $C\left(0;\varphi \right)=D_1\left(0;\varphi \right)=D_2\left(0;\varphi \right)=E_1\left(0;\varphi \right)=E_2\left(0;\varphi \right)=0$ . Similar to the results of previous researchers (Heston, 1993; Christoffersen et al., 2009; Lin & He, 2021), it is not difficult to solve the analytical solutions of the Riccati equations $D_1\left(\tau ;\varphi \right)$ and $D_2\left(\tau ;\varphi \right)$ . Subsequently, $C\left(\tau ;\varphi \right)$ , $E_1\left(\tau ;\varphi \right)$ , and $E_2\left(\tau ;\varphi \right)$ can be integrated from the remaining ODEs. Specifically, $C\left(\tau ;\varphi \right)$ can be solved as

$\begin{array}{c}C\left(\tau ;\varphi \right)=\overline{C}\left(\tau ;\varphi \right)+{\displaystyle {\int }_0^{\tau }\langle {\lambda }_t{E}_1\left(s;\varphi \right)+{\alpha }_t{E}_2\left(s;\varphi \right)+\frac{1}{2}{\epsilon }_t^2{D}_1^2\left(s;\varphi \right)}\\ +\frac{1}{2}{\delta }_t^2{D}_2^2\left(s;\varphi \right)+\frac{1}{2}{\eta }_t^2{E}_1^2\left(s;\varphi \right)+\frac{1}{2}{\beta }_t^2{E}_2^2\left(s;\varphi \right),X_s\rangle \text{d}s,\end{array}$ (11)

where $\overline{C}\left(\tau ;\varphi \right)$ is defined as

$\overline{C}\left(\tau ;\varphi \right)=ri\varphi \tau +{\overline{v}}_1{E}_1\left(\tau ;\varphi \right)+{\overline{v}}_2{E}_2\left(\tau ;\varphi \right).$

Now that we have derived the internal expected value, the next step is to compute the external expectation in order to obtain the characteristic function.

$\begin{array}{l}f\left(\varphi ;\tau ,y,v_1,v_2,{\theta }_1,{\theta }_2,X_t\right)\\ =\mathbb{E}\left[h\left(\varphi ;T,y,v_1,v_2,{\theta }_1,{\theta }_2,X_t|X_T\right)|y_t,v_{1t},v_{2t},{\theta }_{1t},{\theta }_{2t},X_t\right]\\ ={\text{e}}^{C\left(T;\varphi \right)+D_1\left(T;\varphi \right)v_1+D_2\left(T;\varphi \right)v_2+E_1\left(T;\varphi \right){\theta }_1+E_2\left(T;\varphi \right){\theta }_2+i\varphi y}\cdot \mathbb{E}\left[{\text{e}}^{{\displaystyle {\int }_0^{\tau }\Lambda \left(s;\varphi ,X_s\right)\text{d}s}}|X_t\right].\end{array}$ (12)

According to Elliott and Lian (Elliott & Lian, 2013), this expectation can be further derived as

$\mathbb{E}\left[{\text{e}}^{{\displaystyle {\int }_0^{\tau }\Lambda \left(t;\varphi ,X_s\right)\text{d}s}}|X_t\right]=\langle {\text{e}}^M{X}_t,I\rangle ,$ (13)

where

$\begin{array}{c}\Lambda \left(s;\varphi ,X_s\right)=\langle {\lambda }_t{E}_1\left(s;\varphi \right)+{\alpha }_t{E}_2\left(s;\varphi \right)+\frac{1}{2}{\epsilon }_t^2{D}_1^2\left(s;\varphi \right)+\frac{1}{2}{\delta }_t^2{D}_2^2\left(s;\varphi \right)\\ +\frac{1}{2}{\eta }_t^2{E}_1^2\left(s;\varphi \right)+\frac{1}{2}{\beta }_t^2{E}_2^2\left(s;\varphi \right),X_s\rangle ,\end{array}$

and $M=A^{\text{T}}\tau +B$ . Here, $A=\left(\begin{array}{cc}-{\lambda }_{12}& {\lambda }_{12}\\ {\lambda }_{21}& -{\lambda }_{21}\end{array}\right)$ represents the transition rate matrix of the Markov chain $X_t$ and $B=\left(\begin{array}{cc}{p}_1& 0\\ 0& p_2\end{array}\right)$ with $p_1={\displaystyle {\int }_0^{\tau }{\Lambda }_1\left(s;\varphi \right)\text{d}s}$ , $p_2={\displaystyle {\int }_0^{\tau }{\Lambda }_2\left(s;\varphi \right)\text{d}s}$ . The functions ${\Lambda }_1\left(s;\varphi \right)$ and ${\Lambda }_2\left(s;\varphi \right)$ are represented as follows.

$\begin{array}{c}{\Lambda }_1\left(s;\varphi \right)={\lambda }_1{E}_1\left(s;\varphi \right)+{\alpha }_1{E}_2\left(s;\varphi \right)+\frac{1}{2}({\epsilon }_1^2{D}_1^2\left(s;\varphi \right)+{\delta }_1^2{D}_2^2\left(s;\varphi \right)\\ +{\eta }_1^2{E}_1^2\left(s;\varphi \right)+{\beta }_1^2{E}_2^2\left(s;\varphi \right)),\end{array}$

$\begin{array}{c}{\Lambda }_2\left(s;\varphi \right)={\lambda }_2{E}_1\left(s;\varphi \right)+{\alpha }_2{E}_2\left(s;\varphi \right)+\frac{1}{2}({\epsilon }_2^2{D}_1^2\left(s;\varphi \right)+{\delta }_2^2{D}_2^2\left(s;\varphi \right)\\ +{\eta }_2^2{E}_1^2\left(s;\varphi \right)+{\beta }_2^2{E}_2^2\left(s;\varphi \right)).\end{array}$

After deriving the closed-form solution for the new model, we will investigate its properties via comprehensive numerical experiments. Prior to this, it is crucial to conduct a numerical comparison between our analytical results and those obtained from Monte Carlo simulations. This step ensures the absence of algebraic errors and validates the robustness of our conclusions. Additionally, the following section will delve into the implications of integrating two long-term mean-reverting processes into the double Heston model, along with the introduction of a transformation mechanism.

3. Numerical Experiments and Discussions

This section examines various properties of our newly derived formulas through numerical experiments. Initially, we compare European option prices calculated using Monte Carlo simulations with those from our formula to verify the pricing accuracy of the MRSNDH model and the NDH model. Next, by comparing option prices between the NDH model and the double Heston model, we evaluate the effect of including a stochastic long-term mean process. Finally, by contrasting the MRSNDH model with the NDH model, we explore how the state transition model influences option prices.

In this section, the initial parameters are set as follows: the current state is set to 1, the risk-free rate ($r$ ) is 0.05, the mean reversion speeds ($k_1$ and $k_2$ ) are 6 and 8, respectively, the constant parts of the long-term mean (${\overline{v}}_1$ and ${\overline{v}}_2$ ) are specified, the initial values of the volatility process ($v_1$ and $v_2$ ) are 0.2 and 0.1, respectively, the initial values of the stochastic long-term mean (${\theta }_1$ and ${\theta }_2$ ) are 0.12 and 0.08, respectively, the volatility terms (${\sigma }_1$ and ${\sigma }_2$ ) for asset volatility are 0.1 and 0.2, respectively, the correlation coefficients (${\rho }_1$ and ${\rho }_2$ ) are 0.1 and 0.2, respectively, the two transition rates are set to (${\lambda }_{12}=10$ and ${\lambda }_{21}=20$ ), the option strike price ($K$ ) is 10, and the parameter values for ($X_t$ ) related to the Markov chain will be explained in the figure’s title. Consistent with the standard Monte Carlo simulation framework, we generate 500,000 independent sample paths; each path yields a single option price estimate, and the final Monte Carlo estimator is obtained as the arithmetic mean of these 500,000 simulated prices.

The MRSNDH model’s pricing formula, under initial parameters, yields prices closely matching those from Monte Carlo simulations, as illustrated in Figure 1. Specifically, Figure 1(a) shows that both prices rise with increasing underlying asset value, aligning with expected financial behaviour. Figure 1(b) indicates that the relative error between the two methods is below 0.6%. Similarly, Figure 2(a) compares option prices derived from our formula with those from Monte Carlo simulations under the NDH model.

Note: Time-related parameters have been set to: $\tau =1$ , ${\epsilon }_1=0.01$ , ${\epsilon }_2=0.02$ ; ${\delta }_1=0.015$ , and ${\delta }_2=0.01$ ; ${\lambda }_1=0.03$ , ${\lambda }_2=0.01$ ; ${\eta }_1=0.01$ , ${\eta }_2=0.015$ ; ${\alpha }_1=0.01$ , ${\alpha }_2=0.02$ ; ${\beta }_1=0.01$ , ${\beta }_2=0.015$ .

Figure 1. Verification of MRSNDH model accuracy against the benchmark: (a) Absolute price comparison, (b) Relative error.

When the MRSNDH model is simplified to the NDH model, variables $\epsilon$ , $\delta$ , $\lambda$ , $\eta$ , $\alpha$ and $\beta$ assume the values associated with state 1 in Figure 1. Figure 2(a) shows that the prices are very similar, and the price of the European call option increases monotonically with the underlying asset’s value. Figure 2(b) indicates that the relative error between the two does not exceed 0.3%. Figure 1 and Figure 2 demonstrate that the formulas derived in this chapter are accurate, regardless of whether a state transition mechanism is included in the model.

Note: The parameters involved in the NDH model are: $\tau =1$ , $\epsilon =0.01$ , $\delta =0.015$ , $\lambda =0.03$ , $\eta =0.01$ , $\alpha =0.01$ , $\beta =0.01$ .

Figure 2. Verification of the NDH model accuracy against the benchmark: (a) Absolute price comparison, (b) Relative error.

Note: The relevant parameters at this time are: $\tau =1$ , $\epsilon ={\epsilon }_1=0.01$ , ${\epsilon }_2=0.02$ ; $\delta ={\delta }_1=0.015$ , ${\delta }_2=0.01$ ; $\lambda ={\lambda }_1=0.03$ , ${\lambda }_2=0.01$ ; $\eta ={\eta }_1=0.01$ , ${\eta }_2=0.015$ ; $\alpha ={\alpha }_1=0.01$ , ${\alpha }_2=0.02$ ; $\beta ={\beta }_1=0.01$ , ${\beta }_2=0.015$ ; ${\lambda }_{12}=10z$ , ${\lambda }_{21}=20z$ .

Figure 3. MRSNDH model versus NDH model with respect to a scale parameter z.

After the accuracy of the pricing formula has been verified, the impact of introducing the state transition model can be demonstrated. A scaling parameter $z$ is introduced to adjust the transition rates ${\lambda }_{12}=10$ and ${\lambda }_{21}=20$ , as depicted in Figure 3. When $z$ equals 0, due to the lack of state transitions, the model with state transition mechanisms simplifies to one without state transitions, leading to identical option prices between the MRSNDH model and the NDH model. It is important to note that under current parameters, as the scaling parameter $z$ increases, the European call option price of the MRSNDH model decreases, with the rate of decrease slowing down.

In Figure 4 that follows, the European call option prices of the NDH model and the double Heston model are compared. Initially, their prices are similar; however, as the time to maturity increases, the difference between them becomes more pronounced, with the option prices from the NDH model consistently higher than those from the double Heston model.

Figure 4. Comparison of NDH model and the double Heston model option prices at different expiration times. Where, $\epsilon =0.0\text{1}$ ; $\delta =0.0\text{15}$ , ${\delta }_2=0.0\text{1}$ ; $\lambda =0.03$ ; $\eta =0.0\text{1}$ ; $\alpha =0.0\text{1}$ ; $\beta =0.0\text{1}$ .

Figure 5 illustrates the impact of introducing a state transition mechanism into the NDH model. Option prices in the NDH model are consistently higher than those in the MRSNDH model. To investigate potential parameter reversal effects between the MRSNDH and NDH models, we compared versions with reversed parameters to their original counterparts. Our analysis shows that parameter reversal does not cause option prices in the MRSNDH model to exceed those in the NDH model; instead, it reduces them further. Note that as options approach expiration, prices converge across all three models due to decreasing likelihood of state transitions and significant price fluctuations. Conversely, differences between the models become more pronounced with longer time to expiration.

Figure 5. The comparison of option prices with different expiration dates between the NDH model after parameter update and the MRSNDH model. The reset parameters are as follows: ${\sigma }_1=0.5$ , ${\sigma }_2=0.3$ ; $\lambda ={\lambda }_1=0.07$ , ${\lambda }_2=0.01$ ; $\eta ={\eta }_1=0.025$ , ${\eta }_2=0.005$ ; $\alpha ={\alpha }_1=0.035$ , ${\alpha }_2=0.01$ ; $\beta ={\beta }_1=0.005$ , ${\beta }_2=0.03$ ; $\epsilon ={\epsilon }_1=0.001$ , ${\epsilon }_2=0.002$ ; $\delta ={\delta }_1=0.002$ , ${\delta }_2=0.001$ . When the parameters are reversed, only the parameters with Markov state transitions will be swapped between 1 and 2.

4. Empirical Studies

In this section, we conduct an empirical study to assess the performance of the proposed model relative to both the double Heston model (Christoffersen et al., 2009) and the He-Chen model (He & Chen, 2021). The purpose is to demonstrate the importance of integrating a series of regime-switching factors and stochastic long-term means, as well as additional stochastic volatility components, within the double Heston framework and to assess the performance of this new model when applied to real financial data. We first describe the dataset and the key screening conditions used in model calibration, and then detail the parameter estimation method. Finally, we present the empirical results, which clearly illustrate the relative performance of the three models in fitting real financial data.

4.1. Data Description

A dataset of XSP European call options on the S&P 500 Index is used for empirical analysis, covering the period from July to September 2024. As is standard practice, the mid-price, calculated as the average of the bid and ask prices, was used as the proxy for the option price in this study. However, it is important to note that raw data cannot be directly used in parameter estimation due to the presence of market noise. To address this issue, appropriate filtering techniques were applied to preprocess the data prior to use in the calibration process. To characterize market regimes, this study classifies trading days in the sample period (July-September 2024) into two distinct states—“high-volatility regime” (state 1) and “low-volatility regime” (state 2)—based on implied volatility (IV). Specifically, the sample-wide mean IV serves as the classification threshold: days with IV above this mean are assigned to state 1; all others to state 2.

First, only option prices observed on Wednesdays and Thursdays are used in the analysis. Specifically, Wednesday observations are employed for parameter estimation, while Thursday observations serve as market benchmarks against which the theoretical option prices—derived from the estimated parameters—are evaluated. This approach is a well-established practice in model calibration and is justified on two key grounds. On one hand, using only one trading day per week mitigates temporal dependence between calibration windows and yields a longer, more statistically independent time series—thereby enhancing the reliability of the results, especially given the computationally intensive nature of the calibration process. On the other hand, Wednesdays and Thursdays are less likely to coincide with U.S. public holidays and exhibit weaker “day-of-the-week” effects than Mondays and Fridays. Second, options with time to maturity shorter than 30 days or longer than 120 days are excluded from the sample. The former typically exhibits low time value and high bid-ask spreads, which may distort calibration outcomes (Bakshi et al., 1997; He & Chen, 2021). Third, options with absolute moneyness exceeding 10% are also excluded from the sample (Shu & Zhang, 2004). In this way, deep in-the-money and deep out-of-the-money options are excluded due to potential liquidity issues and heightened sensitivity to model misspecification. Note that absolute moneyness is defined as the absolute relative difference between the S&P 500 Index level and the corresponding strike price, i.e., $\text{moneyness}=\frac{\left|S-K\right|}{K}$ . In addition to careful selection of the option data, the risk-free interest rate must be determined in advance. The risk-neutral measure—grounded in the no-arbitrage principle and implemented via equivalent martingale measure transformation—converts pricing under heterogeneous investor risk preferences into a tractable expectation computation discounted at the observable risk-free rate. As such, it serves as a foundational tool in derivative valuation, risk management, and empirical analysis of market behavior. In this empirical study, the framework is rigorously applied: the risk-free rate is treated as a core structural parameter, explicitly estimated within the model’s parameter set and consistently employed throughout calibration, derivative pricing, and statistical inference. Here, we use the daily three-month U.S. Treasury bill yield as a proxy for the risk-free interest rate. This maturity is appropriate because the maximum time to maturity among the selected options is less than 120 days (He & Lin, 2023). Following rigorous data screening, 250 high-quality, analytically suitable observations remain for subsequent parameter estimation and predictive modeling. With all required data available, parameter estimation is performed via numerical optimisation; details are provided in the next subsection.

4.2. Parameter Estimation

In this section, we will first review the parameters of all models that need to be determined, and then introduce a specific global optimisation method to obtain all model parameters.

Recall the dynamics of the double Heston model as presented. It is evident that ten parameters require estimation: the mean-reversion speed (${\kappa }_1$ and ${\kappa }_2$ ), the long-term average level (${\theta }_1$ and ${\theta }_2$ ), the volatility of volatility (${\sigma }_1$ and ${\sigma }_2$ ), the correlation factor (${\rho }_1$ and ${\rho }_2$ ), and the initial value of volatility ($v_1$ and $v_2$ ). In contrast, the dynamics of our newly proposed model indicate that 26 parameters must be determined before model evaluation. Eight are the same as those in the double Heston model, including ${\kappa }_i$ , ${\sigma }_i$ , ${\rho }_i$ , $v_i$ , while the stochastic process governing the long-term mean ${\theta }_i$ and the constant parts of the long-term mean ${\overline{v}}_i$ , $i=\left\{1,2\right\}$ , the remaining six parameters, i.e., $\lambda$ , $\eta$ , $\alpha$ , $\beta$ , $\epsilon$ and $\delta$ . During parameter calibration, we enforced the Feller condition $2\kappa \theta \ge {\sigma }^2$ for both variance processes to ensure square-root diffusion admissibility under the risk-neutral measure, and further constrained all parameters to empirically plausible ranges—consistent with prior literature and observed market dynamics.

Having known all the needed parameters, we can now proceed to the estimation part. One of the most popular methods in determining model parameters is to find a set of “optimal” need to choose an appropriate definition for such a distance. In fact, a common approach is to take the Relative Mean Squared Error (RMSE). However, drawing on the insights from case (He & Chen, 2021), it is evident that the aforementioned error calculation method is not suitable for the research presented in this paper. Alternatively, the Mean Squared Error (MSE) method should be employed to assess the estimation results of the parameters.

$\text{MSE}=\frac{1}{N}{\displaystyle \sum }_{i=1}^N{\left(C_i^{\text{Market}}-C_i^{\text{Model}}\right)}^2,$ (14)

As the objective function to measure the distance, with $C^{\text{Market}}$ and $C^{\text{Model}}$ being the market price of an option and the same option calculated from our pricing formula with a particular set of parameters, respectively. $N$ is the total number of observations selected in a single estimation.

Another issue is to choose an appropriate method to minimise the selected objective function, which is a minimisation problem. In the literature, local minimisation is a first choice as it is easy to implement and fast to produce a result. Unfortunately, the objective function is not necessarily convex and thus there exist several local minima. An appropriate initial guess of the solution is usually very crucial for the local minimisation method to be safely used, as it would otherwise be easily stuck in a local minimum and produce unreliable results. Therefore, in this case, global optimisation is much favoured because a properly designed global optimisation algorithm is able to skip local minima and correctly identify the global minimum in an efficient way.

There are various methods for global optimisation; however, due to our model requiring the estimation of numerous parameters, efficiently estimating these parameters under limited computational resources is challenging. PSO is a population-based intelligent optimisation algorithm that offers simplicity, fast convergence, and strong global search capabilities. APSO is an improved version of PSO that incorporates adaptive mechanisms. By introducing adaptive inertia weights and learning factors, APSO dynamically adjusts its search capability, thereby enhancing convergence speed and accuracy. Therefore, we employed APSO for parameter estimation, with the results presented in Table 1. This time, the Hen-Chen model and the double Heston model were selected for comparative experiments with the MRSNDH model. Using APSO, parameter estimation was carried out based on the actual data set, and the obtained results are shown in Table 1. With all the estimated parameters available, we are now able to assess the performance of our newly proposed model. This will be the main issue of the next subsection.

Table 1. Parameter estimation results for three models.

Parameters	MRSNDH Model	He-Chen Model	Double Heston Model
$V_1$	0.0101	0.0100	0.0100
$V_2$	0.0312	-	0.0100
${\rho }_1$	−0.7712	−0.8990	−0.7942
${\rho }_2$	−0.8296	-	−0.3095
${\kappa }_1$	1.4116	0.5000	8.6853
${\kappa }_2$	9.6448	-	7.3456
${\sigma }_1$	0.4687	0.1762	0.1131
${\sigma }_2$	0.3896	-	0.1000
${\overline{v}}_1$	0.0466	-	-
${\overline{v}}_2$	0.0100	-	-
${\theta }_1$	0.0794	0.3182	0.0500
${\theta }_2$	0.0522	-	0.5000
${\lambda }_1$	0.0294	0.1281	-
${\lambda }_2$	−0.1520	-	-
${\eta }_1$	0.0184	0.1188	-
${\eta }_2$	0.1908	-	-
${\alpha }_1$	0.0283	-	-
${\alpha }_2$	0.2388	-	-
${\beta }_1$	0.0217	-	-
${\beta }_2$	0.0134	-	-
${ϵ}_1$	0.2429	-	-
${ϵ}_2$	0.1920	-	-
${\overline{\delta }}_1$	0.1567	-	-
${\overline{\delta }}_2$	0.0490	-	-
${\lambda }_{12}$	11.7251	-	-
${\lambda }_{21}$	2.3301	-	-

4.3. Empirical Results

In this section, we present the empirical results of our newly proposed model alongside those of the double Heston model and He-Chen model, based on the same set of option data. In particular, Table 2 exhibits the in- and out-of-sample errors of the three models. It can be seen from Table 2 that our newly proposed model is significantly superior to the other two models in terms of both in-sample and out-of-sample errors, while the error results of the He-Chen model and the double Heston model are similar both out-of-sample and in-sample. Specifically, from the perspective of in-sample error, the daily average MSE of our model is 45.4209, which is only approximately 40% of that of the double Heston model. On the other hand, when considering the out-of-sample error, the errors of the three models are all greater than their in-sample errors. This is in line with financial intuition. After all, in actual predictions, the out-of-sample error is mostly higher than the in-sample error. Furthermore, the out-of-sample error difference between our model and the other two models has further widened. The MSE of our model is approximately 36% of that of the He-Chen model, which also indicates that our model is superior to the He-Chen model. If the in-sample and out-of-sample errors of one model are both lower than those of another model, it can be considered that the model is better. Therefore, it can be concluded that our newly proposed model is at least a better choice for the selected dataset.

Table 2. Comparison of in-sample and out-of-sample errors.

Model	In-sample error	Out-of-sample error
MRSNDH	45.4209	46.1136
He-Chen	112.8932	127.1055
Double Heston	111.2830	125.3227

A further point of interest is the comparative performance of the three models across varying levels of moneyness. Accordingly, we present out-of-sample pricing errors disaggregated by moneyness category—namely, in-the-money (0.90 < S/K < 0.97), at-the-money (0.97 ≤ S/K ≤ 1.03), and out-of-the-money (1.03 < S/K < 1.10)—with mean absolute errors and root-mean-square errors reported in Table 3. It can be seen from Table 3 that the MRSNDH model is at least a better choice than the double Heston model and the He-Chen model on the adopted datasets. In particular, it can be clearly observed that the out-of-sample error of out-of-the-money options is much larger than that of in-the-money options and at-the-money options. For this category, our model performs better than the Heston model. The improvements of the other two types are also more significant. Compared with the other two models, the error of our model in at-the-money options has decreased by nearly 66%, while in in-the-money options, the error gap is the most obvious. The error of the MRSNDH model is even less than 50% of the other two models. Therefore, we can confidently draw the conclusion that our model can certainly be a good competitor of the Heston model in the actual market.

Table 3. Comparison of different option errors.

Model	In the money	At the money	Out of money
MRSNDH	29.2542	48.3971	117.0683
He-Chen	60.7964	143.1054	168.2299
Double Heston	59.4057	141.1682	168.2299

5. Conclusion

In this paper, we propose the MRSNDH model by incorporating an additional volatility term into the volatility process and modeling the long-term mean of the double Heston framework through another stochastic process. This structure enhances the model’s ability to fit real-world financial data. After deriving a closed-form pricing formula for European options under the MRSNDH model, we numerically validate its accuracy by comparing the results with those obtained from Monte Carlo simulations. Furthermore, we conduct a comparative numerical analysis of option prices generated by the Heston model and our proposed model to highlight their differences. Finally, an empirical study based on S&P 500 index options demonstrates that the MRSNDH model consistently outperforms both the double Heston model and the He-Chen model, particularly for in-the-money options. These findings suggest that the MRSNDH model can serve as a more effective alternative to the double Heston model in practical financial applications. Future research may explore extending the model to other types of derivative instruments and investigating alternative state-switching mechanisms.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Bakshi, G., Cao, C., & Chen, Z. (1997). Empirical Performance of Alternative Option Pricing Models. The Journal of Finance, 52, 2003-2049. [Google Scholar] [CrossRef]
[2]	Bakshi, G., Ju, N., & Ou-Yang, H. (2006). Estimation of Continuous-Time Models with an Application to Equity Volatility Dynamics. Journal of Financial Economics, 82, 227-249. [Google Scholar] [CrossRef]
[3]	Beckers, S. (1983). Variances of Security Price Returns Based on High, Low, and Closing Prices. The Journal of Business, 56, 97-112. [Google Scholar] [CrossRef]
[4]	Black, F., & Scholes, M. (1973). The Pricing of Options and Corporate Liabilities. Journal of Political Economy, 81, 637-654. [Google Scholar] [CrossRef]
[5]	Christoffersen, P., Heston, S., & Jacobs, K. (2009). The Shape and Term Structure of the Index Option Smirk: Why Multifactor Stochastic Volatility Models Work So Well. Management Science, 55, 1914-1932. [Google Scholar] [CrossRef]
[6]	Dumas, B., Fleming, J., & Whaley, R. E. (1998). Implied Volatility Functions: Empirical Tests. The Journal of Finance, 53, 2059-2106. [Google Scholar] [CrossRef]
[7]	Dupire, B. (1994). Pricing with a Smile. Risk, 7, 18-20.
[8]	Elliott, R. J., & Lian, G. (2013). Pricing Variance and Volatility Swaps in a Stochastic Volatility Model with Regime Switching: Discrete Observations Case. Quantitative Finance, 13, 687-698. [Google Scholar] [CrossRef]
[9]	Eraker, B. (2004). Do Stock Prices and Volatility Jump? Reconciling Evidence from Spot and Option Prices. The Journal of Finance, 59, 1367-1403. [Google Scholar] [CrossRef]
[10]	Gong, X., & Zhuang, X. (2016). Option Pricing and Hedging for Optimized Lévy Driven Stochastic Volatility Models. Chaos, Solitons & Fractals, 91, 118-127. [Google Scholar] [CrossRef]
[11]	Hagan, P. S., Kumar, D., Lesniewski, A. S., & Woodward, D. E. (2002). Managing Smile Risk. Wilmott Magazine, 1, 84-108.
[12]	Hamilton, J. D. (1990). Analysis of Time Series Subject to Changes in Regime. Journal of Econometrics, 45, 39-70. [Google Scholar] [CrossRef]
[13]	He, G. (2017). Option Volatility Estimation Based on Particle Swarm Optimization Algorithm. Journal of Sichuan University (Natural Science Edition), 54, 925-928. (In Chinese)
[14]	He, X., & Chen, W. (2021). A Closed-Form Pricing Formula for European Options under a New Stochastic Volatility Model with a Stochastic Long-Term Mean. Mathematics and Financial Economics, 15, 381-396. [Google Scholar] [CrossRef]
[15]	He, X., & Lin, S. (2023). A Closed-Form Pricing Formula for European Options under a New Three-Factor Stochastic Volatility Model with Regime Switching. Japan Journal of Industrial and Applied Mathematics, 40, 525-536. [Google Scholar] [CrossRef]
[16]	Heston, S. L. (1993). A Closed-Form Solution for Options with Stochastic Volatility with Applications to Bond and Currency Options. Review of Financial Studies, 6, 327-343. [Google Scholar] [CrossRef]
[17]	Hull, J., & White, A. (1987). The Pricing of Options on Assets with Stochastic Volatilities. The Journal of Finance, 42, 281-300. [Google Scholar] [CrossRef]
[18]	Jacquier, E., Polson, N. G., & Rossi, P. E. (2004). Bayesian Analysis of Stochastic Volatility Models with Fat-Tails and Correlated Errors. Journal of Econometrics, 122, 185-212. [Google Scholar] [CrossRef]
[19]	Kennedy, J., & Eberhart, R. (1995). Particle Swarm Optimization. In Proceedings of ICNN’95—International Conference on Neural Networks (pp. 1942-1948). IEEE. [Google Scholar] [CrossRef]
[20]	Lin, S., & He, X. (2021). Analytically Pricing European Options under a New Two-Factor Heston Model with Regime Switching. Computational Economics, 59, 1069-1085. [Google Scholar] [CrossRef]
[21]	Mehrdoust, F., Noorani, I., & Hamdi, A. (2022). Calibration of the Double Heston Model and an Analytical Formula in Pricing American Put Option. Journal of Computational and Applied Mathematics, 392, Article ID: 113422. [Google Scholar] [CrossRef]
[22]	Mehrdoust, F., Noorani, I., & Hamdi, A. (2023). Two-Factor Heston Model Equipped with Regime-Switching: American Option Pricing and Model Calibration by Levenberg-Marquardt Optimization Algorithm. Mathematics and Computers in Simulation, 204, 660-678. [Google Scholar] [CrossRef]
[23]	Ratnaweera, A., Halgamuge, S. K., & Watson, H. C. (2004). Self-Organizing Hierarchical Particle Swarm Optimizer with Time-Varying Acceleration Coefficients. IEEE Transactions on Evolutionary Computation, 8, 240-255. [Google Scholar] [CrossRef]
[24]	Scott, L. O. (1987). Option Pricing When the Variance Changes Randomly: Theory, Estimation, and an Application. The Journal of Financial and Quantitative Analysis, 22, 419-438. [Google Scholar] [CrossRef]
[25]	Shi, Y., & Eberhart, R. (1998). A Modified Particle Swarm Optimizer. In 1998 IEEE International Conference on Evolutionary Computation Proceedings. IEEE World Congress on Computational Intelligence (Cat. No.98TH8360) (pp. 69-73). IEEE. [Google Scholar] [CrossRef]
[26]	Shu, J., & Zhang, J. E. (2004). Pricing S&P 500 Index Options under Stochastic Volatility with the Indirect Inference Method. Journal of Derivatives Accounting, 1, 171-186. [Google Scholar] [CrossRef]
[27]	Wiggins, J. B. (1987). Option Values under Stochastic Volatility: Theory and Empirical Estimates. Journal of Financial Economics, 19, 351-372. [Google Scholar] [CrossRef]
[28]	Zhang, S. M., & Feng, Y. (2019). American Option Pricing under the Double Heston Model Based on Asymptotic Expansion. Quantitative Finance, 19, 211-226. [Google Scholar] [CrossRef]