A Fast Computational Scheme for Solving the Temporal-Fractional Black–Scholes Partial Differential Equation

Abstract

In this work, we propose a fast scheme based on higher order discretizations on graded meshes for resolving the temporal-fractional partial differential equation (PDE), which benefits the memory feature of fractional calculus. To avoid excessively increasing the number of discretization points, such as the standard finite difference or meshfree methods, and, at the same time, to increase the efficiency of the solver, we employ discretizations on spatially non-uniform meshes with an attention on the non-smoothness area of the underlying asset. Therefore, the PDE problem is transformed to a linear system of algebraic equations. We perform numerical simulations to observe and check the behavior of the presented scheme in contrast to the existing methods.

1. Introduction

1.1. PDE Formulation

Many of the new practical problems in the literature [1,2] require the derivative not to be of integer order and, due to this, fractional differentiation in several senses in terms of definition must be dealt with in order to model the problems.

Let us assume that a derivative is to be defined by the (right) modified Riemann–Liouville (RL) notion. Then the fractional Black–Scholes (BS) price of an option u, by considering $T>0$ as the maturity time, can be written as follows (forward in time) [1,2]:

$$\begin{array}{cc}\hfill \frac{{\partial }^{\alpha }u(z,t)}{\partial t^{\alpha }}& =\frac{1}{2}{\sigma }^2{z}^2\frac{{\partial }^{2}u(z,t)}{\partial z^2}+(r-q)z\frac{\partial u(z,t)}{\partial z}-ru(z,t),z\ge 0,\hfill \end{array}$$

(1)

where the interest rate is r and $0<\alpha \le 1$. Moreover, $\alpha =2H$, is the real-valued Hurst exponent and $1>H>0$, which describes an exponent for the time series’s memory [3]. Clearly (1) will be simplified to the classic BS partial differential equation (PDE) by assuming $H=1/2$. Additionally, $\sigma$ and q are some positive constants. u is deterministic and stands for the fractional option price in an equity market.

For most settings under (1), we must solve numerically the PDE problem to the present and then find the present value of the financial derivative. Time-fractional financial PDEs describe models which have important applications in memory behaviors of long-time heavy tail decay for computational finance. This raises a need for fractional PDE problems and subsequently finding new effective methods are necessary.

1.2. Initial Conditions

The initial condition for the fractional BS option pricing problem (1) is written for the two common cases of put and call as follows (respectively) ([4] Chapter 1):

$$u(z,0)={(K-z)}^{+},$$

(2)

and

$$u(z,0)={(z-K)}^{+},$$

(3)

where K stands for the strike price.

1.3. Boundary Conditions

For the put and call cases, the side conditions can be given by, respectively, [5]:

$$\underset{z\to \infty }{lim}u(z,t)=0,u(0,t)=Kexp\{-rt\},$$

(4)

$$\underset{z\to \infty }{lim}u(z,t)=z_{max}exp(-qt)-Kexp(-rt),u(0,t)=0,$$

(5)

where $z_{max}$ is a large fixed constant.

1.4. Caputo Fractional Derivative

For the real continuous function g, the modified Riemann–Liouville (RL) fractional derivative in the fractional calculus can be defined as [6]:

$$g^{\left(\alpha \right)}\left(t\right):={\displaystyle \frac{d^{\alpha }g\left(t\right)}{d{t}^{\alpha }}}=\left\{\begin{array}{cc}{\displaystyle {\left(\mathsf{\Gamma }(-\alpha )\right)}^{-1}{\int }_0^t{(t-\zeta )}^{-\alpha -1}(-g\left(0\right)+g\left(\zeta \right))d\zeta ,}\hfill & \alpha <0,\hfill \\ {\displaystyle {\left(\mathsf{\Gamma }(-\alpha +1)\right)}^{-1}{\displaystyle \frac{d}{dt}}{\int }_0^t{(t-\zeta )}^{-\alpha }(-g\left(0\right)+g\left(\zeta \right))d\zeta ,}\hfill & 0<\alpha <1,\hfill \\ {\displaystyle {\displaystyle \frac{d^o\left(g^{(\alpha -o)}\left(t\right)\right)}{d{t}^o}},}\hfill & o\le \alpha <o+1,\hfill \end{array}\right.$$

(6)

where $\mathsf{\Gamma }(·)$ is the Gamma function and o is a constant. Recall that fractional differentiation in the Caputo sense is given as follows [7]:

$${}_0^C{D}_t^{\alpha }g\left(t\right)=\left({\int }_0^t{(t-\zeta )}^{-\alpha }{g}^{\prime }\left(\zeta \right)d\zeta \right){\left(\mathsf{\Gamma }(1-\alpha )\right)}^{-1},0<\alpha <1.$$

(7)

1.5. Time Discretization

Let the time discretization be given as $0=t_0<t_1<\cdots <t_n=T,$ with the equally-spaced step size $\Delta t$ also given. The author in [8] presented the L1 scheme for computing the RL fractional differentiation on the function g as follows:

$${}_0^R{D}_t^{\alpha }\left(g\left(t_n\right)\right)=\Delta t^{-\alpha }\sum _{j=0}^n{f}_{j,n}g\left(t_{n-j}\right)+\mathcal{O}(\Delta t^{-\alpha +2}),$$

(8)

wherein

$$\mathsf{\Gamma }(-\alpha +2)f_{j,n}=\left\{\begin{array}{cc}{\displaystyle {(j-1)}^{-\alpha +1}-j^{-\alpha +1},}\hfill & j=n,\hfill \\ {\displaystyle {(-1+j)}^{-\alpha +1}+{(j+1)}^{-\alpha +1}-2j^{-\alpha +1},}\hfill & 1\le j\le n-1,\hfill \\ {\displaystyle 1,}\hfill & j=0.\hfill \end{array}\right.$$

(9)

For equally-spaced stencils, the L1 scheme can also be written in the form below:

$$\begin{array}{cc}\hfill {}_0^C{D}_t^{\alpha }\left(u(z_i,t_{k+1})\right)& =\frac{\Delta t^{-\alpha }}{\mathsf{\Gamma }[2-\alpha ]}\sum _{j=0}^k\frac{u_{i,j+1}-u_{i,j}}{\Delta t}\left(-{\left(j\right)}^{1-\alpha }+{(1+j)}^{1-\alpha }\right).\hfill \end{array}$$

(10)

1.6. Background on Numerical Methods

Numerical methods for PDEs are necessary for challenging PDEs arising in real applications. Mostly, methods are constructed belonging to various categories such as finite differences (FD) and finite element methods (FEMs). FD methods discussed for such problems mainly in [9,10,11]. The meshfree RBF methods based on a shape parameter were also discussed for resolving (1) very recently in [12].

It is commented that it is possible to first impose a logarithmic transformation on (1) to obtain a PDE formulation which does not have variable coefficients anymore. However, this is not done here since this transformation encounters the logarithm of zero and an unbounded domain on two sides, which is harder to deal with than the original one which has unboundedness only from one side.

The main target of the present investigation is to follow the works [13,14,15] and propose a fast and efficient scheme for solving (1) but with an employment of graded meshes. The motivation behind choosing graded meshes is that the original fractional PDE problem has non-smoothness at the strike price and thus a non-uniform discretization of higher order will help improving the accuracy of the numerical results. The main motivation behind the designed numerical method is to have a fourth-order method for spatial discretization with an adaptation on the financial hot area of the problem. In fact, a main drawback in the already published works to tackle (1) is that they rely on second-order approximations on uniform meshes which can result in losing some accuracies for low number of discretization points specially around the strike price ($z=K$).

1.7. Paper’s Outline

The remaining parts of this manuscript are as follows. Section 2 provides a graded mesh, i.e., a nonuniform mesh that we need to discretize (1)–(3) adaptively with a clear focus on the financially significant area, at which the underlying asset tends to $z=K$. Section 3 describes a high-order discretization along the spatial variable of fourth order by employing five/six adjacent nodes since the fourth order must attain on a graded mesh. Section 4 investigates how the proposed solver can be constructed in sparse arrays and matrix forms as elegantly as possible. Section 5 is devoted to the numerical investigation of our high-order solver. Numerical experiments are given and discussed in detail therein. A summary of the work is furnished in Section 6.

2. Spatial Discretization Nodes

Numerical solution methods generally utilize the discretization methord to perform approximate calculations. When the interval is partitioned more finely, the calculation result is closer to the theoretical solution. To proceed, we first recall an algorithm to produce adaptive meshes with an emphasis on the non-smooth part of the initial condition. This is pursued by employing the inverse of hyperbolic function which and mapping some sets of points based on such function to have more number of discretization nodes on the hot area. Then, we construct the fourth-order FD approximations on these meshes using five adjacent nodes for estimating the first derivative and six adjacent nodes in computing the second derivative of a function.

We now assume that ${\left\{z_i\right\}}_{i=1}^m$ is a partition for $z\in [z_{min},z_{max}]$. Now, a well-resulted mesh can be given as follows [16]:

$$z_i=\mathsf{\Psi }\left({\mu }_i\right),i=1,2,\dots ,m,$$

(11)

where $m\gg 3$ and

$${\mu }_{max}={\mu }_m>\cdots >{\mu }_2>{\mu }_1={\mu }_{min},$$

(12)

are m equidistant nodes and we also use:

$$\begin{array}{cc}\hfill & {\mu }_{\mathrm{int}}=\frac{z_{\mathrm{right}}-z_{\mathrm{left}}}{d_1},\hfill \\ \hfill & {\mu }_{min}={sinh}^{-1}\left(\frac{z_{min}-z_{\mathrm{left}}}{d_1}\right),\hfill \\ \hfill & {\mu }_{max}={\mu }_{\mathrm{int}}+{sinh}^{-1}\left(\frac{z_{max}-z_{\mathrm{right}}}{d_1}\right).\hfill \end{array}$$

(13)

We too take into account that $z_{min}=0$ and $z_{max}=3K$. The parameter ${\upsilon }_1>0$ controls the node density around $z=K$. In addition, one defines

$$\mathsf{\Psi }\left(\mu \right)=\left\{\begin{array}{cc}{z}_{\mathrm{left}}+{\upsilon }_{1}sinh\left(\mu \right),\hfill & {\mu }_{min}\le \mu <0,\hfill \\ z_{\mathrm{left}}+{\upsilon }_1\mu ,\hfill & 0\le \mu \le {\mu }_{\mathrm{int}},\hfill \\ z_{\mathrm{right}}+{\upsilon }_{1}sinh(\mu -{\mu }_{\mathrm{int}}),\hfill & {\mu }_{\mathrm{int}}<\mu \le {\mu }_{max}.\hfill \end{array}\right.$$

(14)

The idea behind the definition of $\mathsf{\Psi }\left(\mu \right)$ is to use hyperbolic functions in order to stretch the points. Using such a function and its application as a zooming function, we can stretch the discretization points. Here in (13), we consider ${\upsilon }_1=\frac{K}{20},$ while $z_{\mathrm{left}}=max\{e^{-0.0025T},$$0.5\}\times K$, $z_{\mathrm{right}}=K$, and $[z_{\mathrm{left}},z_{\mathrm{right}}]\subset [0,z_{max}]$,16].

3. A Fast High-Order Discretization

Let us consider the graded mesh (14) along the spatial independent variable of the PDE problem. To derive and introduce FD approximation on such graded meshes, five and six stencil nodes must be taken into account at each position of the computational discretized domain to finally attain approximated FD formulas which are of higher order of convergence, four, for both the first and second derivatives of the function.

Accordingly, without loss of generality, we take into consideration that $g\left(z\right)$ is a sufficiently differentiable function and there is a given set of discretized points as follows:

$$\{z_1,z_2,\cdots ,z_{m-1},z_m\}.$$

(15)

Here (15) is an arbitrary discretization of points and thus is it a graded mesh. In practice, it is produced by (11)–(14). Now five nodal points are considered as follows:

$$\{\{z_{i-2},g\left(z_{i-2}\right)\},\{z_{i-1},g\left(z_{i-1}\right)\},\{z_i,g\left(z_i\right)\},\{z_{i+1},g\left(z_{i+1}\right)\},\{z_{i+2},g\left(z_{i+2}\right)\}\},$$

(16)

and we calculate the interpolatory polynomial passing through the adjacent points; see e.g., [17]. By employing $z=z_i$ and some symbolic computations, one obtains the following FD formulation on our graded mesh:

$$g^{\prime }\left(z_i\right)=a_{i-2}g\left(z_{i-2}\right)+a_{i-1}g\left(z_{i-1}\right)+a_{i}g\left(z_i\right)+a_{i+1}g\left(z_{i+1}\right)+a_{i+2}g\left(z_{i+2}\right)+\mathcal{O}\left(h^4\right),$$

(17)

wherein the maximum spacing of local grid is h. Finally we attain

$$\begin{array}{cc}\hfill a_{i-2}& =-\frac{{\varrho }_{i-1,i}{\varrho }_{i,i+1}{\varrho }_{i,i+2}}{{\varrho }_{i-2,i-1}{\varrho }_{i-2,i}{\varrho }_{i-2,i+1}{\varrho }_{i-2,i+2}},\hfill \end{array}$$

(18)

$$\begin{array}{cc}\hfill a_{i-1}& =\frac{{\varrho }_{i-2,i}{\varrho }_{i,i+1}{\varrho }_{i,i+2}}{{\varrho }_{i-2,i-1}{\varrho }_{i-1,i}{\varrho }_{i-1,i+1}{\varrho }_{i-1,i+2}},\hfill \end{array}$$

(19)

$$\begin{array}{cc}\hfill a_i& =\frac{1}{{\varrho }_{i-2,i}{\varrho }_{i,i-1}{\varrho }_{i,i+1}{\varrho }_{i,i+2}}{Z}_1,\hfill \end{array}$$

(20)

$$\begin{array}{cc}\hfill a_{i+1}& =\frac{{\varrho }_{i-2,i}{\varrho }_{i,i-1}{\varrho }_{i,i+2}}{{\varrho }_{i-2,i+1}{\varrho }_{i+1,i-1}{\varrho }_{i+1,i}{\varrho }_{i+1,i+2}},\hfill \end{array}$$

(21)

$$\begin{array}{cc}\hfill a_{i+2}& =\frac{{\varrho }_{i-2,i}{\varrho }_{i,i-1}{\varrho }_{i,i+1}}{{\varrho }_{i-2,i+2}{\varrho }_{i+2,i-1}{\varrho }_{i+2,i}{\varrho }_{i+2,i+1}},\hfill \end{array}$$

(22)

using ${\varrho }_{l,q}=z_l-z_q$ and

$$\begin{array}{cc}\hfill Z_1=& z_{i-2}(z_{i-1}({\varrho }_{i+1,i}+{\varrho }_{i+2,i})\hfill \\ \hfill & +3z_i^2-2(z_{i+2}+z_{i+1})z_i+z_{i+1}{z}_{i+2})+z_i(-4z_i^2+3(z_{i+1}+z_{i+2})z_i\hfill \\ \hfill & -2z_{i+1}{z}_{i+2})+z_{i-1}(3z_i^2-2(z_{i+2}+z_{i+1})z_i+z_{i+1}{z}_{i+2}).\hfill \end{array}$$

(23)

The FD formulations (of the sided type) must be derived similarly and employed for the points on the boundaries and only the one near to the boundaries to have a unified fourth order of convergence. It is noted that, in the case of applying these formulas on structured grids, the fourth order of convergence is preserved.

Using a similar approach, we may obtain the weighting coefficients for approximating the 2nd derivative of the function using graded meshes and having fourth order of convergence. Hence, it is considered that:

$$\{\{z_{i-3},g\left(z_{i-3}\right)\},\{z_{i-2},g\left(z_{i-2}\right)\},\{z_{i-1},g\left(z_{i-1}\right)\},\{z_i,g\left(z_i\right)\},\{z_{i+1},g\left(z_{i+1}\right)\},\{z_{i+2},g\left(z_{i+2}\right)\}\},$$

(24)

and the 2nd-derivative interpolatory function $g^{″}\left(z\right)$ is calculated based on z. At this moment, by considering $z=z_i$, one obtains [17]:

$$\begin{array}{cc}\hfill g^{″}\left(z_i\right)=& b_{i-3}g\left(z_{i-3}\right)+b_{i-2}g\left(z_{i-2}\right)+b_{i-1}g\left(z_{i-1}\right)\hfill \\ \hfill & +b_{i}g\left(z_i\right)+b_{i+1}g\left(z_{i+1}\right)+b_{i+2}g\left(z_{i+2}\right)+\mathcal{O}\left(h^4\right),\hfill \end{array}$$

(25)

wherein

$$\begin{array}{cc}\hfill b_{i-3}& =\frac{Z_2}{{\varrho }_{i-3,i-2}{\varrho }_{i-3,i-1}{\varrho }_{i-3,i}{\varrho }_{i-3,i+1}{\varrho }_{i-3,i+2}},\hfill \end{array}$$

(26)

$$\begin{array}{cc}\hfill b_{i-2}& =\frac{Z_3}{{\varrho }_{i-3,i-2}{\varrho }_{i-2,i-1}{\varrho }_{i-2,i}{\varrho }_{i-2,i+1}{\varrho }_{i-2,i+2}},\hfill \end{array}$$

(27)

$$\begin{array}{cc}\hfill b_{i-1}& =\frac{Z_4}{{\varrho }_{i-2,i-1}{\varrho }_{i-1,i-3}{\varrho }_{i-1,i}{\varrho }_{i-1,i+1}{\varrho }_{i-1,i+2}},\hfill \end{array}$$

(28)

$$\begin{array}{cc}\hfill b_i& =\frac{Z_5}{{\varrho }_{i-3,i}{\varrho }_{i,i-2}{\varrho }_{i,i-1}{\varrho }_{i,i+1}{\varrho }_{i,i+2}},\hfill \end{array}$$

(29)

$$\begin{array}{cc}\hfill b_{i+1}& =\frac{Z_6}{{\varrho }_{i-3,i+1}{\varrho }_{i+1,i-2}{\varrho }_{i+1,i-1}{\varrho }_{i+1,i}{\varrho }_{i+1,i+2}},\hfill \end{array}$$

(30)

$$\begin{array}{cc}\hfill b_{i+2}& =\frac{Z_7}{{\varrho }_{i-3,i+2}{\varrho }_{i+2,i-2}{\varrho }_{i+2,i-1}{\varrho }_{i+2,i}{\varrho }_{i+2,i+1}}.\hfill \end{array}$$

(31)

Here, we obtain

$$\begin{array}{cc}\hfill Z_2=& z_{i-1}(-6z_i^2+4(z_{i+2}+z_{i+1})z_i-2z_{i+1}{z}_{i+2})+2z_i(4z_i^2\hfill \\ \hfill & -3(z_{i+2}+z_{i+1})z_i+2z_{i+1}{z}_{i+2})+z_{i-2}(-6z_i^2\hfill \\ \hfill & +4(z_{i+1}+z_{i+2})z_i-2z_{i+1}{z}_{i+2}+z_{i-1}(4z_i-2(z_{i+1}+z_{i+2}))),\hfill \\ \hfill Z_3=& 2(z_{i-3}(z_{i-1}({\varrho }_{i+1,i}+{\varrho }_{i+2,i})+3z_i^2-2(z_{i+1}\hfill \\ \hfill & +z_{i+2})z_i+z_{i+1}{z}_{i+2})+z_i(-4z_i^2+3(z_{i+1}+z_{i+2})z_i-2z_{i+1}{z}_{i+2})\hfill \\ \hfill & +z_{i-1}(3z_i^2-2(z_{i+1}+z_{i+2})z_i+z_{i+1}{z}_{i+2})),\hfill \\ \hfill Z_4=& 2(z_{i-3}(z_{i-2}({\varrho }_{i+1,i}+{\varrho }_{i+2,i})+3z_i^2-2(z_{i+1}\hfill \\ \hfill & +z_{i+2})z_i+z_{i+1}{z}_{i+2})+z_i(-4z_i^2+3(z_{i+1}+z_{i+2})z_i-2z_{i+1}{z}_{i+2})\hfill \\ \hfill & +z_{i-2}(3z_i^2-2(z_{i+1}+z_{i+2})z_i+z_{i+1}{z}_{i+2})),\hfill \\ \hfill Z_5=& 2(z_{i-2}(z_{i-1}({\varrho }_{i+1,i}+{\varrho }_{i+2,i}-z_i)+6z_i^2-3(z_{i+1}\hfill \\ \hfill & +z_{i+2})z_i+z_{i+1}{z}_{i+2})+z_{i-3}\left(z_{i-1}\right({\varrho }_{i+1,i}\hfill \\ \hfill & +{\varrho }_{i+2,i}-z_i)+z_{i-2}({\varrho }_{i+1,i}+{\varrho }_{i+2,i}\hfill \\ \hfill & +z_{i-1}-z_i)+6z_i^2-3z_{i+1}{z}_i-3z_{i+2}{z}_i+z_{i+1}{z}_{i+2})\hfill \\ \hfill & +z_i(2z_i(3z_{i+1}-5z_i)+z_{i-1}(6z_i-3z_{i+1}))\hfill \\ \hfill & +(3z_i(2z_i-z_{i+1})+z_{i-1}(z_{i+1}-3z_i))z_{i+2}),\hfill \end{array}$$

(32)

$$\begin{array}{cc}\hfill Z_6=& 2(z_i(z_{i-1}(2{\varrho }_{i,i+2}+z_i)+z_i(3z_{i+2}-4z_i))+z_{i-2}(z_i(2{\varrho }_{i,i+2}+z_i)\hfill \\ \hfill & +z_{i-1}({\varrho }_{i+2,i}-z_i))+z_{i-3}(2z_i{\varrho }_{i,i+2}\hfill \\ \hfill & +z_{i-1}({\varrho }_{i+2,i}-z_i)+z_{i-2}({\varrho }_{i-1,i}+{\varrho }_{i+2,i})+z_i^2\left)\right),\hfill \\ \hfill Z_7=& 2(z_i(z_{i-1}(2{\varrho }_{i,i+1}+z_i)+z_i(3z_{i+1}-4z_i))+z_{i-2}\left(z_i\right(2{\varrho }_{i,i+1}\hfill \\ \hfill & +z_i)+z_{i-1}({\varrho }_{i+1,i}-z_i))+z_{i-3}(2z_i{\varrho }_{i,i+1}\hfill \\ \hfill & +z_{i-1}({\varrho }_{i+1,i}-z_i)+z_{i-2}({\varrho }_{i-1,i}+{\varrho }_{i+1,i})+z_i^2\left)\right).\hfill \end{array}$$

Here, in fact, the FD approximation (25) based on Taylor expansions of appropriate orders and severe simplifications were performed on graded meshes. Note that, for the nodes $\{z_1,z_2,z_3,z_{m-1},z_m\}$ and to keep the quartical speed, a similar procedure can be performed by the six adjacent nodes as discussed above but for that specific node. The relations (17) and (25) are applicable on both graded and uniform meshes and, in the case of equidistant discretion nodes, can yield in more accurate formulations.

4. Construction of the Solver

By denoting $u_{i,j}\simeq u(z_i,t_j),$ one can find now the whole discretization solver for (1) as follows:

$$\begin{array}{cc}\hfill D_t^{\alpha }u(z_i,t_{j+1})+\mathcal{O}(\Delta t^{2-\alpha })=& \frac{1}{2}{\sigma }^2{z}_i^2[b_{i-3,j}{u}_{i-3,j}+b_{i-2,j}{u}_{i-2,j}+b_{i-1,j}{u}_{i-1,j}\hfill \\ \hfill & +b_{i,j}{u}_{i,j}+b_{i+1,j}{u}_{i+1,j}+b_{i+2,j}{u}_{i+2,j}]+\mathcal{O}\left(h^4\right)\hfill \\ \hfill & +(r-q)z_i[a_{i-2,j}{u}_{i-2,j}+a_{i-1,j}{u}_{i-1,j}+a_{i,j}{u}_{i,j}\hfill \\ \hfill & +a_{i+1,j}{u}_{i+1,j}+a_{i+2,j}{u}_{i+2,j}]+\mathcal{O}\left(h^4\right)\hfill \\ \hfill & -r{u}_{i,j}.\hfill \end{array}$$

(33)

It is now possible to write down

$$D_t^{\alpha }u(z_i,t_{j+1})={\left(\mathsf{\Gamma }[2-\alpha ]\right)}^{-1}\sum _{l=0}^j\frac{u_{i,l+1}-u_{i,l}}{\Delta t_{l+1}}\left({(-t_l+t_{j+1})}^{1-\alpha }-{(t_{j+1}-t_{l+1})}^{1-\alpha }\right).$$

(34)

To continue more formally, now the derivative matrices corresponding to the weights obtained by employing fourth-order approximations in Section 3 are gathered up in some matrices as follows:

$$M_z=\left\{\begin{array}{ccc}{a}_{i,j}\mathrm{by}\left(18\right)\hfill & \hfill & i-j=2,\hfill \\ a_{i,j}\mathrm{by}\left(19\right)\hfill & \hfill & i-j=1,\hfill \\ a_{i,j}\mathrm{by}\left(20\right)\hfill & \hfill & i-j=0,\hfill \\ a_{i,j}\mathrm{by}\left(21\right)\hfill & \hfill & j-i=1,\hfill \\ a_{i,j}\mathrm{by}\left(22\right)\hfill & \hfill & j-i=2,\hfill \\ 0\hfill & \hfill & \mathrm{otherwise},\hfill \end{array}\right.$$

(35)

and

$$M_{zz}=\left\{\begin{array}{ccc}{b}_{i,j}\mathrm{by}\left(26\right)\hfill & \hfill & i-j=3,\hfill \\ b_{i,j}\mathrm{by}\left(27\right)\hfill & \hfill & i-j=2,\hfill \\ b_{i,j}\mathrm{by}\left(28\right)\hfill & \hfill & i-j=1,\hfill \\ b_{i,j}\mathrm{by}\left(29\right)\hfill & \hfill & i-j=0,\hfill \\ b_{i,j}\mathrm{by}\left(30\right)\hfill & \hfill & j-i=1,\hfill \\ b_{i,j}\mathrm{by}\left(31\right)\hfill & \hfill & j-i=2,\hfill \\ 0\hfill & \hfill & \mathrm{otherwise}.\hfill \end{array}\right.$$

(36)

We now denote the time derivative matrix as $M_t$. Hence, to derive our proposed method using matrix notations, we consider the following matrix:

$$\mathrm{Y}=\frac{1}{2}{\sigma }^2{Z}^2(M_{zz}\otimes I_t)+(r-q)Z(M_z\otimes I_t)-r{I}_N,$$

(37)

where ⊗ is the Kronecker product and

$$Z=\mathrm{diag}\left(z_1,z_2,\cdots ,z_m\right)\otimes I_t.$$

(38)

Consider $A_{m\times n}$ and $B_{p\times q}$ and two matrices, then their Kronecker product $C=A\otimes B$ is an $\left(mp\right)\times \left(nq\right)$ matrix as followss [18]:

$$c_{\alpha ,\beta }=a_{i,j}{b}_{k,l},$$

(39)

wherein

$$\alpha =p(i-1)+k,\beta =q(j-1)+l.$$

(40)

Here $I_N=I_z\otimes I_t,$ is an $N\times N$ identity matrix, with $I_z$ and $I_t$ being $m\times m$ and $n\times n$ identity matrices, respectively. Thus, the temporal-fractional PDE (33) is finally fully discretized as follows:

$$(I_z\otimes M_t)U={\rm Y}U,$$

(41)

where $U={(u_{1,1},u_{1,2},\dots ,u_{m,n})}^{*}$. Now we take into account

$$\mathcal{B}=(I_z\otimes M_t)-{\rm Y}.$$

(42)

Thus, one obtains $\mathcal{B}U=0$. By incorporating the boundary and initial conditions, one is able to find a system of linear algebraic equations:

$$\overline{\mathcal{B}}U=\overline{b},$$

(43)

wherein $\overline{b}$ and $\overline{\mathcal{B}}$ are the vector and system matrix, respectively, after incorporating the conditions related to the initial and side.

Based on truncation errors obtained using the approximations used along the spatial and temporal variables of the PDE problem (1), the convergence order is $\mathcal{O}(\Delta t^{2-\alpha })+\mathcal{O}\left(h^4\right)$. Moreover, the consistency of the solving procedure can be obtained when $\Delta t,h$ tend to zero.

5. Numerical Tests

The major point here is to compare and exhibit the efficiency of different schemes to solve (1) on the same numerical domain. The solvers are summarized as follows:

• The quadratically convergent FD method presented in [19] and shown by FDM.
• The RBF-FD method on graded meshes given in [20] and denoted by SM.
• The presented scheme given in this work and shown by PMBS via the non-uniform mesh in Section 3.

Mathematica 12.0 [21] has been used to perform the calculations and write the programs. The computer specifications are Windows 11, Core i7-8750h. Here, $z_{max}=3K$ is considered for the numerical domains and the computational time is given in second. Furthermore, the absolute error can be calculated as:

$$\epsilon =\left|u_{\mathrm{approx}}-u_{\mathrm{ref}}\right|,$$

(44)

wherein $u_{\mathrm{ref}}$ and $u_{\mathrm{approx}}$ are the referenced and numerical solutions, respectively.

Example 1.

Let us consider the parameters below for solving (1) in the call case:

$$\alpha =0.8,\sigma =40\%,K=100\$,q=0,T=1\mathrm{year},r=5\%.$$

(45)

The reference value is $u(K,T)\simeq 18.1883$, [20].

Table 1. Convergence history for Example 1.

m,n	${\mathit{u}}_{\mathit{FDM}}$	${\mathit{\epsilon }}_{\mathit{FDM}}$	T${}_{\mathbf{FDM}}$	${\mathit{u}}_{\mathit{SM}}$	${\mathit{\epsilon }}_{\mathit{SM}}$	T${}_{\mathbf{SM}}$	${\mathit{u}}_{\mathit{PMBS}}$	${\mathit{\epsilon }}_{\mathit{PMBS}}$	T${}_{\mathbf{PMBS}}$
10	16.070	2.11 × 10${}^0$	0.01	17.378	8.1 × 10${}^{-1}$	0.01	17.381	8.0 × 10${}^{-1}$	0.01
20	17.889	2.98 × 10${}^{-1}$	0.02	17.694	4.94 × 10${}^{-1}$	0.02	17.701	4.8 × 10${}^{-1}$	0.02
40	17.913	2.75 × 10${}^{-1}$	0.16	18.003	1.85 × 10${}^{-1}$	0.18	18.091	9.7 × 10${}^{-2}$	0.16
80	18.039	1.48 × 10${}^{-1}$	3.96	18.213	2.46 × 10${}^{-2}$	3.84	18.205	1.6 × 10${}^{-2}$	3.63
120	18.055	1.33 × 10${}^{-1}$	20.02	18.190	1.64 × 10${}^{-3}$	20.96	18.187	1.3 × 10${}^{-3}$	19.17

consists of the convergence history of different solvers. Based on the observations, the proposed high-order method is efficient and gives higher accuracies in lower computational time. Additionally, to check the numerical stability of the presented scheme in practice, the computational solution attained for Example 1 is given in Figure 1 for various numbers of discretization points.

Example 2.

Let us consider the parameters below for solving (1) in the put case:

$$\alpha =0.7,\sigma =30\%,q=20\%,T=2\mathrm{year},r=5\%,K=100\$.$$

(46)

The reference solution is taken from [20] as $u(K,T)\simeq 24.6299$.

Numerical pieces of evidence have been furnished in

Table 2. Convergence history for Example 2.

m,n	${\mathit{u}}_{\mathit{FDM}}$	${\mathit{\epsilon }}_{\mathit{FDM}}$	T${}_{\mathbf{FDM}}$	${\mathit{u}}_{\mathit{SM}}$	${\mathit{\epsilon }}_{\mathit{SM}}$	T${}_{\mathbf{SM}}$	${\mathit{u}}_{\mathit{PMBS}}$	${\mathit{\epsilon }}_{\mathit{PMBS}}$	T${}_{\mathbf{PMBS}}$
11	23.914	7.1 × 10${}^{-1}$	0.01	23.631	9.9 × 10${}^{-1}$	0.03	23.712	9.1 × 10${}^{-1}$	0.02
21	24.316	3.1 × 10${}^{-1}$	0.03	24.117	5.1 × 10${}^{-1}$	0.05	24.239	3.9 × 10${}^{-1}$	0.05
41	24.466	1.6 × 10${}^{-1}$	0.15	24.356	2.7 × 10${}^{-1}$	0.21	24.546	8.3 × 10${}^{-2}$	0.19
81	24.545	8.4 × 10${}^{-2}$	4.14	24.565	6.4 × 10${}^{-2}$	4.69	24.601	2.8 × 10${}^{-2}$	4.37
161	24.580	4.9 × 10${}^{-2}$	63.97	24.636	6.1 × 10${}^{-3}$	64.01	24.624	5.9 × 10${}^{-3}$	62.64

. They too reveal that the presented scheme converges faster and is more effective for solving the temporal fractional BS PDE. For Example 2, the numerical solutions for various numbers of discretization points are given in Figure 2 as well.

Example 3.

In this test, consider the same assumptions as in Example 2 but different $\alpha =\{0.2,0.5,0.8\}$ to observe how the numerical solutions behave.

This is pursued by imposing the proposed high-order procedure of this work with 21 discretization points, which is 20 un-equally spaced discretized nodes. The nodes are obtained via (11). The numerical observations are gathered in Figure 3. It is revealed that by increasing $\alpha$ the price of the option increases especially on the hotzone (around the strike price).

6. Conclusions

One of the fundamental and important problems in mathematical finance is to price options via the PDE of BS. Several extensions to this model were given in the literature, such as the model of Heston. However extensions based on the fractional Brownian motion which leads to fractional BS PDE are of practical importance since the memory feature of fractional derivatives can help obtain more reliable option prices in equity markets. To achieve this goal, we have proposed a high-order discretization scheme for solving the fractional BS PDE using graded meshes, i.e., on non-uniform meshes. The motivation of using this was to concentrate more on the hotzone. In this area, the initial condition has a payoff. We have presented our solver on sparse matrices and by the help of the Kronecker product, we have transformed the whole PDE into a system of algebraic linear equations. Numerical tests have been carried out and confirmed the superiority of the novel high-order method for solving fractional BS PDE. Table 1 and Table 2 show the convergence history of different solvers while revealing that our proposed solver is better than the FDM and the SM schemes. Forthcoming studies can be focused on the use of graded meshes on rough Heston PDE along with high-order discretizations to propose a fast solver for such a problem.