A New Parallel Difference Method for Solving Time Fractional Black-Scholes Model

Xuebin Yang¹, Lifei Wu^2*, Yu Zhang²

¹China Academy of Public Finance and Public Policy, Central University of Finance and Economics, Beijing, China.
²School of Mathematics and Physics, North China Electric Power University, Beijing, China.

DOI: 10.4236/jmf.2022.124036                

Abstract

A modified Black-Scholes (B-S) model with time fractional derivative is studied when the price change of the underlying is considered as a fractal transmission system. It is very practical in application to study the numerical computation of this time fractional B-S model (TFBSM) governing European options. This paper constructs mixed alternative segment Crank-Nicolson (MASC-N) parallel difference scheme, which accuracy is spatially 2 order and temporally 2-α order. In addition, we have provided a theoretical proof for the stability and convergence of the MASC-N scheme, and demonstrated the accuracy and effectiveness of the proposed method through numerical experiments. The numerical results show that the MASC-N scheme is easy to implement when applied to time fractional B-S model.

Keywords

Time Fractional B-S Model, MASC-N Scheme, Stability, Convergence, Numerical Experiments

1. Introduction

Option is one of the most widely used financial derivatives, and it has shown very important significance in both theory and practice to study the option pricing. The Black-Scholes (B-S) model, which was proposed by Black F., Scholes M. and Merton R. (1973) [1] [2], led to a boom in options trading and scientifically legitimized the activities of the Chicago board options exchange and other option markets around the world. The classical B-S model is of seven basic assumptions, sometimes not meet the actual financial market applications. Therefore, the researchers improved the classical B-S model and obtained some new models such as B-S model with transaction cost [3] [4], jump-diffusion model [5], stochastic volatility model [6] and fractional B-S model [7] [8].

Through observation and research on the stock market, it is found that the most fundamental characteristics and basic state of the capital market are random fluctuations, and some complex dynamic systems are usually fractional, such as fractional Brownian motion, continuous time walk, etc. Based on the fractional stochastic differential equation driven by fractional Brownian motion, many scholars have made some progress in the research of fractional B-S equation in recent years. W. Wyss (2000) [9] deduced a time fractional B-S equation governing European call option. Cartea A. and Del-Castillo-Negrete D. (2007) [10] obtained several fractional diffusion models for pricing option in market with jumps. By considering the fractional dynamics driven, Jumarie G. (2008, 2010) [11] [12] derived time fractional B-S model and time-and-space fractional B-S model. Soon after that, Zhang J.R. et al. (2010) [13] assumed that the underlying price change follows a fractional Itô process and the change in option price with time is a fractal transmission system, obtained a bi-fractional Black-Scholes-Merton model. Chen W. et al. (2015) [14] proposed a simplified version of Liang et al.’s model, they assumed that the underlying price change still follows the classical Brownian motion, but considered fractal transmission system. As a result, they obtained double barrier options pricing based on time fractional B-S model.

The analytical solution of the fractional B-S models is usually via integral transform methods [9] [12] [13] [14], homotopy perturbation methods [15], wavelet based hybrid methods [16], or via Lie symmetry transformations [17] and so on. The solutions usually are difficult to calculate, because they usually contain the form of convolution of some special functions or infinite series with an integral. Therefore, studying numerical approximate solutions of fractional B-S models is a very practical and important research project. The existing numerical algorithms for solving the fractional order differential equations mainly include: finite difference method, finite element method, spectral method, moving mesh method, series approximation method, matrix transformation method and so on [18] [19] [20]. Roughly speaking, finite difference method has been well established.

For solving the time fractional B-S equation, a few achievements on the numerical methods are as follow. Song L.N. and Wang W.G. (2013) [21] employed implicit finite difference method to solve the time fractional B-S equation together with the conditions satisfied by the standard put option. Kumar S. and Singh D. (2014) [22] presented a numerical algorithm for the time fractional B-S equation with boundary condition by homotopy perturbation method and homotopy analysis method. Bhowmik S.K. (2014) [23] proposed an explicit-implicit numerical scheme with a low order of convergence for a partial integro-differential equation that arises in option pricing theory by using a finite difference technique. Yang X.Z. et al. (2015) [24] deduced an Explicit-Implicit and Implicit-Explicit difference scheme for time fractional B-S equation. Zhang H. et al. (2016) [25] derived an implicit difference scheme for the time fractional B-S equation governing European options. Phaochoo P. et al. (2016) [26] proposed a meshless local Petrov-Galerkin method, which involves not only a meshless interpolation for the trial functions, but also a meshless integration of the weak-form. Nuugulu S.M. et al. (2021) [27] propose a corresponding robust numerical method which is based on the extension of a Crank-Nicolson (C-N) finite difference method. For solving time-fractional B-S equation, An X. et al. (2021) [28] proposed space-time spectral method employs the Jacobi polynomials for the temporal discretisation and Fourier-like basis functions for the spatial discretisation. Based on the alternating segment technology, Yan R.F. et al. (2021) [29] applied C-N format, and four kinds of Saul’yev asymmetric format to construct the alternating segmented C-N parallel scheme. Sarboland M. and Aminataei A. (2022) [30] applied the multiquadric quasi-interpolation scheme and the integrated radial basis function networks scheme to provide a numerical method to approximate the solution of the time fractional Black-Sholes equation.

With the rapid development of multi-core and cluster technology, it is of great theoretical significance and practical value to construct the parallel difference scheme with good stability for solving time fractional B-S model (TFBSM). In combination with the call option, we construct the mixed alternative segment C-N (MASC-N) scheme for TFBSM. The stability and convergence of the MASC-N scheme are analyzed. Finally, numerical examples demonstrate the effectiveness and accuracy of the MASC-N scheme for solving the TFBSM.

2. MASC-N Parallel Difference Method

2.1. Time Fractional Black-Scholes Model

In the present work, the option price $P\left(S,t\right)$ is suggested to be subject to the time fractional Black-Scholes equation with the following form [25] [31]:

$\frac{{\partial }^{\alpha }P\left(S,t\right)}{\partial t^{\alpha }}+\frac{1}{2}{\sigma }^2{S}^2\frac{{\partial }^{2}P\left(S,t\right)}{\partial S^2}+rS\frac{\partial P\left(S,t\right)}{\partial S}-rP\left(S,t\right)=0.$ (1)

where $t>0$, $0<\alpha \le 1$, P is the price of European call option; S and t are asset price and time; r is risk free interest rate; $\sigma$ represents volatility of underlying asset. The fractional derivative $P_t^{\left(\alpha \right)}\left(S,t\right)$ is Riemann-Liouville time fractional derivative with respect to t; $P_S$ and $P_{SS}$ refer to first derivative and second derivative with respect to S.

The value of European call option is taken as a solution of (1) with the following terminal and boundary conditions [1]:

1) $P\left(S,T\right)=\max \left\{S-K,0\right\}$. The profit and loss when the option expires is its price. Here K is the exercise price;

2) $S\to \infty$, $P\left(S,t\right)~S-K{\text{e}}^{-r\left(T-t\right)}$. When S is sufficiently large, the option price is close to $S-K{\text{e}}^{-r\left(T-t\right)}$, T is the due date of the options;

3) $P\left(0,t\right)=0$, which implies that the option price is approaching to zero when S is zero.

The European call option pricing model is to solve the follow time fractional B-S equation:

$\{\begin{array}{l}\frac{{\partial }^{\alpha }P\left(S,t\right)}{\partial t^{\alpha }}+\frac{1}{2}{\sigma }^2{S}^2\frac{{\partial }^{2}P\left(S,t\right)}{\partial S^2}+rS\frac{\partial P\left(S,t\right)}{\partial S}-rP\left(S,t\right)=0,\\ P\left(S,T\right)=\max \left\{S-K,0\right\}.\end{array}$ (2)

The boundary conditions $P\left(0,t\right)=0$, $\underset{S\to \infty }{\lim }P\left(S,t\right)=S-K{\text{e}}^{-r\left(T-t\right)}$, and the solution region is $\Sigma =\left\{0\le S\le \infty ,0\le t\le T\right\}$.

Equation (2) is an anti-variable coefficient parabolic equation. With the change of variables [2] [25]: $S={\text{e}}^x,t=T-\tau ,V\left(x,\tau \right)=P\left({\text{e}}^x,T-\tau \right)$.

Equation (2) can be transformed into the parabolic equations as follows:

$\{\begin{array}{l}{V}_{\tau }^{\left(\alpha \right)}\left(x,\tau \right)+\left(r-\frac{1}{2}{\sigma }^2\right)V_x\left(x,\tau \right)-\frac{1}{2}{\sigma }^2{V}_{xx}\left(x,\tau \right)-rV\left(x,\tau \right)=0,\\ V\left(x,0\right)=\max \left\{{\text{e}}^x-K,0\right\}.\end{array}$ (3)

The solution region converses into: ${\Sigma }_0=\left\{-\infty \le x\le +\infty ,0\le \tau \le T\right\}$.

To solve the TFBSM numerically it is necessary to truncate the original unbounded domain into a finite interval. Here we truncate the range of variable x in problem (1) to a finite interval $\left[M^{-},M^{+}\right]$. Then the solution region we consider is of the following finite domain: ${\Sigma }_1=\left\{M^{-}\le x\le M^{+},0\le \tau \le T\right\}$.

Meanwhile, boundary conditions are transformed into the form:

$V\left(M^{+},\tau \right)={\text{e}}^{M^{+}}-K{\text{e}}^{-r\tau },V\left(M^{-},\tau \right)=0$.

2.2. Construction of the MASC-N Scheme

Let $h=\left(M^{+}-M^{-}\right)/M$ and $k=T/N$, $\left(M,N\in Z_{+}\right)$ be space and time steps. The computation domain ${\Sigma }_1$ is discretized by uniform grid $\left(x_i,{\tau }_n\right)$ :

$\{\begin{array}{l}{x}_i=M^{-}+\left(i-1\right)h,i=1,2,\cdots ,M+1,\\ {\tau }_n=\left(n-1\right)k,n=1,2,\cdots ,N+1.\end{array}$

$V_i^n=V\left(x_i,{\tau }_n\right)$ denotes an approximate solution (exact solution of MASC-N scheme) of (3) in $x_i$ and ${\tau }_n$. The corresponding initial and boundary conditions are:

$V_i^1=f\left(x_i\right)=\max \left\{{\text{e}}^{x_i}-K,0\right\}$, $V_{M+1}^n=f_2\left({\tau }_n\right)={\text{e}}^{M^{+}}-K{\text{e}}^{-r\left(n-1\right)k}$.

The discrete format of time fractional derivative $P_t^{\left(\alpha \right)}\left(S,t\right)$ in the point $\left(x_i,{\tau }_{n+1}\right)$ using L1 interpolation approximation of the modified Riemann-Liouville fractional derivative is as follows [32] [33]:

$\begin{array}{l}\frac{{\partial }^{\alpha }V\left(x_i,{\tau }_{n+1}\right)}{\partial {\tau }^{\alpha }}\\ =\frac{k^{-\alpha }}{\Gamma \left(2-\alpha \right)}{\displaystyle \underset{j=1}{\overset{n}{\sum }}\left[V\left(x_i,{\tau }_{n+2-j}\right)-V\left(x_i,{\tau }_{n+1-j}\right)\right]\left[j^{1-\alpha }-{\left(j-1\right)}^{1-\alpha }\right]}+O\left(k^{2-\alpha }\right).\end{array}$

Letting $l_j=j^{1-\alpha }-{\left(j-1\right)}^{1-\alpha }$, we can define the operator $L_{\tau }^{\alpha }V\left(x_i,{\tau }_{n+1}\right)$ as:

$L_{\tau }^{\alpha }V\left(x_i,{\tau }_{n+1}\right)=\frac{k^{-\alpha }}{\Gamma \left(2-\alpha \right)}{\displaystyle \underset{j=1}{\overset{n}{\sum }}{l}_j\left[V\left(x_i,{\tau }_{n+2-j}\right)-V\left(x_i,{\tau }_{n+1-j}\right)\right]}.$ (4)

Moreover, the discrete formats of the space derivatives are:

$\begin{array}{l}\frac{\partial V\left(x_i,{\tau }_n\right)}{\partial x}=\frac{V_{i+1}^n-V_{i-1}^n}{2h}+O\left(h^2\right),\\ \frac{{\partial }^{2}V\left(x_i,{\tau }_n\right)}{\partial x^2}=\frac{V_{i+1}^n-2V_i^n+V_{i-1}^n}{h^2}+O\left(h^2\right).\end{array}$ (5)

Substituting (4) and (5) into (3), we can derive the classical explicit scheme, implicit scheme and C-N scheme of time fractional B-S Equation (3):

The classical explicit scheme of time fractional B-S Equation (3) is

$\begin{array}{l}{L}_{\tau }^{\alpha }V\left(x_i,{\tau }_{n+1}\right)\\ =\frac{1}{2}{\sigma }^2\frac{V_{i+1}^n-2V_i^n+V_{i-1}^n}{h^2}+\left(r-\frac{1}{2}{\sigma }^2\right)\frac{V_{i+1}^n-V_{i-1}^n}{2h}-r{V}_i^n+O\left(k^{2-\alpha }+h^2\right).\end{array}$

After combining similar terms and omitting the truncation errors, it can be rewritten as

$V_i^{n+1}=a{V}_{i-1}^n-2b{V}_i^n+c{V}_{i+1}^n+{\displaystyle \underset{j=1}{\overset{n-1}{\sum }}{w}_j{V}_i^{n+1-j}}+l_n{V}_i^1.$ (6)

The classical implicit scheme of time fractional B-S Equation (3) gives

$\begin{array}{c}{L}_{\tau }^{\alpha }V\left(x_i,{\tau }_{n+1}\right)=\frac{1}{2}{\sigma }^2\frac{V_{i+1}^{n+1}-2V_i^{n+1}+V_{i-1}^{n+1}}{h^2}+\left(r-\frac{1}{2}{\sigma }^2\right)\frac{V_{i+1}^{n+1}-V_{i-1}^{n+1}}{2h}\\ -r{V}_i^{n+1}+O\left(k^{2-\alpha }+h^2\right).\end{array}$

It can be rewritten as

$-a{V}_{i-1}^{n+1}+\left(1+2b\right)V_i^{n+1}-c{V}_{i+1}^{n+1}={\displaystyle \underset{j=1}{\overset{n-1}{\sum }}{w}_j{V}_i^{n+1-j}}+l_n{V}_i^1.$ (7)

The classical C-N scheme of time fractional B-S Equation (3) gives

$\begin{array}{l}-\frac{1}{2}a{V}_{i-1}^{n+1}+\left(1+b\right)V_i^{n+1}-\frac{1}{2}c{V}_{i+1}^{n+1}\\ =\frac{1}{2}a{V}_{i-1}^n-b{V}_i^n+\frac{1}{2}c{V}_{i+1}^n+{\displaystyle \underset{j=1}{\overset{n-1}{\sum }}{w}_j{V}_i^{n+1-j}}+l_n{V}_i^1.\end{array}$ (8)

where

$\begin{array}{l}a=\frac{m{\sigma }^2}{2h^2}-\frac{m\left(r-0.5{\sigma }^2\right)}{2h},b=\frac{m{\sigma }^2}{2h^2}+\frac{mr}{2},c=\frac{m{\sigma }^2}{2h^2}+\frac{m\left(r-0.5{\sigma }^2\right)}{2h},\\ m=k^{\alpha }\Gamma \left(2-\alpha \right),w_j=2j^{1-\alpha }-{\left(j+1\right)}^{1-\alpha }-{\left(j-1\right)}^{1-\alpha }.\end{array}$

Simulating the method of alternating explicit and implicit scheme [34], the design of MASC-N scheme we construct is as follows. Let $M-1=\left(2s+1\right)l=ql$, here $q$ and $l$ are positive integers $q,l\ge 3$, we can divide the points on each time level into $q$ sections which are sequentially recorded as $S_1,S_2,\cdots ,S_q$.

On the odd level, we apply explicit scheme (6) to calculate points $x_i\left(i={i^{″}}_1,{i^{\prime }}_2,{i^{″}}_3,{i^{\prime }}_4,\cdots ,{i^{″}}_{2s-1},{i^{\prime }}_{2s}\right)$, at points $x_i\left(i={i^{″}}_2,{i^{\prime }}_3,{i^{″}}_4,{i^{\prime }}_5,\cdots ,{i^{″}}_{2s},{i^{\prime }}_{2s+1}\right)$, we apply explicit scheme (7) to calculate. At the remaining points, we apply C-N scheme (8). When it turns to the even level, we apply explicit scheme (6) to calculate points $x_i\left(i={i^{″}}_2,{i^{\prime }}_3,{i^{″}}_4,{i^{\prime }}_5,\cdots ,{i^{″}}_{2s},{i^{\prime }}_{2s+1}\right)$, at points $x_i\left(i={i^{″}}_1,{i^{\prime }}_2,{i^{″}}_3,{i^{\prime }}_4,\cdots ,{i^{″}}_{2s-1},{i^{\prime }}_{2s}\right)$, we apply explicit scheme (7) to calculate. At the remaining points, we apply C-N scheme (8). The schematic design of the MASC-N scheme is shown in Figure 1.

Figure 1. Schematic design with segment nodes of the MASC-N scheme.

Then the approximate discrete (MASC-N) scheme of Equation (1) can be expressed in the matrix form:

$\{\begin{array}{l}\left(I+A_{1}G\right)V^{n+1}=-A_{2}G{V}^n+w_1{V}^n+\cdots +w_{n-1}{V}^2+l_n{V}^1+b^n,\\ \left(I+A_{2}G\right)V^{n+2}=-A_{1}G{V}^{n+1}+w_1{V}^{n+1}+\cdots +w_n{V}^2+l_{n+1}{V}^1+b^{n+1}.\end{array}$ (9)

where $n=1,3,\cdots$, $V^n={\left(\begin{array}{ccccc}{V}_2^n& V_3^n& \cdots & V_{M-1}^n& V_M^n\end{array}\right)}^{\text{T}}$, $b^n={\left(\begin{array}{ccccc}\frac{a}{2}\left(V_1^n+V_1^{n+1}\right)& 0& \cdots & 0& \frac{c}{2}\left(V_{M+1}^n+V_{M+1}^{n+1}\right)\end{array}\right)}^{\prime }$.

$A_1$ and $A_2$ are (M-1)-order diagonal matrices and they meet $A_1+A_2=I$, $A_1=diag\left({\theta }_1,{\theta }_2,\cdots ,{\theta }_{M-2},{\theta }_{M-1}\right)$.

${\theta }_i=\{\begin{array}{l}0,i={i^{″}}_{2l-1},{i^{\prime }}_{2l},0<l\le s\\ 1,i={i^{″}}_{2l},{i^{\prime }}_{2l+1},0<l\le s\\ \frac{1}{2},\text{elsewhere}\end{array}$, $G={\left(\begin{array}{ccccc}2b& -c& & & \\ -a& 2b& -c& & \\ & \ddots & \ddots & \ddots & \\ & & -a& 2b& -c\\ & & & -a& 2b\end{array}\right)}_{M-1}$.

Here I is an $\left(M-1\right)$ order unit matrix.

3. Theoretical Analysis of MASC-N Scheme

3.1. Existence and Uniqueness of the MASC-N Scheme Solution

Lemma 1. The matrices defined by the MASC-N scheme (9) are the nonnegative real matrices.

Proof.

$G+G^{\text{T}}=\left(\begin{array}{ccccc}\text{4}b& -\left(a+c\right)& & & \\ -\left(a+c\right)& \text{4}b& -\left(a+c\right)& & \\ & \ddots & \ddots & \ddots & \\ & & -\left(a+c\right)& \text{4}b& -\left(a+c\right)\\ & & & -\left(a+c\right)& \text{4}b\end{array}\right)$,

where $2b=\frac{m{\sigma }^2}{h^2}+mr$, $-\left(a+c\right)=-\frac{m{\sigma }^2}{h^2}$, $m=k^{\alpha }\Gamma \left(2-\alpha \right)$, we have $G+G^{\prime }$ as a strictly diagonally dominant three diagonal matrix, since the main diagonal elements are positive real numbers, $G+G^{\prime }$ is a positive definite matrix, $G+G^{\prime }$ is nonnegative definite matrix, G is nonnegative definite matrix.

Combined with lemma 1, ${\left(I+A_{1}G\right)}^{-1}$ and ${\left(I+A_{2}G\right)}^{-1}$ exist, and then we have following theorem:

Theorem 1. Solution of the MASC-N scheme (9) of time fractional Black-Scholes Equation (1) is existing and unique.

3.2. Stability of the MASC-N Scheme

Lemma 2 (Kellogg [35] [36]). If the matrix C is a nonnegative real matrix, that is, $C+C^{\text{T}}$ meet the nonnegative, then for any parameter $\rho \ge 0$, there are estimators $\Vert {\left(\rho I+C\right)}^{-1}\Vert \le {\rho }^{-1}$.

Lemma 3 ( [37]). The coefficients $w_j=l_j-l_{j+1}$ satisfy:

$1>w_1>w_2>\cdots >w_n>0$, ${\displaystyle \underset{j=1}{\overset{n-1}{\sum }}{w}_j}=1-l_n$, $l_1=1$.

Let ${\tilde{V}}_i^n$ be the approximate solution of Equation (1) at mesh point $\left(x_i,t_n\right)$. Defining ${\epsilon }_i^n={\tilde{V}}_i^n-V_i^n$, $E^n=\left({\epsilon }_2^n,{\epsilon }_3^n,\cdots ,{\epsilon }_M^n\right)$, $1\le n\le N+1$, and substituting $V_i^n={\tilde{V}}_i^n-{\epsilon }_i^n$ into MASC-N scheme, we have $E^n$ satisfying:

$\{\begin{array}{l}\left(I+A_{1}G\right)E^{n+1}=-A_{2}G{E}^n+w_1{E}^n+\cdots +w_{n-1}{E}^2+l_n{E}^1,\\ \left(I+A_{2}G\right)E^{n+2}=-A_{1}G{E}^{n+1}+w_1{E}^{n+1}+\cdots +w_n{E}^2+l_{n+1}{E}^1,\end{array}\text{}n=1,3,\cdots$ (10)

When $n\ge 3$, we have

$\begin{array}{c}{E}^{n+2}={\left(I+A_{2}G\right)}^{-1}\left(w_{1}I-A_{1}G\right){\left(I+A_{1}G\right)}^{-1}\left(w_{1}I-A_{2}G\right)E^n\\ +{\left(I+A_{2}G\right)}^{-1}\left(w_{1}I-A_{1}G\right){\left(I+A_{1}G\right)}^{-1}\left(w_2{E}^{n-1}+\cdots +w_{n-1}{E}^2+l_n{E}^1\right)\\ +{\left(I+A_{2}G\right)}^{-1}\left(w_2{E}^n+\cdots +w_n{E}^2+l_{n+1}{E}^1\right).\end{array}$ (11)

Take the norm of both sides of Equation (11)

$\begin{array}{c}\Vert E^{n+2}\Vert =\Vert {\left(I+A_{2}G\right)}^{-1}\left(w_{1}I-A_{1}G\right){\left(I+A_{1}G\right)}^{-1}\left(w_{1}I-A_{2}G\right)\Vert \Vert E^n\Vert \\ +\Vert {\left(I+A_{2}G\right)}^{-1}\left(w_{1}I-A_{1}G\right){\left(I+A_{1}G\right)}^{-1}\Vert \Vert \left(w_2{E}^{n-1}+\cdots +w_{n-1}{E}^2+l_n{E}^1\right)\Vert \\ +\Vert {\left(I+A_{2}G\right)}^{-1}\Vert \Vert \left(w_2{E}^n+\cdots +w_n{E}^2+l_{n+1}{E}^1\right)\Vert .\end{array}$

Due to the growth matrix $T={\left(I+A_{2}G\right)}^{-1}\left(w_{1}I-A_{1}G\right){\left(I+A_{1}G\right)}^{-1}\left(w_{1}I-A_{2}G\right)$ and we define $\tilde{T}=\left(I+A_{2}G\right)T{\left(I+A_{2}G\right)}^{-1}$, assume that the eigenvalue of G is r, then

$\begin{array}{c}\Vert T\Vert =\Vert \tilde{T}\Vert =\Vert \left(w_{1}I-A_{1}G\right){\left(I+A_{1}G\right)}^{-1}\left(w_{1}I-A_{2}G\right){\left(I+A_{2}G\right)}^{-1}\Vert \\ =\max \left\{\left|{\left(\frac{w_1-\theta r}{1+\theta r}\right)}^2\right|\right\}\le 1.\end{array}$

Then we use mathematical induction to prove $\Vert E^n\Vert \le \Vert E^1\Vert$.

When $n=1$, $\left(I+A_{1}G\right)E^2=\left(I-A_{2}G\right)E^1$

$\Vert E^2\Vert =\Vert {\left(I+A_{1}G\right)}^{-1}\left(I-A_{2}G\right)E^1\Vert \le \Vert E^1\Vert .$

When $\max \left\{w_1,\theta r\right\}=w_1$ is established,

$\begin{array}{c}\Vert E^3\Vert \le \Vert {\left(I+A_{2}G\right)}^{-1}\left(w_{1}I-A_{1}G\right){\left(I+A_{1}G\right)}^{-1}\left(I-A_{2}G\right)\Vert \Vert E^1\Vert +\Vert {\left(I+A_{2}G\right)}^{-1}\Vert \Vert l_2{E}^1\Vert \\ \le \max \left\{\left|\frac{w_1-\theta r}{1+\theta r}\right|\right\}\Vert E^1\Vert +\max \left\{\left|\frac{l_2}{1+\theta r}\right|\right\}\Vert E^1\Vert \\ \le \max \left\{\left|\frac{1-\theta r}{1+\theta r}\right|\right\}\Vert E^1\Vert \\ \le \Vert E^1\Vert .\end{array}$

When $\max \left\{w_1,\theta r\right\}=\theta r$ is established,

$\begin{array}{c}\Vert E^3\Vert \le \Vert {\left(I+A_{2}G\right)}^{-1}\left(w_{1}I-A_{1}G\right){\left(I+A_{1}G\right)}^{-1}\left(I-A_{2}G\right)\Vert \Vert E^1\Vert +\Vert {\left(I+A_{2}G\right)}^{-1}\Vert \Vert l_2{E}^1\Vert \\ \le \max \left\{\left|\frac{w_1-\theta r}{1+\theta r}\right|\right\}\Vert E^1\Vert +\max \left\{\left|\frac{l_2}{1+\theta r}\right|\right\}\Vert E^1\Vert \\ \le \max \left\{\frac{\theta r+l_2-w_1}{1+\theta r}\right\}\Vert E^1\Vert \\ \le \Vert E^1\Vert .\end{array}$

Assuming that $n\le k+1$, we have $\Vert E^n\Vert \le \Vert E^1\Vert$. When $n=k+2$, we prove it in two cases:

Case 1: $\max \left\{w_1-\theta r\right\}\le \max \left\{w_1,\theta r\right\}\le w_1$

$\begin{array}{c}\Vert E^{n+2}\Vert =\Vert {\left(I+A_{2}G\right)}^{-1}\left(w_{1}I-A_{1}G\right){\left(I+A_{1}G\right)}^{-1}\left(w_{1}I-A_{2}G\right)\Vert \Vert E^n\Vert \\ +\Vert {\left(I+A_{2}G\right)}^{-1}\left(w_{1}I-A_{1}G\right){\left(I+A_{1}G\right)}^{-1}\Vert \Vert \left(w_2{E}^{n-1}+\cdots +w_{n-1}{E}^2+l_n{E}^1\right)\Vert \\ +\Vert {\left(I+A_{2}G\right)}^{-1}\Vert \Vert \left(w_2{E}^n+\cdots +w_n{E}^2+l_{n+1}{E}^1\right)\Vert \\ \le {\left(\frac{w_1}{1+\theta r}\right)}^2\Vert E^1\Vert +\frac{w_1\left(1-w_1\right)}{{\left(1+\theta r\right)}^2}\Vert E^1\Vert +\frac{1-w_1}{1+\theta r}\Vert E^1\Vert \\ \le w_1\Vert E^1\Vert +\left(1-w_1\right)\Vert E^1\Vert \\ \le \Vert E^1\Vert .\end{array}$

Case 2: $\max \left\{w_1-\theta r\right\}\le \max \left\{w_1,\theta r\right\}\le \theta r$

$\begin{array}{c}\Vert E^{n+2}\Vert =\Vert {\left(I+A_{2}G\right)}^{-1}\left(w_{1}I-A_{1}G\right){\left(I+A_{1}G\right)}^{-1}\left(w_{1}I-A_{2}G\right)\Vert \Vert E^n\Vert \\ +\Vert {\left(I+A_{2}G\right)}^{-1}\left(w_{1}I-A_{1}G\right){\left(I+A_{1}G\right)}^{-1}\Vert \Vert \left(w_2{E}^{n-1}+\cdots +w_{n-1}{E}^2+l_n{E}^1\right)\Vert \\ +\Vert {\left(I+A_{2}G\right)}^{-1}\Vert \Vert \left(w_2{E}^n+\cdots +w_n{E}^2+l_{n+1}{E}^1\right)\Vert \end{array}$

$\begin{array}{c}\le {\left(\frac{\theta r}{1+\theta r}\right)}^2\Vert E^1\Vert +\frac{\theta r\left(1-w_1\right)}{{\left(1+\theta r\right)}^2}\Vert E^1\Vert +\frac{1-w_1}{1+\theta r}\Vert E^1\Vert \\ \le \frac{\theta r}{1+\theta r}\left(\frac{\theta r}{1+\theta r}+\frac{1-w_1}{1+\theta r}\right)\Vert E^1\Vert +\frac{1-w_1}{1+\theta r}\Vert E^1\Vert \\ \le \frac{\theta r}{1+\theta r}\Vert E^1\Vert +\frac{1-w_1}{1+\theta r}\Vert E^1\Vert \\ \le \Vert E^1\Vert .\end{array}$

To sum up, we obtain the following theorem.

Theorem 2. The MASC-N scheme (9) of time-fractional Black-Scholes Equation (1) is unconditional stability.

3.3. Convergence of the MASC-N Scheme

Let $v_i^n=v\left(x_i,t_n\right)\left(i=1,2,\cdots ,M+1;n=1,2,\cdots ,N+1\right)$ be the exact solution of Equation (1) at mesh point $\left(x_i,t_n\right)$. Defining $e_i^n=v_i^n-V_i^n$ and $e^n={\left(e_2^n,e_3^n,\cdots ,e_{M-1}^n,e_M^n\right)}^{\text{T}}$, and substituting $e_i^n=v_i^n-V_i^n$ into MASC-N scheme, we have $e^n$ satisfying:

$\{\begin{array}{l}\left(I+A_{1}G\right)e^{n+1}=-A_{2}G{e}^n+w_1{e}^n+\cdots +w_{n-1}{e}^2+l_n{e}^1+R^n,\\ \left(I+A_{2}G\right)e^{n+2}=-A_{1}G{e}^{n+1}+w_1{e}^{n+1}+\cdots +w_n{e}^2+l_{n+1}{e}^1+R^{n+1},\end{array}\text{}n=1,3,\cdots$ (12)

and $e^1=0,R^n=k^{\alpha }O\left(k^{2-\alpha }+h^2\right)$,which means there exists a positive constant C such that $\Vert R^n\Vert \le C{k}^{\alpha }\left(k^{2-\alpha }+h^2\right)$.

Imitate the proof of stability, when $n\ge 3$, we have

$\begin{array}{c}{e}^{n+2}={\left(I+A_{2}G\right)}^{-1}\left(w_{1}I-A_{1}G\right){\left(I+A_{1}G\right)}^{-1}\left(w_{1}I-A_{2}G\right)e^n\\ +{\left(I+A_{2}G\right)}^{-1}\left(w_{1}I-A_{1}G\right){\left(I+A_{1}G\right)}^{-1}\left(w_2{e}^{n-1}+\cdots +w_{n-1}{e}^2+R^n\right)\\ +{\left(I+A_{2}G\right)}^{-1}\left(w_2{e}^n+\cdots +w_n{e}^2+R^{n+1}\right).\end{array}$

When $n=1$, it is straightforward to see that:

$\Vert e^2\Vert \le \Vert {\left(I+A_{1}G\right)}^{-1}\Vert \Vert R^1\Vert \le l_1^{-1}{k}^{\alpha }C\left(k^{2-\alpha }+h^2\right).$

When $\max \left\{w_1,\theta r\right\}=w_1$ is established,

$\begin{array}{c}\Vert e^3\Vert \le \Vert {\left(I+A_{2}G\right)}^{-1}\left(w_{1}I-A_{1}G\right){\left(I+A_{1}G\right)}^{-1}\Vert \Vert R^1\Vert +\Vert {\left(I+A_{2}G\right)}^{-1}\Vert \Vert R^2\Vert \\ \le \left|\frac{w_1-\theta r}{{\left(1+\theta r\right)}^2}\right|C{k}^{\alpha }\left(k^{2-\alpha }+h^2\right)+\frac{1}{1+\theta r}C{k}^{\alpha }\left(k^{2-\alpha }+h^2\right)\\ \le \frac{w_1+1}{{\left(1+\theta r\right)}^2}C{k}^{\alpha }\left(k^{2-\alpha }+h^2\right)\\ \le l_2^{-1}C{k}^{\alpha }\left(k^{2-\alpha }+h^2\right).\end{array}$

When $\max \left\{w_1,\theta r\right\}=\theta r$ is established,

$\begin{array}{c}\Vert e^3\Vert \le \Vert {\left(I+A_{2}G\right)}^{-1}\left(w_{1}I-A_{1}G\right){\left(I+A_{1}G\right)}^{-1}\Vert \Vert R^1\Vert +\Vert {\left(I+A_{2}G\right)}^{-1}\Vert \Vert R^2\Vert \\ \le \left|\frac{\theta r}{1+\theta r}\right|C{k}^{\alpha }\left(k^{2-\alpha }+h^2\right)+\left|\frac{1}{1+\theta r}\right|C{k}^{\alpha }\left(k^{2-\alpha }+h^2\right)\le C{k}^{\alpha }\left(k^{2-\alpha }+h^2\right)\\ \le l_2^{-1}C{k}^{\alpha }\left(k^{2-\alpha }+h^2\right).\end{array}$

Assuming that $n\le k+1$, we have $\Vert e^n\Vert \le l_{n-1}^{-1}C{k}^{\alpha }\left(k^{2-\alpha }+h^2\right)$. When $n=k+2$, we prove it in two cases:

Case 1: $\max \left\{w_1-\theta r\right\}\le \max \left\{w_1,\theta r\right\}\le w_1$

$\begin{array}{c}\Vert e^{n+2}\Vert \le \Vert {\left(I+A_{2}G\right)}^{-1}\left(w_{1}I-A_{1}G\right){\left(I+A_{1}G\right)}^{-1}\left(w_{1}I-A_{2}G\right)\Vert \Vert e^n\Vert \\ +\Vert {\left(I+A_{2}G\right)}^{-1}\left(w_{1}I-A_{1}G\right){\left(I+A_{1}G\right)}^{-1}\Vert \Vert \left(w_2{e}^{n-1}+\cdots +w_{n-1}{e}^2\right)+R^n\Vert \end{array}$

$\begin{array}{c}+\Vert {\left(I+A_{2}G\right)}^{-1}\Vert \Vert \left(w_2{e}^n+\cdots +w_n{e}^2\right)+R^{n+1}\Vert \\ \le {\left(\frac{w_1}{1+\theta r}\right)}^2{l}_{n-1}^{-1}C{k}^{\alpha }\left(k^{2-\alpha }+h^2\right)+\left(\frac{w_1-\theta r}{{\left(1+\theta r\right)}^2}+\frac{1}{1+\theta r}\right)\\ \times \left({\displaystyle \underset{i=2}{\overset{n}{\sum }}{l}_{n+1}^{-1}{w}_i}+1\right)C{k}^{\alpha }\left(k^{2-\alpha }+h^2\right)\end{array}$

$\begin{array}{c}\le {\left(\frac{w_1}{1+\theta r}\right)}^2{l}_{n-1}^{-1}C{k}^{\alpha }\left(k^{2-\alpha }+h^2\right)+\frac{w_1+1}{{\left(1+\theta r\right)}^2}\left({\displaystyle \underset{i=2}{\overset{n}{\sum }}{w}_i}+l_{n+1}\right)l_{n+1}^{-1}C{k}^{\alpha }\left(k^{2-\alpha }+h^2\right)\\ \le w_1^2{l}_{n-1}^{-1}C{k}^{\alpha }\left(k^{2-\alpha }+h^2\right)+\left(1+w_1\right)\left(1-w_1\right)l_{n+1}^{-1}C{k}^{\alpha }\left(k^{2-\alpha }+h^2\right)\\ \le l_{n+1}^{-1}C{k}^{\alpha }\left(k^{2-\alpha }+h^2\right).\end{array}$

Case 2: $\max \left\{w_1-\theta r\right\}\le \max \left\{w_1,\theta r\right\}\le \theta r$

$\begin{array}{c}\Vert e^{n+2}\Vert =\Vert {\left(I+A_{2}G\right)}^{-1}\left(w_{1}I-A_{1}G\right){\left(I+A_{1}G\right)}^{-1}\left(w_{1}I-A_{2}G\right)\Vert \Vert e^n\Vert \\ +\Vert {\left(I+A_{2}G\right)}^{-1}\left(w_{1}I-A_{1}G\right){\left(I+A_{1}G\right)}^{-1}\Vert \Vert \left(w_2{e}^{n-1}+\cdots +w_{n-1}{e}^2\right)+R^n\Vert \\ +\Vert {\left(I+A_{2}G\right)}^{-1}\Vert \Vert \left(w_2{e}^n+\cdots +w_n{e}^2\right)+R^{n+1}\Vert \\ \le {\left(\frac{\theta r}{1+\theta r}\right)}^2{l}_{n-1}^{-1}C{k}^{\alpha }\left(k^{2-\alpha }+h^2\right)+\left(\frac{\theta r-w_1}{{\left(1+\theta r\right)}^2}+\frac{1}{1+\theta r}\right)\\ \times \left({\displaystyle \underset{i=2}{\overset{n}{\sum }}{l}_{n+1}^{-1}{w}_i}+1\right)C{k}^{\alpha }\left(k^{2-\alpha }+h^2\right)\end{array}$

$\begin{array}{c}\le {\left(\frac{\theta r}{1+\theta r}\right)}^2{l}_{n-1}^{-1}C{k}^{\alpha }\left(k^{2-\alpha }+h^2\right)+\frac{2\theta r+1}{{\left(1+\theta r\right)}^2}\left({\displaystyle \underset{i=2}{\overset{n}{\sum }}{l}_{n+1}^{-1}{w}_i}+1\right)C{k}^{\alpha }\left(k^{2-\alpha }+h^2\right)\\ \le {\left(\frac{\theta r}{1+\theta r}\right)}^2{l}_{n-1}^{-1}C{k}^{\alpha }\left(k^{2-\alpha }+h^2\right)+\frac{2\theta r+1}{{\left(1+\theta r\right)}^2}\left({\displaystyle \underset{i=2}{\overset{n}{\sum }}{w}_i}+l_{n+1}\right)l_{n+1}^{-1}C{k}^{\alpha }\left(k^{2-\alpha }+h^2\right)\\ \le {\left(\frac{\theta r}{1+\theta r}\right)}^2{l}_{n-1}^{-1}C{k}^{\alpha }\left(k^{2-\alpha }+h^2\right)+\frac{2\theta r+1}{{\left(1+\theta r\right)}^2}{l}_{n+1}^{-1}C{k}^{\alpha }\left(k^{2-\alpha }+h^2\right)\\ \le l_{n+1}^{-1}C{k}^{\alpha }\left(k^{2-\alpha }+h^2\right).\end{array}$

Due to

$\underset{n\to \infty }{\lim }\frac{l_n^{-1}}{n^{\alpha }}=\underset{n\to \infty }{\lim }\frac{n^{-\alpha }}{n^{1-\alpha }-{\left(n-1\right)}^{1-\alpha }}=\underset{n\to \infty }{\lim }\frac{n^{-1}}{1-{\left(1-\frac{1}{n}\right)}^{1-\alpha }}=\frac{1}{1-\alpha }.$

There is a positive constant C₁, such that

$\Vert e^n\Vert \le n^{\alpha }{C}_1\left(k^2+k^{\alpha }{h}^2\right)={\left(nk\right)}^{\alpha }{C}_1\left(k^{2-\alpha }+h^2\right),n=1,2,\cdots ,N.$

Here $C_1=\frac{C}{1-\alpha }$, it is clear that $nk<T$, then we can know that $\left|v\left(x_i,{\tau }_n\right)-V_i^n\right|\le C_2\left(k^{2-\alpha }+h^2\right)$ and $C_2=T^{\alpha }{C}_1$.

To sum up, we obtain the following theorem.

Theorem 3. The MASC-N scheme (9) of time-fractional Black-Scholes Equation (1) is convergent, and $\left|v\left(x_i,{\tau }_n\right)-V_i^n\right|\le C\left(k^{2-\alpha }+h^2\right)$, whereC is a positive constant.

4. Numerical Experiments

Numerical experiments will be done in Matlab 2013b, based on the Intel Core i3-2330 [e-mail] analytic solution of TFBSM is very complex to be obtained [13]. We use the MASC-N scheme (9), implicit numerical method (INM) [25] to calculate European call option prices governed by the time factional B-S model. To show the accuracy and effectiveness of the schemes numerically, we firstly take example 1 with the given terminal and boundary condition as an example.

Example 1 [25] [31]. Consider the following time fractional B-S model with homogeneous boundary conditions

$\{\begin{array}{l}{}_{0}D{}_{\tau }^{\alpha }U\left(x,\tau \right)=a\frac{{\partial }^{2}U\left(x,\tau \right)}{\partial x^2}+b\frac{\partial U\left(x,\tau \right)}{\partial x}-cU\left(x,\tau \right)+g\left(x,\tau \right),\\ U\left(0,\tau \right)=0,U\left(1,\tau \right)=0,\\ U\left(x,0\right)=x^2\left(1-x\right).\end{array}$

where the source term

$\begin{array}{l}g=\left(\frac{2{\tau }^{2-\alpha }}{\Gamma \left(3-\alpha \right)}+\frac{2{\tau }^{1-\alpha }}{\Gamma \left(2-\alpha \right)}\right)x^2\left(1-x\right)\\ -{\left(t+1\right)}^2\left[a\left(2-6x\right)+b\left(2x-3x^2\right)-c{x}^2\left(1-x\right)\right]\end{array}$

is chosen so that the exact solution of example 1 is $U={\left(t+1\right)}^2{x}^2\left(1-x\right)$. Here we take the parameters values as $r=0.05,\sigma =0.25,\alpha =0.7,a=\frac{1}{2}{\sigma }^2,b=r-a,c=r,T=1$.

Next, we verify the calculation precision and convergence order of MASC-N parallel difference schemes for solving time fractional B-S equation. Convergence order of spatial layer (COS) and convergence order of time horizon (COT) are defined respectively [38] [39]. $V_i^n$ is the exact solution and ${\bar{V}}_i^n$ is the solution of MASC-N.

$E_{\text{2},\Delta t}^m={\Vert V_j^m-{\bar{V}}_j^m\Vert }_2={\left\{{\displaystyle \underset{j=1}{\overset{n}{\sum }}{\left(V_j^m-{\bar{V}}_j^m\right)}^2\Delta t}\right\}}^{\frac{1}{2}},$

$E_{\infty ,\Delta t}^m={\Vert V_j^m-{\bar{V}}_j^m\Vert }_{\infty }=\underset{1\le j\le n}{\max }\left|V_j^m-{\bar{V}}_j^m\right|,$

$\text{COT}\approx {\log }_{\text{2}}\left(E_{\Delta x,\Delta t}^m/E_{\Delta x,\Delta t/2}^m\right),l=2,\infty ,$

$E_{\text{2},\Delta x}^n={\Vert V_i^n-{\bar{V}}_i^n\Vert }_2={\left\{{\displaystyle \underset{i=1}{\overset{m}{\sum }}{\left(V_i^n-{\bar{V}}_i^n\right)}^2\Delta x}\right\}}^{\frac{1}{2}},$

$E_{\infty ,\Delta x}^n={\Vert V_i^n-{\bar{V}}_i^n\Vert }_{\infty }=\underset{1\le i\le m}{\max }\left|V_i^n-{\bar{V}}_i^n\right|,$

$\text{COS}\approx {\log }_{\text{2}}\left(E_{\Delta x,\Delta t}^n/E_{\Delta x/2,\Delta t/4}^n\right),l=2,\infty .$

Table 1 and Table 2 list the numerical errors in 2-norm and the maximum-norm and their corresponding convergence rates. As can be seen from the table, the numerical solutions obtained by the MASC-N scheme are in excellent agreement with exact solutions. We can see that for fixed space step h = 1/100 and some different time steps k = 1/40, 1/80, 1/160 and 1/320, the corresponding orders of convergence for the MASC-N scheme are close to 1.3 when α = 0.7. The space convergence orders of the MASC-N scheme approach to 2. It shows good agreement with the conclusion of Theorem 3.

Example 2 [25]. Assuming that the current price of the underlying stock is 97 dollars, the strike price of option is 50 dollars, the risk-free interest rate is 0.05 per year, the deadline of the option is 12 months, the volatility is 0.25 per year.

First, let $M=101$, $N=100$, $l=200$, $q=5$. The parameters are $S=97$, $K=50$, $T=12$, $r=0.05$, $\sigma =0.25$, $M^{-}=\ln 0.1$, $M^{+}=\ln 100$.

We give the plots of the option price curves $P\left(S,T\right)$ under the INM, MASC-N scheme.

The parameter values are selected according to paper [25]. From Figure 2, we can see that the characteristics of Figure 2 are completely consistent with Figure 2 in paper [25] and the numerical solutions of MASC-N scheme are very close to that of implicit scheme.

Second, in order to verify the stability and the accuracy of MASC-N scheme, we regard the solution $V_i^n$ of INM as the control solution, and the solution ${\bar{V}}_i^n$ of MASC-N is used as the perturbation solution, then define the sum of relative error for every time level (SRET) and difference total energy (DTE):

$\text{SRET}\left(n\right)={\displaystyle \underset{i=1}{\overset{M}{\sum }}\frac{\left|V_i^n-{\bar{V}}_i^n\right|}{V_i^n}}$, $\text{DTE}\left(i\right)=\frac{1}{2}{\displaystyle \underset{n=1}{\overset{N}{\sum }}{\left(V_i^n-{\bar{V}}_i^n\right)}^2}$.

Table 1. Numerical errors and orders of convergence of MASC-N scheme for the example 1 when M = 101.

Table 2. Numerical errors and orders of convergence of MASC-N scheme for the example 1 when N = 100.

Figure 2. Call option prices obtained by MASC-N scheme and INM for α= 0.7.

Let $M=101$, $N=100$, other parameters are the same as above. The results of numerical experiments are as follows:

From Figure 3, we can see that the SRET curves of the MASC-N scheme are larger at the beginning and decrease along with the movement of time step and gradually keep a constant level. This shows that the MASC-N scheme has better stability.

Figure 4 suggests that the DTE of numerical solutions of the MASC-N scheme has the same order of magnitude compared with that of INM, showing that the accuracy of MASC-N scheme is close to that of INM. On the other hand, maximum value of DTE doesn’t exceed 0.014, indicating that the MASC-N scheme to solve the time fractional B-S equation have good accuracy.

From the viewpoint of computational efficiency, INM is the serial format and needs to use the catch-up method to solve three diagonal equations. The advantage of the MASC-N scheme of this paper is that they can divide the large numerical problems $AX=b$ into five small problems $A_{i}X=b,i=1,2,\cdots ,5$. Using “parfor” to implement parallel computing means that the cycle is divided into independent parts and the various parts can carry out parallel. It leads to the improvement of the efficiency of program execution and the effective reduction of computation time [40] [41].

Finally, in order to better reflect the superiority of parallel difference schemes, we choose INM, MASC-N scheme (9) to perform numerical tests. The number of time layers is 1000 (N = 1000), and the spatial grid points gradually increase from 201 to 1001. Next, the number of space grid points is 1001 (M = 1001), and the time layers gradually increase from 200 to 1000. And results are shown in Figure 5 and Table 3.

Obviously, MASC-N parallel scheme has superiority in saving the computing time with the increase of time layer as shown in Figure 5. Compared with the INM, the computing time of the MASC-N scheme can save nearly 90% when the space grid is the same (see in Table 3). With the encryption of the spatial grid, the computational efficiency of the MASC-N scheme is significantly higher than that of the INM. Especially when the number of spatial points is large enough, the MASC-N method has obvious localization characteristics in computing and communication.

Figure 3. The curves of SRET at time layer.

Figure 4. The curve of DTE at space layer.

Table 3. Comparison of two schemes’ computing time at M = 1001.

Example 3. Similarly, we consider the case that the underlying asset is a call option. Here we take the parameters values as $S=97$, $K=50$, $r=0.05$, $\sigma =0.25$ and $T=\text{6}$. Using INM, MASC-N difference scheme to compute the price of European call options with different values of α and the relative error when α = 0.7, the results are calculated as follows:

In general, option price obtained by the classical B-S equation (α = 1) is lower than real financial market price whose deadline is 6 months [42]. From Figure 6 and Table 4, the option price is a decreasing function of α, the option price calculated by time fractional B-S equation is higher than that of traditional B-S equation. We conclude that the TFBSM can more correctly capture the characteristics of a jump or large movement, compared with the classical B-S process. It is visible enough that we can get closer to the real financial market by properly selecting the value of α.

Figure 5. Comparison of the two schemes’ computing time.

Figure 6. Call option prices using MASC-N scheme for different values of α.

Table 4. Comparison of numerical results of several schemes for some values of α.

5. Conclusions

For the time fractional B-S equation with boundary conditions satisfied by standard European call options, this paper constructs the MASC-N parallel difference scheme. This scheme has been shown to be uniquely solved, unconditionally stable and convergent. And convergence rates are temporally 2-α order and spatially 2-order.

Moreover, compared with the implicit difference scheme, the MASC-N difference scheme is easy to realize parallel computing. So this scheme can greatly improve the computational efficiency, and it has obvious advantages in solving the multi-asset option pricing problems. All the results imply that MASC-N difference scheme is feasible to solve the time fractional B-S equation.

Acknowledgements

The work is supported by the Fundamental Research Funds for the Central Universities (No. 2021MS045).

Conflicts of Interest

The authors declare no conflicts of interest.

References

[1]	Kwok, Y. (2008) Mathematical Models of Financial Derivatives. 2nd Edition, Springer-Verlag, Berlin.
[2]	Jiang, L.S. (2008) Mathematical Models and Methods for Option Pricing. 2nd Edition, Higher Education Press, Beijing. (In Chinese)
[3]	Barles, G. and Soner H.M. (1998) Option Pricing with Transaction Costs and a Nonlinear Black-Scholes Equation. Finance and Stochastics, 2, 369-397. https://doi.org/10.1007/s007800050046
[4]	Davis, H.A.M. and Vassilios, G. (1993) European Option Pricing with Transaction Costs. SIAM Journal on Control & Optimization, 31, 470-493. https://doi.org/10.1137/0331022
[5]	Merton, R.C. (1976) Option When Underlying Stock Returns Are Discontinuous. Journal of Financial Economics, 3, 125-144. https://doi.org/10.1016/0304-405X(76)90022-2
[6]	Hull, J.C. and White, A.D. (1987) The Pricing of Options on Asset with Stochastic Volatilities. The Journal of Finance, 42, 281-300. https://doi.org/10.1111/j.1540-6261.1987.tb02568.x
[7]	Bjork, T. and Hult, H. (2005) A Note on Wick Products and the Fractional Black-Scholes Model. Finance and Stochastics, 9, 197-209. https://doi.org/10.1007/s00780-004-0144-5
[8]	Wang, X.T. (2010) Scaling and Long-Range Dependence in Option Pricing I: Pricing European Option with Transaction Costs under the Fractional Black-Scholes Model. Physica A Statistical Mechanics & Its Applications, 389, 438-444. https://doi.org/10.1016/j.physa.2009.09.041
[9]	Wyss, W. (2000) The Fractional Black-Scholes Equations. Fractional Calculus & Applied Analysis, 3, 51-61.
[10]	Cartea, A. and Del-Castillo-Negrete, D. (2007) Fractional Diffusion Models of Option Prices in Markets with Jumps. Physica A Statistical Mechanics & Its Applications, 374, 749-763. https://doi.org/10.1016/j.physa.2006.08.071
[11]	Jumarie, G. (2008) Stock Exchange Fractional Dynamics Defined as Fractional Exponential Growth Driven by Gaussian White Noise. Application to Fractional Black-Scholes Equations. Insurance Mathematics & Economics, 42, 271-287. https://doi.org/10.1016/j.insmatheco.2007.03.001
[12]	Jumarie, G. (2010) Derivation and Solutions of Some Fractional Black-Scholes Equations in Coarse-Grained Space and Time. Application to Merton’s Optimal Portfolio. Computers & Mathematics with Applications, 59, 1142-1164. https://doi.org/10.1016/j.camwa.2009.05.015
[13]	Liang, J.R., Wang, J., Zhang, W.J., et al. (2010) The Solution to a Bi-Fractional Black-Scholes-Merton Differential Equation. International Journal of Pure and Applied Mathematics, 58, 99-112.
[14]	Chen, W., Xu, X. and Zhu, S.P. (2015) Analytically Pricing Double Barrier Options Based on a Time-Fractional Black-Scholes Equation. Computers & Mathematics with Applications, 69, 1407-1419. https://doi.org/10.1016/j.camwa.2015.03.025
[15]	Elbeleze, A.A., Kiliçman, A. and Taib, B.M. (2013) Homotopy Perturbation Method for Fractional Black-Scholes European Option Pricing Equations Using Sumudu Transform. Mathematical Problems in Engineering, 2013, Article ID: 524852. https://doi.org/10.1155/2013/524852
[16]	Hariharan, G. (2013) An Efficient Wavelet Based Approximation Method to Time Fractional Black-Scholes European Option Pricing Problem Arising in Financial Market. Applied Mathematical Sciences, 7, 3445-3456. https://doi.org/10.12988/ams.2013.35261
[17]	Yu, J., Feng, Y. and Wang, X. (2022) Lie Symmetry Analysis and Exact Solutions of Time Fractional Black-Scholes Equation. International Journal of Financial Engineering, 9, Article ID: 2250023. https://doi.org/10.1142/S2424786322500232
[18]	Sun, Z.Z. and Gao, G.H. (2015) A Finite Difference Method for Fractional Differential Equations. Science Press, Beijing. (In Chinese)
[19]	Rezaei, M., Yazdanian, A.R., Ashrafi, A., et al. (2021) Numerical Pricing Based on Fractional Black-Scholes Equation with Time-Dependent Parameters under the CEV Model: Double Barrier Options. Computers & Mathematics with Applications, 90, 104-111. https://doi.org/10.1016/j.camwa.2021.02.021
[20]	She, M., Li, L., Tang, R. and Li, D. (2021) A Novel Numerical Scheme for a Time Fractional Black-Scholes Equation. Journal of Applied Mathematics and Computing, 66, 853-870. https://doi.org/10.1007/s12190-020-01467-9
[21]	Song, L.N. and Wang, W.G. (2013) Solution of the Fractional Black-Scholes Option Pricing Model by Finite Difference Method. Abstract & Applied Analysis, 2013, Article ID: 194286. https://doi.org/10.1155/2013/194286
[22]	Kumar, S. and Singh, D. (2014) Numerical Computation of Fractional Black-Scholes Equation Arising in Financial Market. Egyptian Journal of Basic & Applied Sciences, 1, 177-183. https://doi.org/10.1016/j.ejbas.2014.10.003
[23]	Bhowmik, S.K. (2014) Fast and Efficient Numerical Methods for an Extended Black-Scholes Model. Computers & Mathematics with Applications, 67, 636-654. https://doi.org/10.1016/j.camwa.2013.12.008
[24]	Yang, X.Z., Zhang, X. and Wu, L.F. (2015) A Kind of Difference Method for Time-Fractional Option Pricing Model. Applied Mathematics: A Journal of Chinese Universities, 30, 234-244. (In Chinese) https://doi.org/10.1186/s13662-015-0643-z
[25]	Zhang, H., Liu, F., Turner, I., et al. (2016) Numerical Solution of the Time Fractional Black-Scholes Model Governing European Options. Computers & Mathematics with Applications, 1, 1722-1783. https://doi.org/10.1016/j.camwa.2016.02.007
[26]	Phaochoo, P., Luadsong, A. and Aschariyaphotha, N. (2016) The Meshless Local Petrov-Galerkin Based on Moving Kriging Interpolation for Solving Fractional Black-Scholes Model. Journal of King Saud University-Science, 28, 111-117. https://doi.org/10.1016/j.jksus.2015.08.004
[27]	Nuugulu, S.M., Gideon, F. and Patidar, K.C. (2021) A Robust Numerical Scheme for a Time-Fractional Black-Scholes Partial Differential Equation Describing Stock Exchange Dynamics. Chaos Solitons & Fractals, 145, Article ID: 110753. https://doi.org/10.1016/j.chaos.2021.110753
[28]	An, X., Liu, F., Zheng, M., et al. (2021) A Space-Time Spectral Method for Time-Fractional Black-Scholes Equation. Applied Numerical Mathematics, 165, 152-166. https://doi.org/10.1016/j.apnum.2021.02.009
[29]	Yan, R.F., He, Y. and Zou, Q. (2021) A Difference Method with Parallel Nature for Solving Time-Space Fractional Black-Schole Model. Chaos Solitons Fractals, 151, Article ID: 111280. https://doi.org/10.1016/j.chaos.2021.111280
[30]	Sarboland, M. and Aminataei, A. (2022) On the Numerical Solution of Time Fractional Black-Scholes Equation. International Journal of Computer Mathematics, 99, 1736-1753. https://doi.org/10.1080/00207160.2021.2011248
[31]	Roul, P. and Goura, V.M.K. (2021) Prasad A Compact Finite Difference Scheme for Fractional Black-Scholes Option Pricing Model. Applied Numerical Mathematics, 166, 40-60. https://doi.org/10.1016/j.apnum.2021.03.017
[32]	Zhang, S. and Huang, J. (2022) A HIGH Order MQ Quasi-Interpolation Method for Time Fractional Black-Scholes Model. Acta Mathematica Scientia Series A, 42, 1496-1505. (In Chinese) http://121.43.60.238/sxwlxbA/CN/Y2022/V42/I5/1496
[33]	Guo, B.L., Pu, X.K. and Huang, F.H. (2011) Fractional Partial Differential Equations and Their Numerical Solutions. Science Press, Beijing. (In Chinese)
[34]	Zhang, S.C. (2010) Finite Difference Numerical Computation of Parabolic Equation. Science Press, Beijing. (In Chinese)
[35]	Tavakoli, R. and Davami, P. (2006) New Stable Group Explicit Finite Difference Method for Solution of Diffusion Equation. Applied Mathematics & Computation, 181, 1379-1386. https://doi.org/10.1016/j.amc.2006.03.005
[36]	Zhang, B.L., Yuan, G.X., Liu, X.P., et al. (1994) Parallel Finite Difference Methods for Partial Differential Equation. Science Press, Beijing. (In Chinese)
[37]	Liu, F., Zhuang, P. and Liu, Q. (2015) Numerical Methods and Applications of Fractional Partial Differential Equations. Science Press, Beijing. (In Chinese)
[38]	Gao, G.H. and Sun, Z.Z. (2011) A Compact Finite Difference Scheme for the Fractional Sub-Diffusion Equations. Journal of Computational Physics, 230, 586-595. https://doi.org/10.1016/j.jcp.2010.10.007
[39]	Zhang, Y.N., Sun, Z.Z. and Wu, H.W. (2011) Error Estimates of Crank-Nicolson-Type Difference Schemes for the Sub-Diffusion Equation. SIAM Journal on Numerical Analysis, 49, 2302-2322. https://doi.org/10.1137/100812707
[40]	Chi, X.B., Wang, Y.W., Wang, Y., et al. (2015) Parallel Computing and Implementation Technology. Science Press, Beijing. (In Chinese)
[41]	Liu, W. (2012) Practical Parallel Programming of Matlab. Beijing University of Aeronautics and Astronautics Press, Beijing. (In Chinese)
[42]	Car, P. and Wu, L.R. (2004) Time-Changed Lévy Processes and Option Pricing. Journal of Financial Economics, 71, 113-141. https://doi.org/10.1016/S0304-405X(03)00171-5