An efficient algorithm to solve the geometric Asian power option price PDE under the stochastic volatility model

Abstract

This study focuses on the valuation of geometric Asian power options and presents an efficient numerical algorithm for solving the option price PDE. The analytical methodology utilizes the fractional Ito formula and replicating portfolio techniques to derive a PDE that characterizes the option price. However, due to the lack of an analytical solution for this PDE, a numerical method is proposed to solve it. The numerical solution involves implementing a time-semi-discrete scheme obtained through forward time difference approximation, while the other derivatives in the equation are approximated using cubic B-spline quasi-interpolation approximation. By employing the respective scheme and incorporating the initial and boundary conditions, the numerical solution for the equation is obtained. Subsequently, the stability of the method is investigated, and numerical results are presented. The main advantages of the presented method are its simplicity for computer implementation and its suitability for multi-dimensional problems.

Introduction

Power options belong to the class of exotic options, which exhibit distinctive features and payoffs compared to traditional options. These options derive their name from the fact that the payoff at maturity is determined by raising the underlying asset price to a positive power, denoted as “n.” In other words, the payoff of a power option is a function of the n-th power of the asset price at expiration. The use of power functions in option payoffs introduces additional complexity and uniqueness to these financial instruments. Traditional options, such as European or American options, have linear or non-linear payoffs based on the asset price alone or in combination with other factors like time or volatility. However, power options introduce a distinct non-linearity that can significantly impact the potential gains or losses for option holders.

In recent decades, researchers have used the combination of this option and different types of Asian options to present powerful financial tools to investors. In 2012, Bin Peng and Fei Peng proposed a general framework to calculate the Asian power option under the jump-fraction model. To do this, they applied the Ito formula and a self-financing strategy and derived a partial differential equation (PDE) in the fractional environment with jumps. Then, under specific boundary conditions, they obtained an analytical solution for the PDE to calculate the option price at any time before the maturity time [1]. After that, Jerim Kim et al. used the stochastic Heston model to predict the asset price and derived a semi-analytical formula to compute the price of the power Asian option under the model. They also introduced explicit formulas for the moment-generating function of log price and the pricing of power options, assuming finite moments. However, when the moments become infinite, they presented two numerical methods to calculate the prices of power put and capped power call options [2]. Youngrok Lee et al. considered different versions of the power option and obtained prices of them under the stochastic volatility model. Indeed, they used the stochastic model that presented by Chöbel and Zhu and achieved the option price formulas for different kinds of the power option. In addition, they calculated the effect of changing the values of the model parameters on the option price [3]. After that researchers introduced the geometric Asian power option as a new version of the power option. In 2018, B.L.S. Prakasa Rao studied about the long memory version of the geometric Brownian motion (gBm) model and applied the mixed fractional gBm model for $H\in (\frac{3}{4},1)$ to forecast the stock price. Then, they evaluated the geometric Asian power option and found a closed-form formula to obtain the option price [4]. Sunday O. Edeki et al. studied the Greek parameters concept in option pricing and derived the Greeks of the continuous arithmetic Asian option pricing model [5].

Wei-Guo Zhang et al. presented an analytical formula to calculate the geometric Asian power option under different states of the strike price where they used the mixed gBm model to predict the stock price. Then, since the price of the underlying asset changes under various factors, it is very difficult to calculate the exact price for derivatives. To solve this problem, they expressed the price of the geometric Asian power option as an interval [6]. In 2018, Foad Shokrollahi computed the price of the geometric Asian power option where the underlying asset follows the time-changed mixed fractional Brownian motion. Then, they considered the asset price where the asset pays constant dividends and obtained the option price [7]. Sunday O. Edekl et al. studied the solution of the Asian option and presented a novel semi-analytical approach called the Projected Differential Transform Method (PDTM) to obtain the analytical solution of a continuous arithmetic Asian option [8]. After that, in 2019, Sunday O. Edekl et al. presented a novel application of the Adomian Decomposition Method (ADM), which is a semi-analytical technique, for the analytical solution of a continuous arithmetic Asian option model [9].

In 2021, Wei Wang et al. studied about the sub-mixed fractional Brownian motion process and calculated the geometric Asian power option price by using the sub-fractional Ito lemma [10]. Gifty Malhotra et al. studied the pricing of continuous geometric Asian options (GAOs) within a novel multifactor stochastic volatility model [11]. in 2022, D. Ahmadian et al. presented numerical experiments encompassing arithmetic Asian rainbow options on two, three, four, and ten underlying assets. The outcomes of these experiments shed light on the effectiveness of the proposed control variate technique in reducing the variance of the Monte Carlo estimator [12]. Christian-Oliver Ewald et al. investigated the pricing of Asian options on commodity futures contracts, taking into account the effects of stochastic convenience yield, stochastic interest rates, and jumps in the commodity spot price [13].

To calculate the option price, if the option price cannot be obtained as an analytical formula, researchers use numerical and simulation methods [14, 15–16]. In many cases, high-dimensional differential equations do not have analytical solutions. Therefore, researchers primarily resort to numerical methods such as the finite difference method and B-spline method to solve these equations [17, 18, 19, 20, 21–22]. The main drawback of using the finite difference method in dimensions higher than two is its reduced accuracy and the high complexity of its computer implementation. Consequently, this method becomes impractical for solving problems with higher dimensions. Alternatively, some authors have employed wavelets to tackle problems with dimensions higher than two [23]. Despite its good accuracy, using wavelets comes with the disadvantage of high computational cost and relatively complex computer implementation. As a result, solving a 4-dimensional problem sometimes requires 2 weeks or more to run a computer program. To address this, we propose a method that is simple to implement, offers good accuracy, and can be easily applied to problems with higher dimensions.

In this study, we apply the mixed fractional Heston model to forecast the stock price and calculate the geometric Asian option pricing under the proposed model. By considering that the structure of the proposed model is complex, so one can not find the closed-form solution to calculate the option price. To obtain the option price, we use the delta Hedging strategy and the fractional Ito Lemma, to find the option price PDE. Then, to solve this PDE, we employ a time-semi-discrete scheme obtained through forward time difference approximation, while the other derivatives in the equation are approximated using the cubic B-spline quasi-interpolation approximation. By using this scheme and incorporating the initial and boundary conditions, we obtain the numerical solution for the option price. After that, we investigate the stability of our method and present some numerical results to demonstrate its accuracy and efficiency.

The presented numerical method is based on the approximation of the time derivative using finite differences and the approximation of other derivatives in the equation using quasi-cubic B-spline interpolation. In the application of cubic B-spline quasi-interpolation, the multiplications of specific and sparse matrices in known vectors are used, and the approximations of the derivatives in the equation are simply obtained in vector form. Therefore, this method can be easily implemented in the computer and a simple algorithm is obtained for problems with more than two independent variables. In addition, since this is an explicit method, the execution speed of its algorithm is also acceptable. On the other hand, the challenge that this numerical method faces is that in some problems, the stability range of the method may be small. That is, in this case, it will be necessary for N to be a relatively large number.

The main contribution of our study is to provide a simple and effective numerical method for pricing geometric Asian power options under the mixed fractional stochastic volatility model. Our method has several advantages over the existing methods. First, it does not require any transformation or discretization of the underlying asset price process, which preserves the original stochastic properties of the model. Second, it does not involve any complex integration or inversion techniques, which reduces the computational complexity and error. Third, it can be easily extended to multi-dimensional problems, such as basket or spread options, by using B-splines. Our study also has important implications for practitioners and researchers in the field of financial engineering, as it provides a reliable and robust tool for valuing and hedging geometric Asian power options and other exotic options under realistic market conditions.

The remainder of this paper is organized as follows. In Section 2, we use the fractional Ito formula and the delta hedging strategy to obtain the geometric Asian power option price PDE under the mixed fractional Heston model. In Section 3, we apply the finite difference method for the time discretization of the option price PDE and utilize the cubic B-spline quasi-interpolation approximation to obtain the rest of other derivatives needed to solve the equation. In Section 4, we study the stability and convergence analysis. Finally, in Section 5, we use the presented methods and calculate the option price.

Option price PDE under the mixed fractional Heston model

In this section, we obtain the geometric Asian power option price PDE under the mixed fractional Heston model. In general, the practical implications of the geometric Asian power option in line with fractional derivative theory are significant for both financial practitioners and researchers. By incorporating fractional derivative models into the pricing and risk management of Asian power options, several practical advantages emerge. Firstly, the inclusion of fractional calculus allows for a more accurate representation of the underlying asset’s price dynamics, capturing long-range dependence and memory effects commonly observed in financial markets. This advancement enables more precise pricing and hedging strategies, leading to improved risk management for financial institutions and investors. Additionally, the utilization of fractional derivative models facilitates a better understanding of the option’s sensitivity to changes in market conditions, allowing for enhanced decision-making in a dynamic trading environment. Moreover, the practical implications extend to the development of investment strategies and the optimization of portfolios. Incorporating geometric Asian power options based on fractional derivative theory enables investors to diversify their portfolios effectively and hedge against market volatility with greater precision. The practical relevance of this research lies in its potential to enhance the efficiency and effectiveness of derivative pricing, risk management, and portfolio optimization techniques in the context of geometric Asian power options, ultimately supporting more informed and profitable investment decisions in the financial industry.

Combining a long memory and stochastic volatility properties and use them in financial model can be a powerful tool for asset price prediction. The long memory component captures persistent patterns and dependencies in the data, while the stochastic volatility component takes into account the variability of asset returns over time. This combination allows for a more comprehensive modeling of asset price dynamics, potentially leading to improved predictions. Let $(\mathrm{\Omega },\mathcal{F},P)$ be a probability space and the stock price follows the following model:

1

$$\begin{array}{cc}& d{S}_t=\mu S_{t}dt+\sqrt{V_t}{S}_{t}d{M}_t^H,\hfill \\ \hfill & d{V}_t=\kappa (\theta -V_t)dt+\sigma \sqrt{V_t}d{N}_t^H,\hfill \\ \hfill & d{M}_t^{H}d{N}_t^H=\rho \left(d,t,+,H,t^{2H-1},d,t\right),\hfill \end{array}$$where ${\left\{M_t^H\right\}}_{t\ge 0}$ and ${\left\{N_t^H\right\}}_{t\ge 0}$ for $H\in (\frac{3}{4},1)$ are two mixed fractional Brownian motion processes with correlation $\rho \in (-1,1)$. In the following based on the model, we derive a PDE to obtain the option price. In general, researchers regularly employ risk-neutral valuation as a fundamental technique to derive the PDE governing option pricing. This approach offers a robust framework by eliminating the influence of stochastic and constructing a risk-free portfolio. The foundation of risk-neutral valuation rests on the assumption of a risk-neutral world, where investors are assumed to be devoid of risk aversion. In the constant volatility model, the random component of the option value is solely determined by the underlying asset, which is assumed to be traded in the market. Consequently, it is feasible to hedge the option and achieve market completeness. Conversely, in the stochastic volatility setting, the valuation of a contingent claim relies on both the randomness exhibited by the underlying asset and the inherent variability of the asset’s return volatility. To derive the PDE for option pricing under a stochastic volatility model using a portfolio approach, we can follow the steps outlined below:

1. Construct a portfolio containing an option, a share of the underlying stock, and a share of the volatility process. The weights of the stock and volatility positions in the portfolio are denoted as $\mathrm{{\rm Y}}$ and $\mathrm{\Phi }$, respectively.

2. Apply the fractional Ito’s lemma to the portfolio value and set up a risk-neutral pricing framework.

3. Define a riskless portfolio by eliminating the sources of randomness (i.e., the stochastic terms) from the portfolio value expression obtained in the previous step. This ensures that the portfolio value changes only due to the risk-free rate.

4. Employ self-financing conditions to derive the PDE for the portfolio value. This involves setting up differential equations involving the portfolio value, stock price, and volatility process.

5. Focus on the option value component of the portfolio and derive the PDE for option pricing by eliminating the stock and volatility terms from the portfolio value equation.

Theorem 1

(Fractional Ito Formula) Let $f(t,x):[0,\infty )\times \mathbb{R}⟶\mathbb{R}$ belongs to ${\mathcal{C}}^{1,2}\left(0,,,\infty ,),\times ,\mathbb{R}\right)$ and let the stochastic process $Y_t$ follows the equation $\left(\frac{1}{2},<,H,<,1\right)$ [24]

2

$$\begin{array}{c}\hfill dY(t)=\upsilon (t,Y_t)dt+\nu (t,Y_t)d{B}_t^H.\end{array}$$Then, we have

3

$$\begin{array}{c}\hfill df(t,Y_t)=(\frac{\partial f}{\partial t}(t,Y_t)+\frac{\partial f}{\partial y}(t,Y_t)\upsilon (t,Y_t)+\frac{1}{2}\frac{{\partial }^{2}f}{\partial y^2}(t,Y_t)b^2(t,Y_t))dt+\frac{\partial f}{\partial y}(t,Y_t)d{B}_t^H,\end{array}$$ where

4

$$\begin{array}{c}\hfill b^2(t,y)=2H(2H-1){\int }_0^t{\left(t,-,s\right)}^{2H-2}{\nu }^2(s,y)ds.\end{array}$$

Lemma 1

Assume $C=C(t,S_t,V_t,K_t)$ belongs to ${\mathbb{C}}^{1,2,2,1}([0,\infty )\times \mathbb{R}\times \mathbb{R})$ at each time $t\ge 0$. Then, we have

5

$$\begin{array}{cc}\hfill \frac{\partial C}{\partial t}& +r{S}_t\frac{\partial C}{\partial S_t}+r{V}_t\frac{\partial C}{\partial V_t}+K_t\frac{n^{\ast }ln(S_t)-ln(K_t)}{t}r{S}_t\frac{\partial C}{\partial K_t}\hfill \\ \hfill & +\left(\frac{1}{2},+,H,t^{2H-1}\right)\left(V_t,S_t^2,\frac{{\partial }^{2}C}{\partial {S_t}^2},+,{\sigma }^2,V_t,\frac{{\partial }^{2}C}{\partial {V_t}^2}\right)\hfill \\ \hfill & +\left(1,+,2,H,t^{2H-1}\right)\rho \sigma V_t{S}_t\frac{{\partial }^{2}C}{\partial S_t\partial V_t}-rC=0,\hfill \end{array}$$where $\frac{3}{4}<H<1$.

Proof

The Payoff function of the geometric Asian power call and put option are as follows:

6

$$\begin{array}{c}\hfill \mathcal{C}(t,S_t,K_t)=max\left\{S_t^{n^{\ast }},-,K_t,,,0\right\},\end{array}$$and

7

$$\begin{array}{c}\hfill P(t,S_t,K_t)=max\left\{K_t,-,S_t^{n^{\ast }},,,0\right\},\end{array}$$where for any $t\ge 0$

8

$$\begin{array}{c}\hfill K_t=exp\left\{\frac{n^{\ast }}{t},{\int }_0^t,ln,S_u,d,u\right\},\end{array}$$is the floating strike price and $n^{\ast }\in \left\{1,,,2,,,\cdots \right\}$ is the power index. We consider a portfolio X contains a geometric Asian power option $C(t,S_t,V_t,K_t)$, $\mathrm{\Phi }$ shares of volatility asset V, and $\mathrm{{\rm Y}}$ shares of the asset S. Thus, the value of the portfolio at time t is

9

$$\begin{array}{c}\hfill X_t=C(t,S_t,V_t,K_t)-{\mathrm{\Phi }}_t{V}_t-{\mathrm{{\rm Y}}}_t{S}_t.\end{array}$$After the time interval $\delta t$, the evaluation of the portfolio $X_t$ at time t is given by

10

$$\begin{array}{c}\hfill \delta X_t=\delta C(t,S_t,V_t,K_t)-{\mathrm{{\rm Y}}}_t\delta S_t-{\mathrm{\Phi }}_t\delta V_t,\end{array}$$where $\delta C$, $\delta S_t$, and $\delta V_t$ are the change of the option price, the asset price, and the volatility amount, respectively. By using fractional Ito lemma for (10), we can write

11

$$\begin{array}{cc}\hfill dX(t)& =[\frac{\partial \mathcal{C}}{\partial t}+\mu S_t\frac{\partial \mathcal{C}}{\partial S_t}+\left(\kappa ,-,\theta ,V,(,t,)\right)\frac{\partial \mathcal{C}}{\partial V_t}+\frac{\partial \mathcal{C}}{\partial K_t}\frac{d{K}_t}{\mathit{dt}}\hfill \\ \hfill & +\left(\frac{1}{2},+,H,t^{2H-1}\right)\left(V,S^2,\frac{{\partial }^{2}C}{\partial {S_t}^2},+,{\sigma }^2,V,\frac{{\partial }^{2}C}{\partial {V_t}^2}\right)\hfill \\ \hfill & +\left(1,+,2,H,t^{2H-1}\right)\rho \sigma V_t{S}_t\frac{{\partial }^{2}C}{\partial S_t\partial V_t}-\mu S_t{\mathrm{{\rm Y}}}_t-\left(\kappa ,-,\theta ,V_t\right){\mathrm{\Phi }}_t]dt\hfill \\ \hfill & +\sqrt{V_t}{S}_t\left(\frac{\partial \mathcal{C}}{\partial S_t},-,{\mathrm{{\rm Y}}}_t\right)d{M}_t^H+\sigma \sqrt{V_t}\left(\frac{\partial \mathcal{C}}{\partial V_t},-,{\mathrm{\Phi }}_t\right)d{N}_t^H.\hfill \end{array}$$To achieve the option price PDE, the random parts of the (11) must be removed. To do so, if ${\mathrm{{\rm Y}}}_t=\frac{\partial \mathcal{C}}{\partial S_t}$ and ${\mathrm{\Phi }}_t=\frac{\partial \mathcal{C}}{\partial V_t}$, we obtain

12

$$\begin{array}{cc}\hfill d{X}_t& =[\frac{\partial \mathcal{C}}{\partial t}+\frac{\partial \mathcal{C}}{\partial K_t}\frac{d{K}_t}{\mathit{dt}}+\left(\frac{1}{2},+,H,t^{2H-1}\right)\left(V_t,S_t^2,\frac{{\partial }^{2}C}{\partial {S_t}^2},+,{\sigma }^2,V_t,\frac{{\partial }^{2}C}{\partial {V_t}^2}\right)\hfill \\ \hfill & +\left(1,+,2,H,t^{2H-1}\right)\rho \sigma V_t{S}_t\frac{{\partial }^{2}C}{\partial S_t\partial V_t}]dt.\hfill \end{array}$$The return of a risk-free portfolio is expected to equal the interest rate r due to arbitrage considerations, the risk-free nature of the portfolio, and the opportunity cost of investing in alternative assets. Therefore, we have

13

$$\begin{array}{c}\hfill d{X}_t=r{X}_{t}dt=r\left({\mathcal{C}}_t,-,{\mathrm{{\rm Y}}}_t,S_t,-,{\mathrm{\Phi }}_t,V_t\right)dt=r\left({\mathcal{C}}_t,-,\frac{\partial \mathcal{C}}{\partial S_t},S_t,-,\frac{\partial \mathcal{C}}{\partial V_t},V_t\right)dt.\end{array}$$Thus,$$\begin{array}{cc}& [\frac{\partial \mathcal{C}}{\partial t}+\mu S_t\frac{\partial \mathcal{C}}{\partial S_t}+\left(\kappa ,-,\theta ,V_t\right)\frac{\partial \mathcal{C}}{\partial V_t}+\frac{\partial \mathcal{C}}{\partial K_t}\frac{d{K}_t}{\mathit{dt}}+\left(\frac{1}{2},+,H,t^{2H-1}\right)\left(V_t,S_t^2,\frac{{\partial }^{2}C}{\partial {S_t}^2},+,{\sigma }^2,V_t,\frac{{\partial }^{2}C}{\partial {V_t}^2}\right)\hfill \\ \hfill & +\left(1,+,2,H,t^{2H-1}\right)\rho \sigma V_t{S}_t\frac{{\partial }^{2}C}{\partial S_t\partial V_t}-\mu S_t{\mathrm{{\rm Y}}}_t-\left(\kappa ,-,\theta ,V_t\right){\mathrm{\Phi }}_t]dt=r\left({\mathcal{C}}_t,-,\frac{\partial \mathcal{C}}{\partial S_t},S_t,-,\frac{\partial \mathcal{C}}{\partial V_t},V_t\right)dt.\hfill \end{array}$$Therefore, the final version of the PDE is derived as follows:$$\begin{array}{cc}\hfill \frac{\partial C}{\partial t}& +r{S}_t\frac{\partial C}{\partial S_t}+r{V}_t\frac{\partial C}{\partial V_t}+K_t\frac{n^{\ast }ln(S_t)-ln(K_t)}{t}r{S}_t\frac{\partial C}{\partial K_t}\hfill \\ \hfill & +\left(\frac{1}{2},+,H,t^{2H-1}\right)\left(V_t,S_t^2,\frac{{\partial }^{2}C}{\partial {S_t}^2},+,{\sigma }^2,V_t,\frac{{\partial }^{2}C}{\partial {V_t}^2}\right)\hfill \\ \hfill & +\left(1,+,2,H,t^{2H-1}\right)\rho \sigma V_t{S}_t\frac{{\partial }^{2}C}{\partial S_t\partial V_t}-rC=0.\hfill \end{array}$$$\square$

Numerical method

In this section, the finite difference method is employed for the time discretization of the partial differential (5). Subsequently, the scheme will be complemented by utilizing the cubic B-spline quasi-interpolation approximation.

Let $\mathrm{\Gamma }=[S_{min},S_{max}]\times [V_{min},V_{max}]\times [K_{min},K_{max}]$ and define ${\mathrm{\Omega }}_T=[0,T)\times \mathrm{\Gamma }$. We consider a class of fully nonlinear parabolic PDEs of the form

14

$$\begin{array}{cc}& C_t+rS{C}_S+rV{C}_V+f(S,K,t)C_K+g(t)\left(V{S}^2{C}_{\mathit{SS}}+{\sigma }^{2}V{C}_{\mathit{VV}}\right)\hfill \\ \hfill & +h(S,V,t)C_{\mathit{SV}}-rC=0,(t,x)\in {\mathrm{\Omega }}_T,\hfill \\ \hfill & C(t,x)=\mathbf{g}(t,x),(t,x)\in (0,T)\times \partial \mathrm{\Gamma },\hfill \\ \hfill & C(T,x)=e^{-rT}max(S^{n^{\ast }}-K,0),x\in \mathrm{\Gamma },\hfill \end{array}$$where $f(S,K,t)=K\frac{n^{\ast }lnS-lnK}{t}rS$, $g(t)=H{t}^{2H-1}+\frac{1}{2}$, $h(S,V,t)=\left(1+2H{t}^{2H-1}\right)\rho \sigma SV$ and $x=(S,V,K)$. Also, the Dirichlet boundary conditions are defined as follows$$\begin{array}{cc}& C(t,S_{min},V,K)=0,C(t,S_{max},V,K)=e^{-r(T-t)}max(S_{max}-K,0),\hfill \\ \hfill & C(t,S,V_{min},K)=e^{-r(T-t)}max\left(S_{\mathit{min}}^{n^{\ast }}\left(e^{rt{n}^{\ast }}-e^{\frac{1}{2}rt{n}^{\ast }}\right),0\right),\hfill \\ \hfill & C(t,S,V_{max},K)=e^{-r(T-t)}max(S_{max}-K,0),\hfill \\ \hfill & C(t,S,V,K_{min})=e^{-r(T-t)}max(S-K_{min},0),C(t,S,V,K_{max})=0.\hfill \end{array}$$We use the change of variable $\tau =T-t$ and the assumption $C(T-t,S,V,K)=C(\tau ,S,V,K)$ to derive the following initial-boundary value problem

15

$$\begin{array}{cc}& C_{\tau }=rS{C}_S+rV{C}_V+f(S,K,t)C_K+g(t)\left(V{S}^2{C}_{\mathit{SS}}+{\sigma }^{2}V{C}_{\mathit{VV}}\right)\hfill \\ \hfill & +h(S,V,t)C_{\mathit{SV}}-rC,(\tau ,x)\in (0,T]\times \mathrm{\Gamma },\hfill \\ \hfill & C(\tau ,x)=\mathbf{g}(T-\tau ,x),(\tau ,x)\in (0,T)\times \partial \mathrm{\Gamma },\hfill \\ \hfill & C(0,x)=e^{-rT}max(S^{n^{\ast }}-K,0),x\in \mathrm{\Gamma },\hfill \end{array}$$We set ${\tau }_n=n\mathrm{\Delta }\tau$, $n=0,1,\cdots ,N$ where $\mathrm{\Delta }\tau =\frac{T}{N}$. Now by assumption $C({\tau }_n,S,V,K)=C^n(S,V,K)$ and using the simple forward finite difference, we obtain the following temporal semi-discrete scheme

16

$$\begin{array}{cc}\hfill C^{n+1}=C^n& +\mathrm{\Delta }\tau (rS{C}_S^n+rV{C}_V^n+f(S,K,t_n)C_K^n\hfill \\ \hfill & +g(t_n)\left(V{S}^2{C}_{\mathit{SS}}^n+{\sigma }^{2}V{C}_{\mathit{VV}}^n\right)\hfill \\ \hfill & +h(S,V,t_n)C_{\mathit{SV}}^n-r{C}^n),\hfill \end{array}$$where $t_n=T-{\tau }_n$.

Assume that an interval $I=[S_{min},S_{max}]$ is given, and let $\mathcal{B}({\mathcal{S}}_{M_S})$ be the space of cubic splines on uniform partition ${\mathcal{S}}_{M_S}=\{S_{m_s}=S_{min}+m_s\mathrm{\Delta }S,m_s=0,1,\cdots ,M_S\}$ with meshlength $\mathrm{\Delta }S=\frac{S_{max}-S_{min}}{M_S}$. Also, let a basis of $\mathcal{B}({\mathcal{S}}_{M_S})$ be $\{B_{j,\mathbf{r}},j=1,2,\cdots ,M_S+3\}$ where $B_{j,\mathbf{r}}$ is the jth cubic B-spline for the knot sequence $\mathbf{r}={\left(r_i\right)}_{i=-3}^{M_S+3}$ where $r_{-3}=r_{-2}=r_{-1}=S_{min}$, $S_{max}=r_{M_S+1}=r_{M_S+2}=r_{M_S+3}$ and $r_i=S_i$, $0\le i\le M_S$. We use $B_j$ instead of $B_{j,\mathbf{r}}$ which is defined using the recursive relation [25]:

17

$$\begin{array}{c}\hfill b_{j,p}(S)=\frac{S-r_j}{r_{j+p}-r_j}{b}_{j,p-1}(S)+\frac{r_{j+p+1}-S}{r_{j+p+1}-r_{j+1}}{b}_{j+1,p-1}(S),\end{array}$$starting from$$\begin{array}{c}\hfill b_{j,0}(S)=\left\{\begin{array}{c}1,S_j\le S<S_{j+1},\hfill \\ 0,\mathrm{otherwise},\hfill \end{array}\right)\end{array}$$where $B_k(S)=b_{k-4,3}(S)$, $k=1,2,\cdots ,M_S+3$, and under the convention that fractions with zero denominator have value zero. By considering the above definition, all the B-splines take the value zero at the end point $S_{max}$. Therefore, in order to avoid asymmetry over the interval $[S_{min},S_{max}]$, it is common to assume the B-splines to be left continuous at $S_{max}$. With these notations, the support of $B_j$ is $supp(B_j)=[r_{j-4},r_j]$.

The univariate cubic spline quasi-interpolants (abbreviation QIs [26]) can be defined as operators of the form

18

$$\begin{array}{c}\hfill Qf=\sum _{j=1}^{M_S+3}{\mu }_j(f)B_j,\end{array}$$where the coefficient functionals ${\mu }_j$ are

19

$$\begin{array}{cc}& {\mu }_1(f)=f_0,{\mu }_2(f)=\frac{1}{18}\left(7,f_0,+,18,f_1,-,9,f_2,+,2,f_3\right),\hfill \\ \hfill & {\mu }_j(f)=\frac{1}{6}\left(-,f_{j-3},+,8,f_{j-2},-,f_{j-1}\right),3\le j\le M_S+1,\hfill \\ \hfill & {\mu }_{M_S+2}=\frac{1}{18}\left(2,f_{M_S-3},-,9,f_{M_S-2},+,18,f_{M_S-1},+,7,f_{M_S}\right),{\mu }_{M_S+3}=f_{M_S}.\hfill \end{array}$$Let ${\mathrm{\Pi }}_3$ denotes the space of polynomials of total degree at most 3. Then, for any subinterval $I_k=[r_{k-1},r_k]$, $1\le k\le M_S$, and any function f [27],$$\begin{array}{c}\hfill {∥f,-,Q,f∥}_{\infty ,I_k}\le \left(1,+,{∥Q∥}_{\infty }\right)w_{\infty ,I_k}\left(f,,,{\mathrm{\Pi }}_3\right),\end{array}$$where the distance of f to polynomials is defined by $w_{\infty ,I_k}\left(f,,,{\mathrm{\Pi }}_3\right)=inf$$\left\{{∥f,-,p∥}_{\infty ,I_k},,,p,\in ,{\mathrm{\Pi }}_3\right\}$. Here ${∥f,-,p∥}_{\infty ,I_k}=\underset{r\in I_k}{max}\left|f(r)-p(r)\right|$. As a consequence of this property, the approximation order is $O\left({(\mathrm{\Delta }S)}^4\right)$ on smooth functions. In [27], one can find that for $f\in C^4[S_{min},S_{max}]$

20

$$\begin{array}{c}\hfill {∥f^{(k)},-,{(Qf)}^{(k)}∥}_{\infty }=O\left({(\mathrm{\Delta }S)}^{4-k}\right),k=0,1,2.\end{array}$$For approximate derivatives of f by derivatives of Qf up to the order ${\left(\mathrm{\Delta },S\right)}^3$ at $S_{m_s}$, we can evaluate the value of f at $S_{m_s}$ by ${(Qf(S_{m_s}))}^{\prime }={\sum }_{j=1}^{M_S+3}{\mu }_j(f)B_j^{\prime }(S_{m_s})$ and ${(Qf(S_{m_s}))}^{\prime \prime }={\sum }_{j=1}^{M_S+3}{\mu }_j(f)B_j^{\prime \prime }(S_{m_s})$. We set $U={\left(f_0,,,f_1,\cdots ,,,f_{M_S}\right)}^T$, $U_S={\left(f_0^{\prime },,,f_1^{\prime },\cdots ,,,f_{M_S}^{\prime }\right)}^T$ and $U_{\mathit{SS}}={\left(f_0^{\prime \prime },,,f_1^{\prime \prime },\cdots ,,,f_{M_S}^{\prime \prime }\right)}^T$, where $f_{m_s}^{\prime }={\left(Q,f_{m_s}\right)}^{\prime }$ and $f_{m_s}^{\prime \prime }={\left(Q,f_{m_s}\right)}^{\prime \prime }$. By solution of the linear system:

21

$$\begin{array}{cc}\hfill {(Qf(S_{m_s}))}^{\prime }=& \sum _{j=1}^{M_S+3}{\mu }_j(f)B_j^{\prime }(S_{m_s}),m_s=0,1,\cdots ,M_S,\hfill \end{array}$$

22

$$\begin{array}{cc}\hfill {(Qf(S_{m_s}))}^{\prime \prime }=& \sum _{j=1}^{M_S+3}{\mu }_j(f)B_j^{\prime \prime }(S_{m_s}),m_s=0,1,\cdots ,M_S,\hfill \end{array}$$we obtain $U_S=D_{M_S}^{1}U$ and $U_{\mathit{SS}}=D_{M_S}^{2}U$ where

23

$$\begin{array}{cc}& D_{M_S}^1=\frac{1}{\mathrm{\Delta }S}{\left(\begin{array}{ccccccccc}-\frac{11}{6}& 3& -\frac{3}{2}& \frac{1}{3}& & & & & \\ -\frac{1}{3}& -\frac{1}{2}& 1& -\frac{1}{6}& & & & & \\ \frac{1}{12}& -\frac{2}{3}& 0& \frac{2}{3}& -\frac{1}{12}& & & & \\ & \frac{1}{12}& -\frac{2}{3}& 0& \frac{2}{3}& -\frac{1}{12}& & & \\ & & \ddots & \ddots & \ddots & \ddots & \ddots & & \\ & & & \frac{1}{12}& -\frac{2}{3}& 0& \frac{2}{3}& -\frac{1}{12}& \\ & & & & \frac{1}{12}& -\frac{2}{3}& 0& \frac{2}{3}& -\frac{1}{12}\\ & & & & & \frac{1}{6}& -1& \frac{1}{2}& \frac{1}{3}\\ & & & & & -\frac{1}{3}& \frac{3}{2}& -3& \frac{11}{6}\end{array}\right)}_{(M_S+1)\times (M_S+1)}\hfill \end{array}$$

24

$$\begin{array}{c}\hfill D_{M_S}^2=\frac{1}{{\left(\mathrm{\Delta },S\right)}^2}{\left(\begin{array}{ccccccccc}2& -5& 4& -1& & & & & \\ 1& -2& 1& & & & & & \\ -\frac{1}{6}& \frac{5}{3}& -3& \frac{5}{3}& -\frac{1}{6}& & & & \\ & -\frac{1}{6}& \frac{5}{3}& -3& \frac{5}{3}& -\frac{1}{6}& & & \\ & & \ddots & \ddots & \ddots & \ddots & \ddots & & \\ & & & -\frac{1}{6}& \frac{5}{3}& -3& \frac{5}{3}& -\frac{1}{6}& \\ & & & & -\frac{1}{6}& \frac{5}{3}& -3& \frac{5}{3}& -\frac{1}{6}\\ & & & & & & 1& -2& 1\\ & & & & & -1& 4& -5& 2\end{array}\right)}_{(M_S+1)\times (M_S+1)}.\end{array}$$Now, we define uniform partitions ${\mathcal{V}}_{M_V}=\{V_{m_v}=V_{min}+m_v\mathrm{\Delta }V,m_v=0,1,\cdots ,M_V,\mathrm{\Delta }V=\frac{V_{max}-V_{min}}{M_V}\}$ and ${\mathcal{K}}_{M_K}=\{K_{m_k}=K_{min}+m_k\mathrm{\Delta }K,m_k=0,1,\cdots ,M_K,\mathrm{\Delta }K=\frac{K_{max}-K_{min}}{M_K}\}$ for intervals $[V_{min},V_{max}]$ and $[K_{min},K_{max}]$, respectively.

Using cubic spline quasi-interpolation (18) and the above method, the following approximations are obtained for derivatives in (16):

25

$$\begin{array}{cc}& (C_S^n(S_0,V_{m_v},K_{m_k}),C_S^n(S_1,V_{m_v},K_{m_k}),\dots ,C_S^n(S_{M_S},V_{m_v},K_{m_k}){)}^T\hfill \\ \hfill & =D_{M_S}^1(C^n(S_0,V_{m_v},K_{m_k}),C^n(S_1,V_{m_v},K_{m_k}),\dots ,C^n(S_{M_S},V_{m_v},K_{m_k}){)}^T,\hfill \\ \hfill & (C_{\mathit{SS}}^n(S_0,V_{m_v},K_{m_k}),C_{\mathit{SS}}^n(S_1,V_{m_v},K_{m_k}),\dots ,C_{\mathit{SS}}^n(S_{M_S},V_{m_v},K_{m_k}){)}^T\hfill \\ \hfill & =D_{M_S}^2(C^n(S_0,V_{m_v},K_{m_k}),C^n(S_1,V_{m_v},K_{m_k}),\dots ,C^n(S_{M_S},V_{m_v},K_{m_k}){)}^T,\hfill \\ \hfill & (C_V^n(S_{m_s},V_0,K_{m_k}),C_V^n(S_{m_s},V_1,K_{m_k}),\dots ,C_V^n(S_{m_s},V_{M_V},K_{m_k}){)}^T\hfill \\ \hfill & =D_{M_V}^1(C^n(S_{m_s},V_0,K_{m_k}),C^n(S_{m_s},V_1,K_{m_k}),\dots ,C^n(S_{m_s},V_{M_V},K_{m_k}){)}^T,\hfill \\ \hfill & (C_{\mathit{VV}}^n(S_{m_s},V_0,K_{m_k}),C_{\mathit{VV}}^n(S_{m_s},V_1,K_{m_k}),\dots ,C_{\mathit{VV}}^n(S_{m_s},V_{M_V},K_{m_k}){)}^T\hfill \\ \hfill & =D_{M_V}^2(C^n(S_{m_s},V_0,K_{m_k}),C^n(S_{m_s},V_1,K_{m_k}),\dots ,C^n(S_{m_s},V_{M_V},K_{m_k}){)}^T,\hfill \\ \hfill & (C_K^n(S_{m_s},V_{m_v},K_0),C_K^n(S_{m_s},V_{m_v},K_1),\dots ,C_K^n(S_{m_s},V_{m_v},K_{M_K}){)}^T\hfill \\ \hfill & =D_{M_K}^1(C^n(S_{m_s},V_{m_v},K_0),C^n(S_{m_s},V_{m_v},K_1),\dots ,C^n(S_{m_s},V_{m_v},K_{M_K}){)}^T,\hfill \end{array}$$for each $m_s=1,2,\dots ,M_S$, $m_v=0,1,\dots ,M_V$, and $m_k=0,1,\dots ,M_K$. Let $D_{M_S}^1=\left(d_0^{M_S},d_1^{M_S},\dots ,d_{M_S}^{M_S}\right)$. According to (25), we have

26

$$\begin{array}{cc}\hfill \left(\begin{array}{c}{C}_S^n(S_0,V,K_{m_k})\\ C_S^n(S_1,V,K_{m_k})\\ ⋮\\ C_S^n(S_{M_S},V,K_{m_k})\end{array}\right)& =C^n(S_0,V,K_{m_k})d_0^{M_S}+C^n(S_1,V,K_{m_k})d_1^{M_S}\hfill \\ \hfill & +\cdots +C^n(S_{M_S},V,K_{m_k})d_{M_S}^{M_S}.\hfill \end{array}$$By taking the derivative from both sides of (26) with respect to V, we have

27

$$\begin{array}{cc}\hfill \left(\begin{array}{c}{C}_{\mathit{SV}}^n(S_0,V,K_{m_k})\\ C_{\mathit{SV}}^n(S_1,V,K_{m_k})\\ ⋮\\ C_{\mathit{SV}}^n(S_{M_S},V,K_{m_k})\end{array}\right)& =C_V^n(S_0,V,K_{m_k})d_0^{M_S}+C_V^n(S_1,V,K_{m_k})d_1^{M_S}\hfill \\ \hfill & +\cdots +C_V^n(S_{M_S},V,K_{m_k})d_{M_S}^{M_S}.\hfill \end{array}$$From (27) and using (25), we can calculate $C_{\mathit{SV}}^n(S_{m_s},V_{m_v},K_{m_k})$ in each point. Note that $C^0(S_{m_s},V_{m_v},K_{m_k})$ can be evaluated from initial condition for $0\le m_s\le M_S$, $0\le m_v\le M_V$ and $0\le m_k\le M_K$. In the following, we present Algorithm 1 to describe the essential procedures for implementing the proposed methodology.

Error analysis

In this section, first, we examine the stability and then we discuss about the convergence of the presented numerical method.

Stability analysis

According to (21) and (22) and matrices (23) and (24), we obtain the partial derivatives of $C^n$ for $1<m_s<M_S$, $1<m_v<M_V$ and $1<m_k<M_K$ as

28

$$\begin{array}{c}\hfill {(C_S)}_{m_s,m_v,m_k}^n=\frac{1}{3\mathrm{\Delta }S}(\frac{1}{4}{C}_{m_s-2,m_v,m_k}^n-2C_{m_s-1,m_v,m_k}^n+2C_{m_s+1,m_v,m_k}^n-\frac{1}{4}{C}_{m_s+2,m_v,m_k}^n),\end{array}$$

29

$$\begin{array}{c}\hfill {(C_V)}_{m_s,m_v,m_k}^n=\frac{1}{3\mathrm{\Delta }V}(\frac{1}{4}{C}_{m_s,m_v-2,m_k}^n-2C_{m_s,m_v-1,m_k}^n+2C_{m_s,m_v+1,m_k}^n-\frac{1}{4}{C}_{m_s,m_v+2,m_k}^n),\end{array}$$

30

$$\begin{array}{cc}\hfill {(C_{\mathit{SS}})}_{m_s,m_v,m_k}^n& =\frac{1}{{(\mathrm{\Delta }S)}^2}(-\frac{1}{6}{C}_{m_s-2,m_v,m_k}^n+\frac{5}{3}{C}_{m_s-1,m_v,m_k}^n-3C_{m_s,m_v,m_k}^n\hfill \\ \hfill & +\frac{5}{3}{C}_{m_s+1,m_v,m_k}^n-\frac{1}{6}{C}_{m_s+2,m_v,m_k}^n),\hfill \end{array}$$

31

$$\begin{array}{cc}\hfill {(C_{\mathit{VV}})}_{m_s,m_v,m_k}^n& =\frac{1}{{(\mathrm{\Delta }V)}^2}(-\frac{1}{6}{C}_{m_s,m_v-2,m_k}^n+\frac{5}{3}{C}_{m_s,m_v-1,m_k}^n-3C_{m_s,m_v,m_k}^n\hfill \\ \hfill & +\frac{5}{3}{C}_{m_s,m_v+1,m_k}^n-\frac{1}{6}{C}_{m_s,m_v+2,m_k}^n),\hfill \end{array}$$

32

$$\begin{array}{cc}\hfill {(C_{\mathit{SV}})}_{m_s,m_v,m_k}^n& =\frac{1}{9\mathrm{\Delta }S\mathrm{\Delta }V}(\frac{1}{16}[C_{m_s-2,m_v-2,m_k}^n-C_{m_s-2,m_v+2,m_k}^n\hfill \\ \hfill & -C_{m_s+2,m_v-2,m_k}^n+C_{m_s+2,m_v+2,m_k}^n]\hfill \\ \hfill & +\frac{1}{2}[-C_{m_s-2,m_v-1,m_k}^n+C_{m_s-2,m_v+1,m_k}^n\hfill \\ \hfill & -C_{m_s-1,m_v-2,m_k}^n+C_{m_s-1,m_v+2,m_k}^n\hfill \\ \hfill & +C_{m_s+1,m_v-2,m_k}^n-C_{m_s+1,m_v+2,m_k}^n\hfill \\ \hfill & +C_{m_s+2,m_v-1,m_k}^n-C_{m_s+2,m_v+1,m_k}^n]\hfill \\ \hfill & +4[C_{m_s-1,m_v-1,m_k}^n-C_{m_s-1,m_v+1,m_k}^n\hfill \\ \hfill & -C_{m_s+1,m_v-1,m_k}^n+C_{m_s+1,m_v+1,m_k}^n]),\hfill \end{array}$$According to (16) and (28)–(32), we consider the following numerical scheme

33

$$\begin{array}{cc}\hfill C_{m_s,m_v,m_k}^{n+1}& ={\rho }_n^{m_s,m_v}{C}_{m_s,m_v,m_k}^n+S_s^{m_s,m_v}{C}_{m_s-2,m_v,m_k}^n+Q_s^{m_s,m_v}{C}_{m_s-1,m_v,m_k}^n\hfill \\ \hfill & +R_s^{m_s,m_v}{C}_{m_s+1,m_v,m_k}^n-W_s^{m_s,m_v}{C}_{m_s+2,m_v,m_k}^n+S_v^{n,m_v}{C}_{m_s,m_v-2,m_k}^n\hfill \\ \hfill & +Q_v^{n,m_v}{C}_{m_s,m_v-1,m_k}^n+R_v^{n,m_v}{C}_{m_s,m_v+1,m_k}^n-W_v^{n,m_v}{C}_{m_s,m_v+2,m_k}^n\hfill \\ \hfill & +\frac{{\lambda }_k^{m_s,m_k}}{12}{C}_{m_s,m_v,m_k-2}^n-\frac{2{\lambda }_k^{m_s,m_k}}{3}{C}_{m_s,m_v,m_k-1}^n+\frac{2{\lambda }_k^{m_s,m_k}}{3}{C}_{m_s,m_v,m_k+1}^n\hfill \\ \hfill & -\frac{{\lambda }_k^{m_s,m_k}}{12}{C}_{m_s,m_v,m_k+2}^n+\frac{{\lambda }_n^{m_s,m_v}}{144}[C_{m_s-2,m_v-2,m_k}^n-C_{m_s-2,m_v+2,m_k}^n\hfill \\ \hfill & -C_{m_s+2,m_v-2,m_k}^n+C_{m_s+2,m_v+2,m_k}^n]\hfill \\ \hfill & +\frac{{\lambda }_n^{m_s,m_v}}{18}[-C_{m_s-2,m_v-1,m_k}^n+C_{m_s-2,m_v+1,m_k}^n-C_{m_s-1,m_v-2,m_k}^n\hfill \\ \hfill & +C_{m_s-1,m_v+2,m_k}^n+C_{m_s+1,m_v-2,m_k}^n-C_{m_s+1,m_v+2,m_k}^n+C_{m_s+2,m_v-1,m_k}^n\hfill \\ \hfill & -C_{m_s+2,m_v+1,m_k}^n]+\frac{4{\lambda }_n^{m_s,m_v}}{9}[C_{m_s-1,m_v-1,m_k}^n-C_{m_s-1,m_v+1,m_k}^n\hfill \\ \hfill & -C_{m_s+1,m_v-1,m_k}^n+C_{m_s+1,m_v+1,m_k}^n],\hfill \end{array}$$where$$\begin{array}{cc}& {\rho }_n^{m_s,m_v}=1-r\mathrm{\Delta }\tau -\frac{3\mathrm{\Delta }\tau }{{(\mathrm{\Delta }S)}^2}{S}_{m_s}^{2}g(t_n)V_{m_v}-\frac{3\mathrm{\Delta }\tau }{{(\mathrm{\Delta }S)}^2}{\sigma }^{2}g(t_n)V_{m_v}\hfill \\ \hfill & {\lambda }_s^{m_s}=\frac{\mathrm{\Delta }\tau }{\mathrm{\Delta }S}{S}_{m_s},{\lambda }_v^{m_v}=\frac{\mathrm{\Delta }\tau }{\mathrm{\Delta }V}{V}_{m_v},{\lambda }_k^{m_s,m_k}=\frac{\mathrm{\Delta }\tau }{\mathrm{\Delta }K}f(S_{m_s},K_{m_k},t_n),\hfill \\ \hfill & {\lambda }_n^{m_s,m_v}=\frac{h(S_{m_s},V_{m_v},t_n)\mathrm{\Delta }\tau }{\mathrm{\Delta }S\mathrm{\Delta }V},{\delta }_n^{m_s,m_v}=\frac{S_{m_s}{V}_{m_v}g(t_n)}{\mathrm{\Delta }S},{\delta }_n=\frac{{\sigma }^{2}g(t_n)}{\mathrm{\Delta }V},\hfill \\ \hfill & S_s^{m_s,m_v}=\frac{{\lambda }_s^{m_s}}{6}\left(\frac{r}{2}-{\delta }_n^{m_s,m_v}\right),Q_s^{m_s,m_v}=\frac{{\lambda }_s^{m_s}}{3}\left(-2r+5{\delta }_n^{m_s,m_v}\right),\hfill \\ \hfill & R_s^{m_s,m_v}=\frac{{\lambda }_s^{m_s}}{3}\left(2r+5{\delta }_n^{m_s,m_v}\right),W_s^{m_s,m_v}=\frac{{\lambda }_s^{m_s}}{6}\left(\frac{r}{2}+{\delta }_n^{m_s,m_v}\right),\hfill \\ \hfill & S_v^{n,m_v}=\frac{{\lambda }_v^{m_v}}{6}\left(\frac{r}{2}-{\delta }_n\right),Q_v^{n,m_v}=\frac{{\lambda }_v^{m_v}}{3}\left(-2r+5{\delta }_n\right),\hfill \\ \hfill & R_v^{n,m_v}=\frac{{\lambda }_v^{m_v}}{3}\left(2r+5{\delta }_n\right),W_v^{n,m_v}=\frac{{\lambda }_v^{m_v}}{6}\left(\frac{r}{2}+{\delta }_n\right).\hfill \end{array}$$Let $\mathrm{\Delta }\mathrm{S}\le \mathrm{\Delta }\mathrm{V}$ and $|C_{s^{\ast },v^{\ast },k^{\ast }}^n|=\underset{m_s,m_v,m_k}{sup}|C_{m_s,m_v,m_k}^n|$. By taking the absolute value and supremum from both sides of (33), we have

34

$$\begin{array}{c}\hfill \underset{m_s,m_v,m_k}{sup}\left|C_{m_s,m_v,m_k}^{n+1}\right|\le P_{m_s,m_v,m_k}^n\left|C_{s^{\ast },v^{\ast },k^{\ast }}^n\right|,\end{array}$$where$$\begin{array}{cc}\hfill P_{m_s,m_v,m_k}^n& =\left|{\rho }_n^{m_s,m_v}\right|+\left|S_s^{m_s,m_v}\right|+\left|Q_s^{m_s,m_v}\right|+\left|R_s^{m_s,m_v}\right|+\left|W_s^{m_s,m_v}\right|\hfill \\ \hfill & +\left|S_v^{n,m_v}\right|+\left|Q_v^{n,m_v}\right|+\left|R_v^{n,m_v}\right|+\left|W_v^{n,m_v}\right|+\frac{3}{2}\left|{\lambda }_k^{m_s,m_k}\right|+\frac{9}{4}{\lambda }_n^{m_s,m_v}.\hfill \end{array}$$Therefore, if $P_{m_s,m_v,m_k}^n\le 1$, the presented method demonstrates stability. In the case of alternative values for $m_s$, $m_v$, and $m_k$, this result can be easily derived with minor adjustments.

Convergence analysis

In (33), by substituting the Taylor expansion for all sentences $C_{m_s+i,m_v,m_k}^n$, $C_{m_s,m_v+i,m_k}^n$, $C_{m_s,m_v,m_k+i}^n$ and $C_{m_s+i,m_v+j,m_k}^n$, $i,j=-2,-1,1,$ and 2, we have

35

$$\begin{array}{cc}\hfill C_{m_s,m_v,m_k}^{n+1}& =(1-r\mathrm{\Delta }\tau +3\mathrm{\Delta }\tau {\sigma }^{2}g(t_n)V_{m_v}\left(\frac{1}{{(\mathrm{\Delta }V)}^2}-\frac{1}{{(\mathrm{\Delta }S)}^2}\right))C_{m_s,m_v,m_k}^n\hfill \\ \hfill & +\mathrm{\Delta }\tau (S_{m_s}r{(C_S)}_{m_s,m_v,m_k}^n+V_{m_v}r{(C_V)}_{m_s,m_v,m_k}^n\hfill \\ \hfill & +f(S_{m_s},K_{m_k},t_n){(C_K)}_{m_s,m_v,m_k}^n\hfill \\ \hfill & +g(t_n)V_{m_v}\left(S_{m_s}^2{(C_{\mathit{SS}})}_{m_s,m_v,m_k}^n+{\sigma }^2{(C_{\mathit{VV}})}_{m_s,m_v,m_k}^n\right)\hfill \\ \hfill & +h(S_{m_s},V_{m_v},t_n){(C_{\mathit{SV}})}_{m_s,m_v,m_k}^n)\hfill \\ \hfill & +O\left({(\mathrm{\Delta }S)}^3+{(\mathrm{\Delta }S)}^2(\mathrm{\Delta }V)+{(\mathrm{\Delta }V)}^2(\mathrm{\Delta }S)+{(\mathrm{\Delta }V)}^3+{(\mathrm{\Delta }K)}^3\right).\hfill \end{array}$$By rewriting (35), we get

36

$$\begin{array}{cc}& TE:=\frac{C_{m_s,m_v,m_k}^{n+1}-C_{m_s,m_v,m_k}^n}{\mathrm{\Delta }\tau }-3{\sigma }^{2}g(t_n)V_{m_v}\left(\frac{1}{{(\mathrm{\Delta }V)}^2}-\frac{1}{{(\mathrm{\Delta }S)}^2}\right)C_{m_s,m_v,m_k}^n\hfill \\ \hfill & -r{S}_{m_s}{(C_S)}_{m_s,m_v,m_k}^n-r{V}_{m_v}{(C_V)}_{m_s,m_v,m_k}^n-f(S_{m_s},K_{m_k},t_n){(C_K)}_{m_s,m_v,m_k}^n\hfill \\ \hfill & -g(t_n)V_{m_v}\left(S_{m_s}^2{(C_{\mathit{SS}})}_{m_s,m_v,m_k}^n+{\sigma }^2{(C_{\mathit{VV}})}_{m_s,m_v,m_k}^n\right)\hfill \\ \hfill & -h(S_{m_s},V_{m_v},t_n){(C_{\mathit{SV}})}_{m_s,m_v,m_k}^n+r{C}_{m_s,m_v,m_k}^n\hfill \\ \hfill & +O\left({(\mathrm{\Delta }S)}^3+{(\mathrm{\Delta }S)}^2(\mathrm{\Delta }V)+{(\mathrm{\Delta }V)}^2(\mathrm{\Delta }S)+{(\mathrm{\Delta }V)}^3+{(\mathrm{\Delta }K)}^3\right).\hfill \end{array}$$From (15), we obtain

37

$$\begin{array}{cc}\hfill {(C_{\tau })}_{m_s,m_v,m_k}^n& =r\left(S_{m_s}{(C_S)}_{m_s,m_v,m_k}^n+V_{m_v}{(C_V)}_{m_s,m_v,m_k}^n\right)\hfill \\ \hfill & +f(S_{m_s},K_{m_k},t_n){(C_K)}_{m_s,m_v,m_k}^n\hfill \\ \hfill & +g(t_n)V_{m_v}(S_{m_s}^2{(C_{\mathit{SS}})}_{m_s,m_v,m_k}^n+{\sigma }^2{(C_{\mathit{VV}})}_{m_s,m_v,m_k}^n)\hfill \\ \hfill & +h(S_{m_s},V_{m_v},t_n){(C_{\mathit{SV}})}_{m_s,m_v,m_k}^n-r{C}_{m_s,m_v,m_k}^n.\hfill \end{array}$$By substituting (37) in (36), we get

38

$$\begin{array}{cc}\hfill \mathit{TE}& =\frac{C_{m_s,m_v,m_k}^{n+1}-C_{m_s,m_v,m_k}^n}{\mathrm{\Delta }\tau }-{(C_{\tau })}_{m_s,m_v,m_k}^n\hfill \\ \hfill & -3{\sigma }^{2}g(t_n)V_{m_v}\left(\frac{1}{{(\mathrm{\Delta }V)}^2}-\frac{1}{{(\mathrm{\Delta }S)}^2}\right)C_{m_s,m_v,m_k}^n\hfill \\ \hfill & +O\left({(\mathrm{\Delta }S)}^3+{(\mathrm{\Delta }S)}^2(\mathrm{\Delta }V)+{(\mathrm{\Delta }V)}^2(\mathrm{\Delta }S)+{(\mathrm{\Delta }V)}^3+{(\mathrm{\Delta }K)}^3\right).\hfill \end{array}$$Therefore, if $\mathrm{\Delta }V=\mathrm{\Delta }S=\mathrm{\Delta }K$, then $TE=O(\mathrm{\Delta }\tau +{(\mathrm{\Delta }S)}^3)$ and the method is convergent.

Fig. 1 [Images not available. See PDF.]

The values of C(T, S, 0.3, K) for some K values

Fig. 2 [Images not available. See PDF.]

The values of C(T, S, V, 3) for $n^{\ast }=1$ (a), $n^{\ast }=2$ (b), and $n^{\ast }=3$ (c)

Fig. 3 [Images not available. See PDF.]

The values of C(T, S, V, K) for $K=3$, $K=5$, $K=7$, $K=9$ with $n^{\ast }=1$ (a) and for $K=3$, $K=9$ with $n^{\ast }=2$ (b)

Fig. 4 [Images not available. See PDF.]

The values of C(T, S, 0.3, K) for some K values

Numerical experiments

In this section, we obtain the numerical solution to PDE (14) using the presented method in MATLAB programming software. To estimate the solution of the (14), first, we consider the following assumptions: $r=0.1$, $H=0.8$, $\rho =0.6$, $\sigma =0.1$, $T=1$, $[S_{min},S_{max}]=[0.1,5.1]$, $[V_{min},V_{max}]=[0.1,0.4]$, $[K_{min},K_{max}]=[1,6]$, $M_S=M_K=25$, $M_V=6$, $\mathrm{\Delta }t={10}^{-5}$ and $m_1=m_2=0.1$. In Fig. 1, the values of C in $T=1$, $V=0.3$ and different values of K with $n^{\ast }=1$ are demonstrated. The values of C in $\tau =T$ and $K=3$ and different values of $n^{\ast }$ are illustrated in Fig. 2. Also, we consider $r=0.1$, $H=0.9$, $\rho =0.4$, $\sigma =0.01$, $T=1$, $S_{min}=0.5$, $S_{max}=10.5$, $V_{min}=0.1$, $V_{max}=5.1$, $K_{min}=1$, $K_{max}=11$, $MS=MK=50$, $MV=25$ and $\mathrm{\Delta }t={10}^{-7}$. In Fig. 3, the values of C(T, S, V, K) for different values of K and $n^{\ast }$ are demonstrated. Finally, the following assumptions are considered $r=0.1$, $H=0.9$, $\rho =0.4$, $\sigma =0.01$, $T=1$, $S_{min}=0.5$, $S_{max}=20.5$, $V_{min}=0.1$, $V_{max}=6.1$, $K_{min}=1$, $K_{max}=21$, $n^{\ast }=1$, $MS=MK=60$, $MV=18$ and $\mathrm{\Delta }t={10}^{-7}$. The values of C(T, S, V, K) for different values of K and $n^{\ast }=1$ are shown in Fig. 4.

Conclusions

In this paper, we examined the pricing of the geometric Asian power option under the mixed stochastic volatility Heston model. The model possessed properties of long memory and stochastic volatility, allowing for the prediction of stock prices and various financial assets. However, calculating the price of the geometric Asian power option under this model proves to be a complicated endeavor. The option price partial differential equation (PDE) in this model takes on a multidimensional form, with no analytical solution available. Therefore, we presented an efficient method to approximate the solution. To achieve this, we proposed a discrete-time scheme using the Euler method and employed the cubic B-Spline pseudo-interpolation approximation to approximate the other derivatives in the equation. Also, we obtained the numerical solution of the equation by applying the final scheme, along with the initial and boundary conditions. Furthermore, we discussed the stability of the method and presented numerical results to support our findings. The main advantage of the proposed method lies in its simplicity for computer implementation. Additionally, the method can be easily applied to multi-dimensional problems with a greater number of independent variables.

Author contribution

Abdulaziz Alsenafi and Fares Alazemi wrote the main manuscript text. Also, they have completed the numerical experiments. Javad Alavi provided assistance in proving the theorems in this manuscript. All authors reviewed the manuscript.

Data availibility

The datasets generated during the current study are available.

Declarations

Ethical approval

The manuscript is not submitted to more than one journal for simultaneous consideration. The manuscript is original and is not published elsewhere in any form or language. The manuscript is not divided into several parts to increase the quantity of submissions and submitted to various journals or to one journal over time. Results are presented clearly, honestly, and without fabrication, falsification, or inappropriate data manipulation (including image-based manipulation). No data, text, or theories by others are not presented without references.

Competing interests

The authors declare no competing interests.

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

References

1. Peng, B; Peng, F. Pricing Asian power options under jump-fraction process. J. Econ. Finance Adm. Sci.; 2012; 17, 33 pp. 2-9.1265.91062

2. Kim, J; Kim, B; Moon, KS; Wee, IS. Valuation of power options under Heston’s stochastic volatility model. J. Econ. Dyn. Control.; 2012; 36, 11 pp. 1796-1813.2965068[DOI: https://dx.doi.org/10.1016/j.jedc.2012.05.005] 1345.91073

3. Lee, Y; Kim, Y; Lee, J. Pricing various types of power options under stochastic volatility. Symmetry.; 2020; 12, 11 1911. [DOI: https://dx.doi.org/10.3390/sym12111911] 07551093

4. Rao, BP. Pricing geometric Asian power options under mixed fractional Brownian motion environment. Phys. A: Stat. Mech.; 2016; 446, pp. 92-99.3437666[DOI: https://dx.doi.org/10.1016/j.physa.2015.11.013] 1400.91607

5. Edeki, S.O., Motsepa, T., Khalique, C.M., Akinlabi, G.O.: The Greek parameters of a continuous arithmetic Asian option pricing model via Laplace Adomian decomposition method. Open Phys. 16(1), 780–785 (2018)

6. Zhang, WG; Li, Z; Liu, YJ. Analytical pricing of geometric Asian power options on an underlying driven by a mixed fractional Brownian motion. Phys. A: Stat. Mech.; 2018; 490, pp. 402-418.3716370[DOI: https://dx.doi.org/10.1016/j.physa.2017.08.070] 1493.91130

7. Shokrollahi, F. The evaluation of geometric Asian power options under time changed mixed fractional Brownian motion. J. Comput. Appl. Math.; 2018; 344, pp. 716-724.3825546[DOI: https://dx.doi.org/10.1016/j.cam.2018.05.042] 1395.91467

8. Sunday O. Edeki, Member, Iaeng, Olabisi O. Ugbebor, Paul O. Ogundile.: Analytical solutions of a sontinuous arithmetic Asian model for option pricing using projected differential transform method. Eng. Lett. 27, 2 (2019)

9. Edeki, S.O., Akinlabi, G. O., González-Gaxiola, O.: Adomian decomposition method for analytical solution of a continuous arithmetic Asian option pricing model. Telkomnika. 17(2), 866–872 (2019)

10. Wang, W; Cai, G; Tao, X. Pricing geometric Asian power options in the sub-fractional Brownian motion environment. Chaos Solit. Fractals.; 2021; 145, 10754.4217794[DOI: https://dx.doi.org/10.1016/j.chaos.2021.110754] 1494.91164

11. Malhotra, G; Srivastava, R; Taneja, HC. Pricing of the geometric Asian options under a multifactor stochastic volatility model. J. Comput. Appl. Math.; 2022; 406, 4358134[DOI: https://dx.doi.org/10.1016/j.cam.2021.113986] 1490.91219113986.

12. Ahmadian, D., Ballestra, L.V., Shokrollahi, F.: A Monte-Carlo approach for pricing arithmetic Asian rainbow options under the mixed fractional Brownian motion. Chaos Solit. Fractals. 158, 112023 (2022)

13. Ewald, C.O., Wu, Y., Zhang, A.: Pricing Asian options with stochastic convenience yield and jumps. Quant. Finance. 23(4), 677–692 (2023)

14. Kaur, J., Natesan, S.: A novel numerical scheme for time-fractional Black-Scholes PDE governing European options in mathematical finance. Numer. Algorithms. 1–31 (2023)

15. An, X., Wang, Q., Liu, F., Anh, V.V., Turner, I.W.: Parameter estimation for time-fractional Black-Scholes equation with S &P 500 index option. Numer. Algorithms. 1–30 (2023)

16. Taghipour, M., Aminikhah, H.: Numerical solution of distributed-order fractional Korteweg-de Vries equation via fractional Zigzag rising diagonal functions. Numer. Algorithms. 1–38 (2023)

17. Terekhov, AV. A three-dimensional Laguerre one-way wave equation solver. Appl Numer Math.; 2022; 173, pp. 380-394.4364763[DOI: https://dx.doi.org/10.1016/j.apnum.2022.01.001] 1486.65131

18. Dehghan, M; Abbaszadeh, M. An efficient technique based on finite difference/finite element method for solution of two-dimensional space/multi-time fractional Bloch-Torrey equations. Appl Numer Math.; 2018; 131, pp. 190-206.3807178[DOI: https://dx.doi.org/10.1016/j.apnum.2018.04.009] 1395.65074

19. Safaei, M; Neisy, A; Nematollahi, N. New splitting scheme for pricing American options under the Heston model. Comput. Econ.; 2018; 52, pp. 405-420. [DOI: https://dx.doi.org/10.1007/s10614-017-9686-4] 1499.91175

20. Haentjens, T. Efficient and stable numerical solution of the Heston-Cox-Ingersoll-Ross partial differential equation by alternating direction implicit finite difference schemes. Int. J. Comput. Math.; 2023; 90, 11 pp. 2409-2430.3171907[DOI: https://dx.doi.org/10.1080/00207160.2013.777710] 1281.91184

21. Safaei, M; Neisy, A; Nematollahi, N. Generalized componentwise splitting scheme for option pricing under the Heston-Cox-Ingersoll-Ross model. Journal of Statistical Theory and Applications.; 2019; 18, 4 pp. 425-438. [DOI: https://dx.doi.org/10.2991/jsta.d.191209.001] 1463.91201

22. Jebreen, HB. A Gaussian radial basis function-finite difference technique to simulate the HCIR equation. J. Comput. Appl. Math.; 2019; 347, pp. 181-195.3867496[DOI: https://dx.doi.org/10.1016/j.cam.2018.08.019] 1405.65132

23. Biazar, J; Goldoust, F. Multi-dimensional Legendre wavelets approach on the Black-Scholes and Heston Cox Ingersoll Ross equations. AIMS math.; 2019; 4, 4 pp. 1046-1064.4137019[DOI: https://dx.doi.org/10.3934/math.2019.4.1046] 1484.91514

24. Longjin, L; Ren, FY; Qiu, WY. The application of fractional derivatives in stochastic models driven by fractional Brownian motion. Phys. A: Stat. Mech. Appl.; 2010; 389, 21 pp. 4809-4818.2684667[DOI: https://dx.doi.org/10.1016/j.physa.2010.06.016] 1260.82064

25. Alavi, J; Aminikhah, H. Numerical study of the inverse problem of generalized Burgers-Fisher and generalized Burgers-Huxley equations. Adv. Math. Phys.; 2021; 2021, pp. 1-15.4235842[DOI: https://dx.doi.org/10.1155/2021/6652108] 1481.65176

26. Aminikhah, H., Alavi, J.: Applying cubic B-Spline quasi-interpolation to solve 1D wave equations in polar coordinates. nt. sch. res. notices. (2013)

27. Sablonniere, P.: Univariate spline quasi-interpolants and applications to numerical analysis. math/0504022. (2005)