1. Introduction
A derivative contract is a financial instrument whose value is determined by the performance of an underlying asset. Options are a specific type of derivative that grants the holder the right, but not the obligation, to purchase or sell a fixed quantity of the underlying asset. The underlying asset may be bought (a call option) or sold (a put option) at a strike price and within a specified time (expiry time).
While the strike price and expiry time are predetermined in vanilla options, exotic options may be based on multiple underlying assets, the strike price may depend on several market factors, and numerous exercise dates may be possible. The two most prominent options are European and American, each with distinct exercise time conditions. There are many more options, including Asian and Bermudan options.
European options can only be exercised when the options contract expires, which is predetermined when the contract is issued. On the other hand, American options can be exercised on any trading day before the options contract expires. Options pricing has been a significant topic of theoretical exploration and practical implementation.
Fischer Black, Myron Scholes, and Robert Merton made a significant contribution to this field. Black and Scholes [Reference Black and Scholes1] proposed a model to price the European options with the assumption that the underlying asset follows the geometric Brownian motion. The analytical formula is derived for European option prices under this model. Despite its widespread use, the simplified assumptions made for analytical simplicity and tractability may not fully capture real-world complexities, potentially leading to some degree of mispricing.
Numerous empirical studies have shown that volatility clustering is a distinctive feature of the equity return rates [Reference Fama8]. The distribution of equity return rates exhibits high peaks and fat tails [Reference Cont2]. A negative correlation exists between share prices and volatility [Reference Engle and Ng6]. Lastly, mean reversion in return rates is observed [Reference Engle and Patton7]. The Black–Scholes model does not account for these additional characteristics.
Two types of nonconstant volatility models are widely studied in the literature to relax the assumption of constant volatility: local volatility and stochastic volatility. In the local volatility model, volatility is considered a deterministic function of the underlying asset price and time. In the stochastic volatility models, volatility is considered a different stochastic process, allowing for dynamic changes. Numerous authors have addressed the options pricing problem under stochastic volatility models. However, finding a closed-form pricing formula for European options becomes complicated after introducing an additional stochastic process, and often numerical methods are implemented for existing models, such as Hull and White [Reference Hull and White23].
Heston [Reference Heston22] proposed a stochastic volatility model where the volatility process follows a square root mean-reverting process, ensuring non-negativity. The Heston model accurately captures the asset dynamics, assuming a constant interest rate and the absence of sudden market shocks. It demonstrates the implied volatility smile effect, a well-known inconsistency in the Black–Scholes model. Recently, He and Chen [Reference He and Chen15] introduced a novel stochastic volatility model by incorporating a stochastic process in the long-term mean of volatility within the Heston framework. Lin and He [Reference Lin and He26] proposed an additional factor in the stochastic volatility process, which follows regime-switching. Pasricha and Goel [Reference Pasricha and Goel34] considered the Heston model for pricing exchange options. He and Lin [Reference He and Lin18] considered the long-term mean of the volatility process to follow a regime-switching stochastic process.
The assumption of a constant interest rate is not realistic, as the real market experiences fluctuations in interest rates. Therefore, the stochastic behavior of interest rates is important in modeling asset prices. To address this, various short-rate models have been introduced in the literature to accurately capture the term structure of interest rates. Vasicek [Reference Vasicek37], Cox et al. [Reference Cox, Ingersoll and Ross3], and Hull and White [Reference Hull and White24] are some examples.
Hamilton [Reference Hamilton13] played a pivotal role in advancing the regime-switching models by proposing a class of discrete-time, Markov regime-switching, and autoregressive time-series models. Empirically, Markov regime-switching models can describe various essential characteristics of economic and financial time series, such as the asymmetry and heavy-tailedness of asset returns, time-varying conditional volatility, volatility clustering, regime-switching, nonlinearity, and other complex features.
Markov regime-switching models offer a framework for identifying variations in economic conditions. Elliott and Nishide [Reference Elliott and Nishide4] derived a closed-form expression for bond prices by considering that the short rate follows the Cox–Ingersoll–Ross (CIR) model with a Markov regime-switching. Elliott et al. [Reference Elliott, Nishide and Osakwe5] extended their work by incorporating regime-switching into the Heston stochastic volatility model and derived the price of European options. Shen and Siu [Reference Shen and Siu36] derived bond options price under the Markov regime-switching Hull–White model.
To account for the fluctuations caused by external factors, Merton [Reference Merton32] introduced the jumps in the diffusion model, and the options price is derived. Additionally, Kou [Reference Kou25] studied the jump-diffusion model, particularly when the jump size follows an asymmetric double exponential distribution. For further applications of the jump process in modeling stock price process, refer, Wang et al. [Reference Wang, Wang and Liu38], Han [Reference Han14], Guo and Bai [Reference Guo and Bai12], and Lin et al. [Reference Lin, Miao, Lee and Zheng29].
Grzelak et al. [Reference Grzelak, Oosterlee and Van Weeren10] considered a hybrid stochastic volatility and stochastic interest rate model for pricing European options. Guo [Reference Guo11] investigated the pricing of the European options using the Heston–Vasicek model. The hybrid stochastic rate and stochastic volatility model are considered to price various financial instruments (see, e.g., Grzelak and Oosterlee [Reference Grzelak and Oosterlee9], Wu et al. [Reference Wu, Jia, Yang and Liu39]). lyu et al. [Reference lyu, Ma and Sun30] evaluated option prices within the Markov regime-switching double Heston stochastic volatility model. The model included a stochastic interest rate and Poisson jumps. Shan et al. [Reference Shan, Yi, Zhang and Shu35] determined options price when the Markov-modulated Merton jump-diffusion process describes the discrete dividends. Ma et al. [Reference Ma, Chen and Lyu31] priced the European options by incorporating the Heston volatility within a Markov regime-switching framework. The model considered the jump process to follow a Poisson process with stochastic intensity. The CIR model is used to model the interest rate, and the size of the jumps follows a double exponential distribution.
Lin and He [Reference Lin and He27] considered the CIR volatility process having stochastic equilibrium levels with regime-switching. He and Lin [Reference He and Lin17] extended the stochastic volatility model [Reference He and Chen15] by considering regime-switching stochastic volatility. He et al. [Reference He, Pasricha and Lin21] considered the regime-switching volatility and stochastic market liquidity risk for pricing European options. Lin and He [Reference Lin and He28] considered the CIR stochastic volatility with market liquidity risks for pricing variance and volatility swaps. He and Lin [Reference He and Lin20] priced the European options considering Heston stochastic volatility and market liquidity risks. He and Lin [Reference He and Lin19] proposed a stochastic interest rate model with long-run mean as another stochastic process. He and Lin [Reference He and Lin16] considered the correlation between asset price and risk-free rate using two-factor hybrid Heston Hull–White model and derived the European options price.
While the CIR model ensures positive interest rates, it does not align well with the negative interest rates seen during the 2007 financial crisis and subsequent quantitative easing policies. For more details on the significance of negative interest rates and drawbacks of the CIR model during low interest rates, refer Orlando and Bufalo [Reference Orlando and Bufalo33]. Therefore, the market requires a model that can deal with negative interest rates. Vasicek model is one such model, and the Hull–White model generalizes it, giving a more general framework for interest rate dynamics. Hence, instead of the CIR model, the Hull–White model is considered.
In this article, the pricing of European options within the framework of a Heston stochastic volatility integrated with a Hull–White stochastic interest rate model is discussed. The long-run mean of the volatility process, the mean reversion level, and the volatility of the risk-free rate alternate between a finite number of regimes. These regimes correspond to the finite state space of a continuous-time Markov chain. The jump occurs based on a Poisson process with stochastic intensity, and the jump sizes follow an asymmetric double exponential distribution. To demonstrate the accuracy of the derived formula, the prices are computed numerically and compared with Monte Carlo simulation (MCS) results. The comparison of this model to the previous models and the sensitivity of price to various parameters are also demonstrated.
The paper is structured as follows: Section 2 explains the proposed model. Section 3 presents the derivation of the characteristic function of this model without considering regime-switching and then considering regime-switching. Further, the derivation of the pricing formula for the European options is presented. Section 4 consists of numerical illustrations. The final section makes concluding remarks with possible extensions of the model.
2. Model description
Consider a filtered probability space 
$(\Omega, \mathcal{F}, \mathcal{F}_{t \in [0,T]} , P)$ which models the uncertainty in the economy. In the absence of arbitrage, an equivalent martingale or risk-neutral measure 
$\mathcal{Q}$ exists. In this model, at time t, St represents the price of the underlying asset, νt is the volatility of the price, rt is the risk-free rate, Nt is the number of jumps up to and including time t, and Ji is the size of the 
$i\mathrm{th}$ jump. Let κ denotes the speed of mean reversion of volatility. The long-run mean of volatility is denoted by θt, and the volatility of the volatility is represented by ην. Assume that α represents the rate of mean reversion, σt symbolizes the volatility, and βt denotes the long-run mean of rt . The jump arrival process 
$\{N_t, t \in [ 0, T ]\}$ is a Poisson process with intensity process 
$\{\lambda_t, t \in [ 0, T ]\}$. Then, under the risk-neutral measure 
$\mathcal{Q}$, the asset price, the volatility, the risk-free rate, and the intensity processes are governed by the following stochastic differential equation:
\begin{align}
d S_{t}&=\left(r_{t}-\lambda_{t} m\right) S_{t} d t+\sqrt{\nu_{t}} S_{t} dB_{1,t}+\left(e^{J}-1\right) S_{t-} d N_{t},
\end{align}
\begin{align}
d \nu_{t}&=\kappa\left(\theta_{t}-\nu_{t}\right) d t+\eta_{\nu} \sqrt{\nu_{t}} dB_{2, t},
\end{align}
\begin{align}
d r_{t}&=\alpha \left(\beta_t-r_{t}\right) d t+\sigma_t d B_{3, t},
\end{align}
\begin{align}
d \lambda_{t}&= \kappa_{\lambda}\left(\theta_{\lambda}-\lambda_{t}\right) d t+\eta_{\lambda} \sqrt{\lambda_{t}} d B_{4, t},
\end{align} where κλ represents the rate of mean reversion, θλ is the long-run mean, and ηλ is the volatility of the intensity. Here 
$\{B_{1,t}, t \in [0,T]\}$, 
$\{B_{2, t}, t \in [0,T]\}$, 
$\{B_{3, t}, t \in [0,T]\}$, and 
$\{B_{4, t}, t \in [0,T]\}$ are standard Weiner processes. Furthermore, 
$dB_{1,t} dB_{2,t} = \rho dt$ and all other processes are independent. The probability density function of the jump size process, 
$J = \{J_i, i = 1, 2, \ldots, N_t\} $ is defined as:
\begin{equation}
g(x)= \frac{p e^{\frac{-x}{\eta_{1}} }}{\eta_{1}} 1_{\{x \geq 0\}}+ \frac{q e^{\frac{x}{\eta_{2}} }}{\eta_{2}} 1_{\{x \lt 0\}},
\end{equation} where 
$0 \lt \eta_{1} \lt 1, \eta_{2} \gt 0 ~ \text{and}~~ p+q=1.$ Let 
$m=\mathbb{E}_{\mathcal{Q}}[e^{J_1}-1]$ then
\begin{equation*}
m =\frac{p}{1-\eta_{1}}+\frac{q}{1+\eta_{2}}-1.
\end{equation*} Let θt, βt, and σt take values 
$\{\theta_1, \theta_2, \ldots, \theta_n\}$, 
$\{\beta_1, \beta_2, \ldots, \beta_n\}$, and 
$\{\sigma_1, \sigma_2, \ldots, \sigma_n\}$, respectively, according to a continuous time Markov chain 
$y = \{y_t: t \geq 0\}$ with n states whose unique stationary distribution exists. The states of y can be interpreted as different states of the economy or different phases of a business cycle. The business cycle alternates between phases of relatively fast growth in economic activity and phases of decrease in output. Without loss of generality, following [Reference Shen and Siu36], assume that the state space of y is identified with a set of unit vectors 
$\{e_1,e_2, \ldots,e_n\}$ where ei is a unit vector in 
$\mathbb{R}^n$ with jth component as Kronecker delta δij. Let 
$\langle \cdot,\cdot \rangle$ be the inner product on 
$\mathbb{R}^n$ where 
$\langle a,b \rangle = a^Tb$. For 
$\hat{\theta} = (\theta_1, \theta_2, \ldots, \theta_n)^T$, 
$\hat{\beta} = (\beta_1, \beta_2, \ldots, \beta_n)^T$, and 
$\hat{\sigma} = (\sigma_1, \sigma_2, \ldots, \sigma_n)^T$, then
\begin{equation*}
\theta_t = \langle \hat{\theta}, y_t\rangle,
~~\beta_t = \langle \hat{\beta}, y_t\rangle,
\,\text{and} ~~ \sigma_t = \langle \hat{\sigma}, y_t\rangle.\end{equation*} Let 
$\mathbb{Q} = [q_{ij}]$ be the transition rate matrix of Markov chain y where qij is the transitional intensity of the chain from regime i to regime j or from state ei to ej. Then,
\begin{equation}
dy = \mathbb{Q}^T y dt +dM,
\end{equation} where M is a vector martingale process. Let yt be independent of 
$B_{1,t}, B_{2, t}, B_{3, t}, B_{4, t}$, and Nt. Assume the initial values 
$S_0 \gt 0$, 
$\nu_0 \gt 0$, 
$r_0 \gt 0$, and 
$\lambda_0 \gt 0$ are given.
Under the risk-neutral measure, the price of the European options is given by:
\begin{equation*}
C(0):=\mathbb{E}_{\mathcal{Q}}\left[e^{-\int_{0}^{T} r(s) d s}\left(S_T-K\right)^{+} {\bigg|} \mathcal{F}_0\right].
\end{equation*} To solve this expression, a new measure is introduced. Define the forward measure 
$\tilde{\mathcal{Q}}$ by
\begin{equation}
\frac{d \tilde{\mathcal{Q}}}{d \mathcal{Q}}=\frac{e^{-\int_{0}^{T} r_u d u}}{\mathbb{E}_
\mathcal{Q}\left[e^{-\int_{0}^{T} r_{u} d u} \right] }.
\end{equation} Let 
$ B(t, T)$ be the price of a bond at time t maturing at time T. Then, 
$ B(t, T)= \mathbb{E}_\mathcal{Q} \left[e^{-\int_{t}^{T} r_{u} d u} {\bigg|} \mathcal{F}_t\right]$. According to Section 3 of [Reference Shen and Siu36], the value 
$B(0, T)$ is given by
\begin{equation}
B(0, T)= R(T,y_0)e^{-P(T) r_0},
\end{equation} where P(T) and 
$R( T, y_0)$ are given by
\begin{align}
P(T) &= \frac{1}{\alpha}[1-e^{-\alpha T}],
\end{align}
\begin{align}
R( T,y_0) &= \mathbb{E}_\mathcal{Q}\left[ \exp{\left\{\int_{0}^{T}\left(\alpha \beta_s P(s)-\frac{1}{2} \sigma_s^2 P^{2}(s)\right) d s\right\}} \bigg| y_0 \right].
\end{align} Let 
$\tilde{R}(t)= \operatorname{exp}\left\{\int_{0}^{t} \langle\alpha \hat\beta P(u)-\frac{1}{2} \hat\sigma^2 P^{2}P(u), y_u\rangle du \right\} $ and 
$\tilde{\mathbf{R}}(t) = \tilde{R}(t) y_t. $ Then,
\begin{align*}
d \tilde{\mathbf{R}}(t) &= \tilde{R}(t) dy_t + d\tilde{R}(t) y_t,\\
&= \tilde{R}(t)\left(\mathbb{Q}^Ty_tdt +d M_t\right) + \left(\langle \alpha \hat\beta P(t)-\frac{1}{2} \hat\sigma^2 P^{2}(t) , y_t \rangle\right) \tilde{R}(t) y_t dt,\\
&= \left(\mathbb{Q}^T + \langle \alpha \hat\beta P(t)-\frac{1}{2} \hat\sigma^2 P^{2}(t) , y_t \rangle\right) \tilde{R}(t) y_tdt + \tilde{R}(t)d M_t, \\
&= \left(\mathbb{Q}^T + \operatorname{diag}\left[\alpha \hat\beta P(t)-\frac{1}{2} \hat\sigma^2 P^{2}(t)\right]\right) \tilde{R}(t) y_t dt + \tilde{R}(t)d M_t.
\end{align*} Let 
$\hat{\mathbf{R}}(t) = \mathbb{E}_{\mathcal{Q}}[\tilde{\mathbf{R}}(t)|y_0]$ then 
$\hat{\mathbf{R}}(t)$ will satisfy,
\begin{equation}
\frac{d\hat{\mathbf{R}}(t)}{d t} =\left(\mathbb{Q}^{T}+\ \operatorname{diag}\left[\alpha \hat\beta P(s)-\frac{1}{2} \hat \sigma^2 P^{2}(s) \right] \right)\hat{\mathbf{R}}(t).
\end{equation}Following the similar method mentioned in [Reference Elliott and Nishide4] gives
\begin{equation*}
\mathbb{E}_{\mathcal{Q}}[\tilde{{R}}(t)|y_0] =
\langle\hat{\mathbf{R}}(t), 1_n \rangle
= \langle \operatorname{e}^U y_0 , 1_n \rangle,
\end{equation*} where 
$1_n = \{1,1,\ldots,1\}^T$ and
\begin{equation}
U= \int_0 ^ t \left[\mathbb{Q}^{T}+\ \operatorname{diag}\left(\alpha \hat\beta P(s)-\frac{1}{2} \hat \sigma^2 P^{2}(s) \right) \right] ds.
\end{equation} The price of the call options under the forward measure 
$\tilde{\mathcal{Q}}$ will be
\begin{equation}
C(0) = \ B(0, T) \left(\mathbb{E}_{\tilde{\mathcal{Q}}}\left[ \left(S_T-K\right)^{+}{\bigg|}\mathcal{F}_0\right]\right) .
\end{equation} The dynamics of 
$X_t,$ νt, and λt remain same under the forward measure 
$\tilde{\mathcal{Q}}$ as in the risk neutral measure 
$\mathcal{Q}$ whereas the dynamics of rt under 
$\tilde{\mathcal{Q}}$ will be given by
\begin{align*}
d r_t&=\left(\alpha(\beta_t - r_t) -\sigma_t^{2} P(t, T)\right) d t+\sigma_t d \tilde{B}_{3, t}, \\
d \tilde{B}_{3, t}&= d B_{3, t}+ \sigma_t P(t, T) d t , \end{align*} where 
$d\tilde{B}_{3, t}$ is Brownian motion under measure 
$\tilde{\mathcal{Q}}$. Additionally, assume that the rate matrix of y under measure 
$\tilde{\mathcal{Q}}$ is given by 
$\tilde{\mathbb{Q}} = [\tilde{q_{ij}}] $ then
\begin{equation*} \tilde{q_{ij}} =
\begin{cases}
q_{ij} \frac{\operatorname{exp}(R(0,T,e_j))}{\operatorname{exp}(R(0,T,e_i))}& i \neq j \\
-\sum_{k \neq i} q_{ki}\frac{\operatorname{exp}(R(0,T,e_k))}{\operatorname{exp}(R(0,T,e_i))} & i=j
\end{cases}.
\end{equation*}3. Options pricing
This section develops a pricing formula for the European options within the framework outlined in Section 2. Proposition 3.1 provides the characteristic function of the log process of asset price without considering the regime-switching of the volatility process. Let 
$X_t = \operatorname{log}S_t$, then
\begin{equation}
d X_{t}=\left(r_{t}-\lambda_{t} m-\frac{1}{2} \nu_{t}\right) d t+\sqrt{\nu_{t}} dB_{1,t}+d \sum_{k=1}^{N_{t}} J_{k}.
\end{equation}Proposition 3.1. The conditional characteristic function 
$\varphi(\delta, \tau, x, \nu, r, \lambda )$ of XT under the forward measure 
$\tilde{\mathcal{Q}}$ is given by
\begin{equation}
\varphi\left(\delta,\tau, x, \nu, r, \lambda \right) = \operatorname{exp}\left\{{A(\delta, \tau)+ B(\delta, \tau) \nu+ C(\delta, \tau)\lambda +D(\delta, \tau) r+i \delta x}\right\},
\end{equation}where
\begin{align*}
D (\delta, \tau)&= \frac{i \delta}{\alpha}(1 - e^{-\alpha \tau}),\\
B(\delta, \tau) &= \frac{2}{\eta_{\nu}^{2}}\left[b(\delta) + d(\delta) \left(\frac{d(\delta) \operatorname{sinh}(d(\delta) \tau) - b(\delta) \operatorname{cosh}(d(\delta)\tau)}{d(\delta) \operatorname{cosh}(d(\delta) \tau) - b(\delta) \operatorname{sinh}(d(\delta)\tau)}\right)\right],\\
C(\delta, \tau)&=\frac{2}{\eta_{\lambda}^{2}} \left[ -\kappa_{\lambda} + d'(\delta) \left(\frac{d'(\delta) \operatorname{sinh}(d'(\delta) \tau) - \kappa_{\lambda} \operatorname{cosh}(d'(\delta)\tau)}{d'(\delta) \operatorname{cosh}(d'(\delta) \tau) - \kappa_{\lambda} \operatorname{sinh}(d'(\delta)\tau)}\right)\right],\\
A (\delta, \tau) &= \int _0 ^ \tau \beta_{u} D(\delta, u) du +\kappa\int _0 ^ \tau \theta_u B(\delta,u) du - \int _0 ^ \tau \sigma_{u}^2P(\delta,u)D(\delta,u) du \\
&\quad + \frac{1}{2}\int _0 ^ \tau \sigma_{u}^2 D^2(\delta,u) du
+ \frac{2 \kappa_{\lambda} \theta_{\lambda}}{\eta_{\lambda}^{2}} \left[ -\kappa_{\lambda} \tau + \operatorname{ln}\left(d'(\delta) \operatorname{cosh}(d'(\delta) \tau) - \kappa_{\lambda} \operatorname{sinh}(d'(\delta)\tau)\right)\right],
\end{align*}and
\begin{equation*}
b(\delta) = \frac{(i \rho \eta_{\nu} \delta-\kappa)}{2},
d(\delta) = \frac{\sqrt{b(\delta)^{2} + \eta_{\nu}^{2} \delta(i + \delta)}}{2}, d'(\delta)= \frac{\sqrt{\kappa_{\lambda}^{2} + 2 \eta_{\lambda}^{2} \left(m i \delta - \Lambda(\delta)\right)}}{2}.
\end{equation*}Proof. Let 
$\tau =T-t$. The conditional characteristic function of XT by definition is 
$\varphi(\delta, \tau, x, \nu, r, \lambda ) = \mathbb{E}_{\tilde{\mathcal{Q}}}[e^{i \delta X_T} {\big|} \nu_t =\nu, r_t=r, X_t=x, \lambda_t = \lambda, y_T ]$. Using Feynman–Kac theorem, 
$\varphi(\delta, \tau, x, \nu, r, \lambda )$ will satisfy the following partial differential equation (PDE),
\begin{align*}
\frac{\partial \varphi}{\partial \tau}&=\left(r -\frac{1}{2} \nu -\lambda m\right) \frac{\partial \varphi}{\partial x}+\frac{1}{2}\nu_t\frac{\partial^{2} \varphi}{\partial x^{2}}+\left(\alpha
(\beta - r)-\sigma^{2} P(t,T)\right) \frac{\partial \varphi}{\partial r} \\
&\quad +\frac{1}{2} \sigma^{2} \frac{\partial^{2} \varphi}{\partial r^{2}}+\kappa \left(\theta_t - \nu\right) \frac{\partial \varphi}{\partial \nu}+\frac{1}{2} \eta_{\nu}^{2} \nu \frac{\partial^{2} \varphi}{\partial \nu^{2}}+
\rho \eta_{\nu} \nu \frac{\partial^{2} \varphi}{\partial X \partial \nu}+\kappa_{\lambda}\left(\theta_{\lambda}- \lambda\right) \frac{\partial \varphi}{\partial \lambda}\\ &\quad +\frac{1}{2} \eta_{\lambda}^{2} \lambda\frac{\partial^{2} \varphi}{\partial \lambda^{2}}+\lambda \int_{- \infty}^{+\infty}\left[\varphi\left(\delta, \tau, x+j, \nu, r, \lambda\right)-\varphi\left(\delta, \tau, x, \nu, r, \lambda \right)\right] g(j) d j.
\end{align*} Let 
$\Lambda(\delta) = \frac{p}{1-i\delta \eta_1} + \frac{q}{1+i \delta \eta_2}-1 $, then
\begin{equation}
\lambda \int_{-\infty}^{\infty}\left[\varphi\left(\delta, \tau, x+j, \nu , r, \lambda \right)-\varphi\left(\delta, \tau, x, \nu, r, \lambda\right)\right] g(j) dj =\Lambda(\delta) \varphi\left(\delta; \tau, x, \nu, r,\lambda\right).
\end{equation} Assume that the characteristic function 
$\varphi(\delta, \tau, x, \nu, r, \lambda ) $ has the following affine form:
\begin{equation}
\varphi\left(\delta, \tau, x, \nu, r, \lambda \right)=e^ {A(\delta, \tau)+ B(\delta, \tau) \nu+ C(\delta, \tau)\lambda +D(\delta, \tau) r+i \delta x} .
\end{equation} Substituting 
$\varphi(\cdot)$ in the above PDE, the following set of ordinary differential equations is obtained with the boundary conditions 
$A(\delta,0) = B(\delta,0) = C(\delta,0) = D(\delta,0) = 0$
\begin{equation*}
\begin{aligned}
\frac{d B}{d \tau} &=\frac{1}{2} \eta_{\nu}^{2} B^{2}+\left(i \rho \eta_{\nu} \delta-\kappa \right) B-\frac{i}{2} \delta-\frac{1}{2}\delta^{2}, \\
\frac{d C}{d \tau}&= \frac{1}{2} \eta_{\lambda}^{2} C^{2} - \kappa_{\lambda} C - m i \delta +\Lambda(\delta), \\
\frac{d D}{d \tau}&= -\alpha D +i \delta,\\
\frac{d A}{d \tau} &=\beta \alpha D -\sigma^{2} P D +\frac{1}{2} \sigma^{2} D^{2}+\kappa \theta B+\kappa_{\lambda} \theta_{\lambda} C .
\end{aligned}
\end{equation*} Solving the above ordinary differential equations (ODEs) will give the required values of 
$A(\delta, \tau),$ 
$ B (\delta, \tau),$ 
$ C(\delta, \tau)$, and 
$D(\delta, \tau)$.
Now, to calculate the characteristic function of XT with the regime-switching of the long-run mean of the volatility, mean reversion level of risk-free rate and volatility of risk-free rate, a method similar to solve Eq. (10) is employed.
Let 
$\psi(\delta,u) = (\psi_1(\delta,u), \psi_2(\delta,u), \ldots, \psi_n(\delta,u))^T$ with
\begin{equation*}
\psi_k(\delta,u) = \left(\kappa \theta_k B(\delta, u)+\beta_k \alpha D(\delta, u) -\sigma_k^{2} P(u,T) D(\delta, u) +\frac{1}{2} \sigma_k^{2} D^{2}(\delta, u) \right).
\end{equation*}The characteristic function of XT considering regime-switching is
\begin{align*}
\mathbb{E}_{\tilde{\mathcal{Q}}}\left[e^{i \delta X_T} {\big|} \mathcal{F}_t\right]&= \mathbb{E}_{\tilde{\mathcal{Q}}}\left[\mathbb{E}_{\tilde{\mathcal{Q}}}\left[e^{i \delta X_T} {\big|} \nu_t =\nu, r_t=r, X_t=x, \lambda_t = \lambda, y_T \right]{\bigg|}y_t\right]\\
&=\operatorname{exp}\{ A_1 (\delta, \tau) + B(\delta, \tau) \nu + C(\delta, \tau)\lambda +D(\delta, \tau) r +i \delta x \}\\
& \quad \times \mathbb{E}_{\tilde{\mathcal{Q}}}\left[\operatorname{exp}\{\int _0 ^ \tau \langle\psi(\delta,u),y_u\rangle du\} \bigg| y_t\right],
\end{align*} where 
$A(\delta, \tau) = A_1 (\delta, \tau) + \int _0 ^ \tau \langle\psi(\delta,u),y_u\rangle du$. Then,
\begin{equation*}
\mathbb{E}_{\tilde{\mathcal{Q}}} \left[\operatorname{exp} \left\{\int _0 ^ \tau \langle\psi(\delta, u),y_u\rangle du\right\}\bigg| y_t\right] = \langle \Psi(\delta, 0, \tau)y_t , 1_n \rangle,
\end{equation*} where 
$1_n = \{1,1,\ldots,1\}^T$ and Ψ is a n × n matrix that satisfies the following differential equation
\begin{equation*}
\frac{d \Psi(\delta, u, \tau)}{d \tau} = \left(\mathbb{\tilde{Q}}^{T}+ \operatorname{diag}(\psi(\delta,\tau)) \right) \Psi(\delta, u, \tau),
\end{equation*} with 
$\Psi(\delta, u, u)=I$ for 
$\tau \geq u^3$. Hence, the characteristic function of XT considering regime-switching is
\begin{equation}\varphi(\delta,\tau)=\operatorname{exp}\{A_1(\delta,\tau)+B(\delta,\tau)\nu+C(\delta,\tau)\lambda+D(\delta,\tau)r+i\delta x\}\langle\Psi(\delta,t,\tau)y_t,1_n\rangle.\end{equation}The price of the call options will be
\begin{equation*}
\begin{aligned}
C(0) &= \ B(0, T) \left(\mathbb{E}_{\tilde{\mathcal{Q}}}\left[ \left(S_T-K\right)^{+}{\bigg|}\mathcal{F}_0\right]\right)\\
&= B(0, T) \left(\mathbb{E}_{\tilde{\mathcal{Q}}}\left[S_T 1_{\{X_T \geq k\}}\right]-e^{k} \mathbb{E}_{\tilde{\mathcal{Q}}}\left[1_{\{X_T \geq k\}}\right]\right)\\
&= B(0, T) \left(\mathbb{E}_{\tilde{\mathcal{Q}}}\left[S_T 1_{\{X_T \geq k\}}\right] -e^k \tilde{P}\left( X_T \geq k\right)\right).
\end{aligned}
\end{equation*} Let 
$f(\cdot)$ be probability density function of XT under measure 
$\tilde{\mathcal{Q}}$, then
\begin{equation*}
\begin{aligned}
\tilde{P}(X_T \geq k)&= \int_k ^ \infty f(y) dy ~~= \int_k ^ \infty \left(\int_{- \infty}^\infty \frac{1}{2\pi} e^{-i \delta y}\varphi(\delta, \tau) d \delta \right) dy\\
&= \frac{1}{2}+\frac{1}{2\pi} \int_{-\infty}^{\infty} \frac{e^{-i \delta k} \varphi(\delta, \tau)}{i \delta} d \delta \\
&= \frac{1}{2}+\frac{1}{\pi} \int_0^{\infty} \operatorname{Re}\left(\frac{e^{-i \delta k} \varphi(\delta, \tau)}{i \delta}\right) d \delta,
\end{aligned}
\end{equation*} where 
$i = \sqrt{-1}$. Consider another function 
$f_1(y)$ defined by 
$f_1 (y) = \frac{e^y f(y)}{\varphi(-i, \tau)}$, where 
$\varphi(-i, \tau)= \int_{-\infty} ^ \infty e^y f(y) dy$. It can be verified that it is a probability density function for some random variable. The characteristic function corresponding to this random variable will be 
$\varphi_1(u,\tau) = \frac{\varphi(u-i, \tau)}{\varphi(-i, \tau)}$. Now consider,
\begin{equation*}
\begin{aligned}
\mathbb{E}_{\tilde{\mathcal{Q}}}\left[S_T 1_{\{X_T \geq k\}}\right]
& = \int_k ^ \infty e^y f(y) dy ~~ = \varphi(-i, \tau) \left(\int_k ^ \infty f_1 (y) dy
\right)\\
& = \varphi(-i, \tau) \left(\frac{1}{2}+\frac{1}{\pi} \int_0^{\infty} \operatorname{Re}\left(\frac{e^{-i \delta k} \varphi_1(\delta, \tau)}{i \delta }\right) d \delta\right)\\
& = \varphi(-i, \tau)\left(\frac{1}{2}+\frac{1}{\pi} \int_0^{\infty} \operatorname{Re}\left(\frac{e^{-i \delta k} \varphi(\delta-i, \tau)}{i \delta \varphi(-i, \tau)}\right) d \delta\right)\\
& = \frac{\varphi(-i, \tau)}{2}+ \left(\frac{1}{\pi} \int_0^{\infty} \operatorname{Re}\left(\frac{e^{-i \delta k} \varphi(\delta-i, \tau)}{i \delta}\right) d \delta\right).
\end{aligned}
\end{equation*}4. Numerical and sensitivity analysis
This section compares the price of options obtained using the derived formula with those generated through MCSs. Additionally, to highlight the importance of incorporating regime-switching and stochastic intensity in options pricing, a comparison is presented. The sensitivity of various parameters to options price is also exhibited. Unless otherwise specified, the values of the selected parameters are provided in Table 1. It is assumed that there are two regimes: one representing the state of economic boom (Regime 1) and the other economic crisis (Regime 2). In times of economic growth, the financial markets are driven by high interest rates and low volatility. Conversely, economic recessions are marked by low interest rates and high volatility. This section is divided into three subsections: Subsection 4.1 presents the numerical calculations, Subsection 4.2 provides the comparisons with other models, and Subsection 4.3 examines the sensitivity of option prices to various parameters across different strike prices.
The value of parameters.
4.1. Numerical calculation
To validate the accuracy of the derived formula, a comparison with MCSs is presented in Table 2 for option prices at various strike prices K, assuming that the economy is initially in regime 2. To calculate the options price from the MCS, 
$100,000$ sample paths are generated, each consisting of 500 time steps. It can be observed that the prices obtained from the derived formula (MP) are point-wise close to the corresponding prices from the MCS. The relative error between the two prices is less than 
$0.4\%$ which verifies the accuracy of the formula. Figure 1 illustrates the option price when the economy is in regime 1 (economic growth) and regime 2 (economic recession). It is observed that the price of the call options increases as the strike price decreases. Additionally, Figure 1 shows that option prices are higher during an economic recession than during economic growth. This can be explained by the fact that option prices tend to increase when the underlying asset’s volatility rises and the asset’s returns are more volatile during a recession than economic growth.
Comparison of option prices against strike prices in regime 1 and regime 2.
European call options price calculated through proposed model and Monte Carlo simulation.
Abbreviations: MCS, Monte Carlo simulation; MP, model price; RE, relative error.
4.2. Model comparison
In this section, a comparative analysis of the proposed models with several models is conducted. The primary objective of this comparison is to demonstrate the significance of the proposed model.
4.2.1. Comparison with the Heston-CIR model with stochastic intensity
Figure 2 illustrates the comparison between the options price derived from the proposed model and those derived from the model proposed by [Reference Ma, Chen and Lyu31]. It is assumed that initially the economy is in regime 1. The graphical comparison demonstrates the differences and similarities in the options price for a range of strike prices K. Both models exhibit a decrease in options price as the strike price increases; however, the option prices in the model proposed in this article are higher.
Comparison of option prices between the proposed model and the CIR model.
4.2.2. Comparison with hybrid Heston–Hull–White model with stochastic intensity
The comparison of the European call options price from the proposed model against those derived from the hybrid Heston–Hull–White model with stochastic intensity and without regime-switching is illustrated in Figure 3. If the economy is in regime 1, the markets are driven by high interest rates and low volatility, which results in lower risk, and hence option prices are less. In regime 2, the option prices are higher as compared to regime 1 as in regime 2 market has lower interest rates and higher volatility. Figure 3 illustrates that for regime 1, the option prices calculated by the regime-switching model are higher compared to the model without regime-switching, whereas for regime 2, the option prices are lower for the regime-switching model than without regime-switching model. This is because, in the regime-switching model, there is a positive probability that the economy will shift into another regime leading to difference in prices.
Comparison of option prices against strike prices for the proposed model and the hybrid Heston–Hull–White model with stochastic intensity while considering the initial regime as Regime 1 (on the left) and Regime 2 (on the right).
4.2.3. Comparison with regime-switching hybrid Heston–Hull–White model without stochastic intensity
To illustrate the impact of stochastic intensity, the comparison of the proposed model with the regime-switching hybrid Heston–Hull–White model having constant intensity is depicted in Figure 4. It is observed that, regardless of the initial state of the economy, the options price is higher in stochastic intensity than in constant intensity.
Comparison of option prices against strike prices for the proposed model and the regime-switching hybrid Heston-Hull-White model without stochastic intensity while considering the initial regime as Regime 1 (on the left) and Regime 2 (on the right).
4.3. Sensitivity analysis
This subsection depicts the impact of the sensitivity of various parameters of the risk-free rate, volatility, and intensity processes on option prices for different strike prices. For all of these figures, it is assumed that initially, the market is in regime 2.
Initial volatility (ν 0): Figure 5 exhibits a plot of the option prices and strike prices for various values of initial volatility. It is observed that as volatility increases, the options price at all strike prices tends to rise. However, at lower strike prices, the difference in options price becomes less significant. One possible explanation is that the option price depends on the intrinsic values of the option. The option’s intrinsic value is higher at a lower strike price. Hence, the effect of volatility is not very significant.
Option prices against strike prices for different initial volatility levels.
Initial risk-free rate (r 0): Figure 6 presents how varying the values of r 0 affect the option prices by plotting option prices for different r 0 values against various strike prices K. The risk-free rate affects option price in two different ways. First, it impacts the underlying asset’s price. Second, it is used as a discounting factor to determine the present value of the options price. The options price increases with the increase in value of r 0. The impact becomes less significant as the strike prices increase.
Option prices against strike prices for different initial risk-free rates.
Mean reversion level of risk-free rate (βt): To examine the impact of mean reversion level, it is assumed that 
$\beta_1 = \beta_2 + 0.02$. Option prices at various strike prices K corresponding to different values of the βt are depicted in Figure 7. As the value of βt increases, the options price also increases across all strike prices. This is because the higher value of βt, which is the long-run mean of the risk-free rate, implies an expectation of increasing risk-free rates.
Option prices against strike prices for different values of βt.
Mean reversion speed of risk-free rate (α): The effect of varying α on option prices is investigated in Figure 8. This figure analyzes the option prices at various strike prices corresponding to different values of mean reversion speed (α). The figure shows that the higher value of α leads to an increased options price for every strike price.
Option prices against strike prices for different values of α.
Volatility of risk-free rate (σt): Figure 9 displays the options price against the parameter σt at various strike prices K. To plot this figure, it is assumed that 
$\sigma_2 = $ 
$ \sigma_1 + 0.02$. It is observed that option prices are positively correlated with the volatility of risk-free rate. A higher σt increases option prices at every strike price.
Option prices against strike prices for different values of σt.
Long-run mean of intensity process (θλ): Figure 10 illustrates the options price at different values of long-run mean of intensity process. From the figure, it can be observed that as the value of θλ increases, there is a corresponding increase in options price for all strike prices. This is because as the value of θλ increases, the expectation of jump intensity increases.
Option prices against strike price for different values of θλ.
5. Conclusions and future work
This study proposes a pricing formula for European options using the regime-switching framework with stochastic intensity. The closed-form formula for options price is derived using the Feynman–Kac theorem with asset volatility and risk-free rate controlled by the economic environment. The results highlight the importance of incorporating stochastic intensity in option pricing and its significant impact. Furthermore, a comparison is made with models without regime-switching framework to depict the difference in option prices. This divergence occurs because the regime-switching model adjusts for fluctuating market conditions and volatility patterns over time, resulting in more accurate pricing that reflects economic movements. As a result, including regime-switching can result in more accurate pricing, better reflecting the complexity of real-world financial markets.
An interesting extension may involve studying the valuation of American options within this framework, considering their flexibility to be exercised at any time before or on the maturity date. Considering the relationship between the intensity process and the jump arrival process allows for a more realistic scenario. Hence, a more general jump arrival process like the Hawkes process may be considered for asset price modeling.
Conflict of interest
The authors declare that they have no competing interests.
Open access
