An analytical approximation of European option prices under a hybrid GARCH-Vasicek model with double exponential jump in the bid-ask price economy

1.   Introduction

Most of the existing option pricing models assume that the market is frictionless and completely liquid, which leads to the situation that the bid and ask prices are treated as the same when pricing various options. However, there are two kinds of prices listed in real markets, i.e., bid and ask prices. The former refers to the price at which traders sell options, while the latter represents the one at which traders buy options. This prompted the development of conic finance theory ^[1,2,3], which acknowledges the non-uniqueness of market prices. It basically assumes that the market acts as a central counterparty; it sells and buys assets at ask price and bid price, respectively. The difference between ask and bid prices is usually called the bid-ask spread, which is an indication of the market liquidity.

In fact, the issue of liquidity has already become a hot topic in the area of risk management and has attracted a lot of research attention on the bid-ask economy. Madan and Cherny ^[1] solved bid and ask prices of European options analytically under a market liquidity model after introducing a single market stress level. However, a number of authors used market data to further study the model and found that the implied market liquidity is not constant, which is far from the assumption in Madan and Cherny's conic option pricing model ^[4,5]. It should be pointed out that the underlying asset price in the above literature is modeled by the simplest geometric Brownian motion (GBM) in order to obtain the explicit forms of the distribution function so that bid and ask prices can be further derived by numerical methods.

It is well known that the GBM is not suitable for describing the price process of the underlying asset since the so-called implied volatility smile or skew has been widely observed. Many studies have extended the GBM model. Mehrdoust and Najafi ^[6] studied European option pricing under a fractional Black-Scholes model with a weak payoff function. Lin and He ^[7] proposed a regime switching fractional Black-Scholes model and obtained the European option pricing formula. Hassanzadeh and Mehrdoust ^[8] investigated option pricing under a multifactor uncertain volatility model. Nevertheless, one of the most famous modifications is the Heston stochastic volatility model ^[9], and there are also various modifications to the Heston model ^[10]. It should be noted that the Heston model with the square root specification cannot describe the nonlinear characteristics of financial time series well in practical applications, although it brings great convenience due to the analytical pricing formula of options. Thus, alternative models have been established, among which the GARCH diffusion model has received increasing attention. This is because the GARCH diffusion model is shown to be able to better describe financial time series ^[11,12]. We refer interested readers to Kaeck and Alexander ^[13] and the references therein for more details on the results associated with the GARCH diffusion model.

On the other hand, the jump-diffusion model represents a refinement of the GBM process, as it effectively captures the discontinuous changes in the underlying stock returns ^[14]. Moreover, the double exponential jump-diffusion model proposed by Kou ^[15] is able to reflect high levels of skewness and leptokurtosis exhibited by financial data. Due to its analytical tractability, this model has received great attention from academia and industry ever since it was put forward. For example, Mehrdoust et al. ^[16] presented the valuation for European options by adding jumps into the Bi-Heston model. Huang and Guo ^[17] got the semi-analytic solution of vulnerable options by assuming that the price process of the underlying asset follows non-affine stochastic volatility with double exponential jump. Considering the stochastic behaviour and jump risks, Hu et al. ^[18] investigated the pricing of European crude oil options. Moreover, we should be aware that the spot interest rate plays a decisive role in modern financial industry and it changes stochastically in the market. Grzelak et al. ^[19], Recchioni and Sun ^[20], and Chen et al. ^[21] proved that it performs much better if the option pricing model replaces the constant interest rate with a stochastic one. Although He and Zhu ^[22] presented a closed-form series solution to European option prices when the volatility and interest rate are both stochastic, which is appealing, their assumption of the CIR interest rate model prevents the interest rate from going negative. This is inconsistent with actual situations, as the short-term government bond markets of the USA and Europe have already witnessed negative interest rates ^[23]. It was even claimed by Recchioni et al. ^[24] that models allowing interest rates to take negative values are able to improve the performance of option pricing and implied volatility forecasting.

Considering all the features discussed above that are able to help improve model performance, we incorporate the Vasicek stochastic interest rate and double exponential jump into the GARCH diffusion model when pricing options in the bid-ask price economy. However, the probability distribution of the log-price process cannot be analytically obtained due to the complexity of the adopted dynamic processes. Fortunately, with the help of the Fourier-cosine method, or the COS method^{* [25,26]}, we are able to derive the probability density function by making use of the corresponding characteristic function. Once the density function is obtained, bid and ask prices can be straightforwardly computed with some numerical schemes, including Gaussian quadrature. The accuracy of our proposed approach is verified via numerical comparison with the Monte Carlo simulation, and the sensitivity analysis is also performed so that the effect of the market liquidity parameter on bid and ask prices is clear.

^*The COS method has been widely used in European and American options pricing due to its simplicity, high pricing accuracy, and high computational efficiency. The central idea of the COS method is to reconstruct the density function by the characteristic function of log-price.

The remainder of the paper is as follows. Section 2 presents a hybrid option pricing model combining the Vasicek stochastic interest rate, GARCH diffusion volatility model, and double exponential jump together. An approximation to the characteristic function is derived in Section 3. In Section 4, using the Fourier cosine method, we obtain the density function of the underlying log-price, and, further, we derive bid-ask prices of European options. Results of numerical experiments are provided in Section 5, with the last section concluding the article.

2.   Model specification

Consider a filtered probability space $\left\{\Omega, \mathcal{F}, \left\{\mathcal{F}_{t}\right\}_{0 \leq t \leq T}, \mathbb{Q}\right\}$, with $\mathbb{Q}$ being a risk-neutral probability measure. We assume that the market is frictionless, i.e., there are no transaction costs and the assets can be traded continuously. The underlying asset price process $S_{t}$, the volatility process $v_{t}$, and the stochastic interest rate process $r_{t}$ under $\mathbb{Q}$ are

with $\mathrm{d}W_{s}(t)\mathrm{d}W_{v}(t) = \rho\mathrm{d} t$. It should be pointed out that the model can be more sophisticated if one incorporates the correlation between the underlying price and interest rate. However, introducing such correlation could break down the analytical tractability according to a number of different literature ^{[19,27,28,29,30]}. Thus, it remains an open question on how to effectively price options with the incorporation of such correlation, and we would like to leave this to future work. The mean reversion speed is denoted by $\kappa_{v}$, while $\theta_v$ and $\sigma_{v}$ respectively represent the long-run mean and instantaneous volatility of volatility. $\lambda$ is the constant intensity of the Poisson process $N_{t}$. We have that $m = E^{\mathbb Q}(e^J-1)$, where $J$ is the jump size following an asymmetric double exponential distribution whose density function can be presented as

where $p$ and $q$ respectively represent the probability of jumping upward and that of jumping downward, with $p+q = 1$. This further implies that $m = \frac{p \eta_{1}}{\eta_{1}-1}+\frac{q \eta_{2}}{\eta_{2}+1}-1$.

Due to the existence of the stochastic interest rate, it is a natural treatment to make a measure transform to convert the price dynamics of the underlying asset under the original measure $\mathbb{Q}$ into those under the T-forward measure $\mathbb{Q}^T$. In order to achieve this, we need to choose the T-discount bond price as the numeraire, which will be provided below. In particular, if we denote $P(t, T)$ as the price of a risk-free zero-coupon bond maturing at time $T$, when the evolution of $r_{t}$ follows Eq (2.3), the $P(t, T)$ can be formulated as

where

and

The measure changing from the risk-neutral probability measure $\mathbb{Q}$ to the T-forward measure $\mathbb{Q}^T$ can be established by the following Radon-Nikodym derivative:

where

Thus, we can express the target model dynamics under $\mathbb{Q}^T$ as

where $\mathrm{d}W_{v}(t)\mathrm{d}W_{v}^{\bot}(t) = 0$.

Letting $x_{t} = lnS_{t}$, we can transform Eqs (2.4)–(2.6) into the following form:

3.   The pricing of European options in bid-ask price economy

3.1.   The COS method

In this subsection, the COS method is briefly introduced for the completeness of the paper. It is well-known that the price $P(x, t)$ of a European option at time $t$ is an expectation under the risk neutral measure according to the classic option pricing theory ^[31], which does not take bid-ask spread into account, i.e.,

where $x = lnS_{t}$, $y = lnS_{T}$, $r$ is the risk-free rate, and $T$ is the maturity time. $f(y|x)$ is the probability density function of the underlying process, and $P(y, T)$ is the payoff function of the option at maturity.

Without significantly losing accuracy, given the special choice of $[a, b]$, we can obtain

The key point of the COS method is that the density $f(y|x)$, which is unknown in most cases, is approximated by a Fourier-cosine series expansion on $[a, b]$, i.e.,

where $A_k(x) = \frac{2}{b-a} \int_a^b f(y \mid x) \cos \left(k \pi \frac{y-a}{b-a}\right) d y$, for $k = 0, 1, \ldots, N-1.$

As a result,

where $V_k = \frac{2}{b-a} \int_a^b P(y, T) \cos \left(k \pi \frac{y-a}{b-a}\right) d y$, and the prime of $\sum^{\prime}$ is used to indicate that the first term of the summation should be multiplied by a weight of $1/2$.

Meanwhile, $A_k(x)\approx F_k(x)$. Since

we have

where $\operatorname{Re}\{\cdot\}$ is an operator to take the real part.

In summary,

where $\Phi$ is the characteristic function. Detailed contents about the COS method can be found in Fang and Oosterlee ^[25].

3.2.   The joint characteristic function

This subsection presents the joint characteristic function of the underlying log-price, stochastic volatility, and stochastic interest rate. The following theorem provides the solution to the joint characteristic function under the T-forward measure $\mathbb{Q}^{T}$.

Theorem 3.1. Given that the underlying asset price follows the dynamics in Eqs (2.7)–(2.9), the characteristic function for $x_{T}$ is given by

where

with

and $T \geq t, \tau = T-t, i = \sqrt{-1}$, ${\Lambda(u) = \frac{p \eta_{1}}{\eta_{1}-i u}+\frac{q \eta_{2}}{\eta_{2}+i u}-1}$.

Proof. By applying the Feynman-Kac theorem, $\Phi(x, v, r, \tau; u)$ satisfies the following partial integral-differential equation (PIDE):

The boundary condition for Eq $(3.1)$ is given by

According to several different studies ^[32,33,34], this PIDE has an exponential-affine solution of the form

with boundary conditions

Moreover,

where $\Lambda(u) = \frac{p \eta_{1}}{\eta_{1}-i u}+\frac{q \eta_{2}}{\eta_{2}+i u}-1$.

Equation $(3.1)$ is very difficult to solve since it is a nonlinear PDE, and thus we first linearize it approximately. The idea is to approximate $v^{\frac{3}{2}}, v^{2}$ in the PIDE using Taylor expansions around the long-run mean of variance as follows:

With the substitution of Eqs (3.2)–(3.4) into PDE $(3.1)$, we obtain

Denote

Then, by matching coefficients, we can derive the following three ordinary differential equations:

According to the boundary condition $D(0, u) = 0$, we obtain

Equation $(3.5)$ is a Riccati equation, whose general solution can be presented as

where

From Eq $(3.5)$, we can rewrite Eq (3.7) as

Integrating both sides of the above-mentioned equation, we have the following result:

where

Combining (13) and the above two expressions, Theorem 3.1 follows.

3.3.   Option pricing in the bid-ask price economy

According to conic finance theory, there are two types of market prices. The minimal acceptable price for selling a claim $X$ is known as the ask price $a_{\gamma}(X)$, while the maximal acceptable price for purchasing the claim is referred to as the bid price $b_{\gamma}(X)$. Then, the following lemma should be introduced. For more detailed information, one can refer to ^[1,2,3].

Lemma 3.1. If we assume that $X$ represents the cashflow of the claim at its expiry time $T$, the ask and bid prices of the claim are respectively determined by

and

where $r$ is the risk-free interest rate, $\Psi$ is the probability distortion function, $\gamma$ denotes the degree of the distortion, and $F_{X}(x)$ is the distribution of random $X$. In addition, according to Lemma 3.1, the ask and bid prices for European call options will be shown in the following theorem.

Theorem 3.2. Given the MINMAXVAR distortion function, i.e., $\Psi_{\gamma}(w) = 1-\left(1-w^{\frac{1}{1+\gamma}}\right)^{1+\gamma}, \quad\gamma\geqslant 0, \quad w \in[0, 1]$, it follows from Eqs (2.1)–(2.3) that the ask and bid prices for European call options can be expressed as

where

with $\widetilde{A}_{0} = \frac{2}{b-a}\operatorname{Re}\{\Phi(0)\}$. $K$ denotes the strike price of the option, while $a_{\gamma}(C)$ and $b_{\gamma}(C)$ represent the ask and bid prices for the option, respectively.

Proof. We can reconstruct the density function $f_{\ln S}$ of the random variable, i.e., the log-asset price $\ln S_T$ with a truncated region $[a, b]$ using the results provided by Fang and Oosterlee ^[25] as

where $\Phi\{\cdot\}$ is the characteristic function of the density function $f(x)$, which is obtained in Theorem 3.1, and $\operatorname{Re}\{\cdot\}$ is an operator to take the real part. It should be noted that the prime of $\sum^{\prime}$ is used to alert that the first term of the summation should be multiplied by a weight of $1/2$.

Meanwhile, let

where $\Phi\{\cdot\}$ is the characteristic function of the density function $f(x)$ which is obtained in Theorem 3.1. Then, the distribution function $F_{\ln S}$ of the log-asset price $\ln S_T$ is given by

As a result, the distribution function $F_S$ of the underlying price $S_T$ can be directly obtained through

By employing Lemma 3.1 and the MINMAXVAR distortion function, European call ask and bid prices can be respectively expressed as

where the formula of $F_{S}(x)$ is provided in Eq (3.9).

Following Theorem 3.2, we can obtain the ask and bid prices of European put options as the following corollary.

Corollary 3.1 If the European put option is also controlled by the stochastic differential equations (2.1)–(2.3), given the same distortion function in Theorem 3.2 and utilizing a similar derivation, the ask and bid prices of European put options can be expressed as

It should be remarked that the integrals involved in ask and bid prices given by Theorem 3.2 and Corollary 3.1 can be numerically computed with the Gaussian quadrature, which is one of the best quadratures with a high degree of accuracy and efficiency. This can be very easily implemented by the MATLAB built-in function $quadgk$, which ensures the speed for the implementation of the proposed approximation method when calculating ask and bid prices.

4.   Numerical analysis

4.1.   Accuracy analysis

In this subsection, we benchmark our results by making use of a Monte Carlo (hereafter, MC) simulation. All the computation is implemented using MATLAB 2016a on a computer equipped with an Intel Core $i3$ CPU @ $2.53$ GHz.

First, we follow Fang and Oosterlee [18, 19] to select $[a, b]$ as

when computing our approximation formula, where $x_{0} = \ln S_{0}, L = 10$, and $c_{n}$ is the $n$-th cumulant of $lnS_{T}$.

Moreover, each sample path in the MC simulation is generated with the time interval being uniformly divided containing $M_1 = 252$ points, and we use $M_2 = 100,000$ as the number of sample paths. If the payoff produced by the $i$-th sample path is denoted by payoff$(i)$, $i = 1, 2, \cdots, M_2$, we can obtain

and

In addition, we utilize a discrete scheme of stochastic differential equations (2.1)–(2.3) as follows:

where $\Delta t = \frac{T}{252}$, $\varepsilon_{i} \sim N(0, 1), i = 1, \ldots, 3$, and $\varepsilon_{1}, \varepsilon_{2}, \varepsilon_{3}$ are independent of each other. $N_{(t)}$ is a Poisson process with parameter $\lambda$. Then,

where $S_{i, T}$ is the price of the underlying asset which is produced by the $i$-th sample path at the maturity time $T$.

With the values of other parameters provided in Table 1, bid and ask prices of European call options are shown in Tables 2 and 3 with different values of strike prices, respectively. The absolute relative error (Abs.R.E.) is defined as

One should note that the results derived with our analytical approximation and those from the MC simulation are close to each other, which is a clear indication that the approximation is of high accuracy. Meanwhile, Figure 1 shows how the bid-ask prices of European call options change with different sample paths $M_2$. It shows that the results of MC simulation usually fluctuate within a confidence interval. It also indicates that the results of MC simulation and our approach get close to each other as the number of paths increases.

4.2.   Sensitivity analysis

In this subsection, the impact of parameter changes on bid and ask prices will be investigated, with a focus on: ($i$) The long-run mean level $\gamma$; ($ii$) The jump intensity $\lambda$ and time to maturity $T-t$; ($iii$) The long-run mean level $\theta_v$ and $\theta_r$.

The effect of different $\gamma$ on the bid-ask prices of European call options are shown in Figure 2. We find that, as the market liquidity indicator $\gamma$ increases, or equivalently the liquidity of the market becomes lower, the bid-ask spread becomes larger and cannot be ignored. As $\gamma$ approaches $0$, the bid and ask prices will converge to a single price. This shows it is reasonable to consider the market liquidity in option pricing.

Figure 3 demonstrates how European option prices are affected by different $\lambda$ and $T-t$. Both bid and ask prices increase with the time to maturity, which is consistent with financial intuition since longer time implies larger time values of options. Larger jump intensity also contributes to greater option prices, since the underlying asset price becomes more volatile when there are possibly more jumps, which leads to higher risks and larger option premiums.

What is plotted in Figure 4 is the influence of $\theta_v$ and $\theta_r$ on European option prices. It is clear that both $\theta_v$ and $\theta_r$ have a positive impact on European call option prices. This is reasonable since bigger $\theta_v$ means that the underlying asset prices are more volatile in the long run. Also, bigger $\theta_r$ indicates that the expected return of the underlying asset under the risk-neutral world is larger. Both will lead to higher option premiums.

5.   Conclusions

This article aims to develop a pricing framework for European options within a hybrid GARCH-Vasicek model incorporating double exponential jumps in the bid-ask price economy. We derive analytical formulas for calculating bid and ask prices of European options utilizing the COS method, employing an approximation approach to obtain the characteristic function. We have also shown how both bid and ask prices vary with different values of parameters.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This work is supported by the Creative Research Groups of the National Natural Science Foundation of Guizhou of China under Grants Qian Jiao HE KY ZI (No. 2019069), the National Natural Science Foundation of China (No. 12101554) and the Fundamental Research Funds for Zhejiang Provincial Universities (No. GB202103001).

The authors would like to gratefully acknowledge the anonymous referees' constructive comments and suggestions, which greatly help to improve the readability and clarity of the manuscript.

Conflict of interest

We confirm that there are no known conflicts of interest associated with this publication, and there has been no significant financial support for this work that could have influenced its outcome.