Two high-order compact difference schemes with temporal graded meshes for time-fractional Black-Scholes equation

1.   Introduction

Fractional differential equations, which can be regarded as the generalization of classical differential equations that involve non-integer order derivatives, have emerged as a powerful tool in modeling and describing many complex phenomena that may not be captured by classical differential equations. These significant applications in a wide range of scientific and engineering fields include anomalous diffusion process, viscoelastic materials, chemistry, economics and finance. Especially for the application in finance, e.g., the modeling of stock prices and financial derivatives, it seems that fractional differential equations have more potential for the description of memory and hereditary properties in the market, see ^[4,12] and the references therein for more discussions.

In this paper, we study efficient numerical schemes for solving one of the fractional models applied in finance, namely the time-fractional Black-Scholes (B-S) equation ^[4], which is described as follows. For the option price $\mathscr{C}(S, \tilde{t})$ with $S$ being the underlying stock price at the current time $\tilde{t}$, we have the time-fractional B-S equation:

where $r$ is the risk-free rate and $\sigma$ is the volatility of the underlying asset. The fractional operator $\frac{\partial^\alpha}{\partial \tilde{t}^\alpha}$ is defined by

We remark that the model (1.1) recovers the classical B-S equation when $\alpha = 1$.

The time-fractional B-S equation (1.1), which can be derived by assuming that the change in the option price with time is a fractal transmission system, has been extensively studied in recent years ^[7]. One can see that the time-fractional B-S Eq (1.1) becomes degenerate at the underlying asset price $S = 0$, which will lead to the failure of applying classical numerical methods. There are some numerical studies on this topic, see ^[24] and references therein for the corresponding treatment and discussion. So, instead of directly solving the original form of Eq (1.1), we first consider an equivalent form by using the following transformation technique and then solve the resulting equation in a truncated space interval.

Letting $x = \ln{S}$ and $t = T - \tilde{t}$ in model (1.1), and denoting $\mathscr{V}(x, t) = \mathscr{C}(e^x, T-t)$, we get the following equivalent form of the time-fractional B-S equation after simple calculation:

with the boundary and initial value conditions given by $\mathscr{V}(-\infty, t) = \varphi(t)$, $\mathscr{V}(+\infty, t) = \nu(t)$ and $\mathscr{V}(x, 0) = h(x)$, respectively. Without loss of generality, we have added the source term $f$ to the above Eq (1.2). The Caputo derivative ${}_CD_{0, t}^\alpha \mathscr{V}(x, t)$ is defined by

For more details of the above procedure, the readers are referred to our recent paper ^[7].

Generally speaking, finding an analytical solution to the fractional differential equation is hard or even impossible. So, one needs to resort to a numerical solution for solving fractional differential equation. Until now, there have been extensive studies on numerical methods for fractional partial differential equations, see ^{[6,9,18,25,27,28]}, just to name a few. Some main numerical studies for the time-fractional B-S Eq (1.1) are listed below. In ^[29], Zhang et al. constructed a discrete implicit numerical scheme with convergent order $O(\tau^{2-\alpha}+h^2)$ for solving Eq (1.2). Here and in what follows, unless otherwise noting, the notations $\tau$ and $h$ denote the temporal and spatial step sizes, respectively. Roul designed a high order numerical approach which was shown to be unconditionally stable and $O(\tau^{2-\alpha} + h^4)$ accurate ^[19]. In ^[10], Kazmi proposed a stable and accurate finite difference scheme for an equivalent integro-differential form of Eq (1.2) and proved that the convergent order is $O(\tau^{2}+h^2)$ by Fourier analysis. Other relevant results can be found in ^[15,23].

The above numerical schemes are obtained under the assumption that the problem solution is sufficiently smooth. Generally speaking, the solution of the time-fractional model always has weak singularity near $t = 0$. Thus, the accuracy and numerical theory of those above mentioned numerical schemes may not hold when solving Eq (1.2) with non-smooth initial conditions. There are some strategies to deal with such issue, e.g., the non-uniform meshes ^[2,11,14,21], the adding of correction terms ^[16,17,26], etc. We will focus on the first one in this paper. Besides, to the best of our knowledge, it seems that the numerical studies with temporal non-uniform meshes for the time-fractional B-S Eq (1.1), specially with high-order accuracy, are still limited.

We note that Abdi et. al. introduced two high-order finite difference schemes with accuracies of $O(\tau^{3-\alpha} + h^6)$ and $O(\tau^{3-\alpha} + h^8)$ based on the L1-2 formula and compact difference schemes ^[1]. However, the accuracy of their numerical schemes will be heavily affected by the insufficient smoothness of the problem solution. This motivates us to derive high-order and stable numerical schemes for solving the time-fractional B-S Eq (1.1) efficiently.

The main purpose of this paper is to develop the efficient finite difference schemes with temporal graded meshes for solving the model (1.1). The contributions are listed below. First, we apply the variable transformation technique to overcome the degeneration of the original Eq (1.1). Then based on a equivalent form, we succeed in designing two high-order and stable finite difference schemes with convergence orders of $O(N^{-2}+h^6)$ and $O(N^{-2}+h^8)$ for solving time-fractional B-S Eq (1.1) by choosing the grading parameter $\gamma = 1/\alpha$. Here, $N$ denotes the total number of temporal grid points. Second, the corresponding stability and error estimates are derived rigorously based on Fourier stability analysis. Third, the numerical examples, including the numerical comparison with the existing methods, are presented to support our numerical finding. It is worth mentioning that the proposed two high-order compact difference schemes can efficiently solve the fractional model (1.2) with weak singularity of solution near $t = 0$, and can achieve higher accuracy by applying fewer grid points. This means that our derived schemes are more competitive compared to existing numerical schemes, since one can obtain a high-order accuracy and reliable numerical solution of the time-fractional B-S Eq (1.1) using less computational time and cost.

The rest of the paper is organized as follows. In Section 2, we show how to obtain the fully discrete high-order compact difference scheme with temporal graded meshes. The stability and error estimates of the two proposed schemes are given by the Fourier method in Sections 3 and 4, respectively. In Section 5, two numerical examples are given to demonstrate the accuracy and adaptability of our proposed schemes. The conclusions are given in Section 6.

2.   The derivation of the schemes

For the purpose of numerical simulation, the time-fractional B-S Eq (1.2) is considered to be solved on a finite domain $(x_l, x_r) = (0, L)$ instead of the whole domain $\mathbb{R}$. In order to get spatial sixth and eighth order accuracy for the schemes, we use variable substitution $\mathscr{V}(x, t) = \mathrm{e}^{\frac{1}{2}\rho x}u(x, t)$ with $\rho = 1-\frac{2r}{\sigma^2}$ to obtain the following equivalent form of Eq (1.2):

where $p = \frac{1}{2}\sigma^2$, $q = r+\frac{1}{8}\rho^2\sigma^2$ and $F(x, t) = \mathrm{e}^{-\frac{1}{2}\rho x}f(x, t)$.

It is well-known that the time-fractional Eq (2.1) is equivalent to the following Volterra integro-differential equation:

where $\mathcal{L}_x u(x, s) = p\frac{\partial^2 u(x, s)}{\partial x^2}-qu(x, s)$.

In order to perform the derivation of the numerical schemes for the above Eq (2.2), here we make the following assumption on the solution as follows. From here on, we denote the capital letter $C$ as a positive constant which is independent of the temporal and spatial stepsizes $\tau$ and $h$ in this paper.

Assumption 2.1. Assume that the solution in Eq (2.2) satisfies $\vert\frac{\partial^n u}{\partial t^n}(x, t)\vert\leq C(1+t^{\alpha-n})$ for $n = 0, 1, 2$.

We remark that a similar assumption is mentioned in some existing literature, e.g., see [22, Theorem 2.1], [13, Assumption 1] or [30, Theorem A.2].

Next, we discretize the Eq (2.2) in space to get spatial sixth and eighth order accuracy.

For a given positive integer $M$, let $\omega_h \equiv\left\{x_m = x_l+m h \mid 0 \leqslant m \leqslant M\right\}$ be a uniform mesh on the finite interval $[x_l, x_r]$, where the spatial step size $h = (x_r-x_l)/M$. By utilizing the Taylor expansion, one can obtain the following sixth and eighth order compact difference approximations for the second-order spatial partial derivative $\frac{\partial^2 u}{\partial x^2}$ in Eq (2.2)(see, e.g. [5, Lemma 1]):

and

Here,

and the central difference operator $\delta_x^2 u(x_m, t) = u(x_{m-1}, t) -2 u(x_{m}, t) + u(x_{m+1}, t)$.

So, applying the two compact difference formulas (2.3) and (2.4) to (2.2), we obtain

where $u_m(s) = u(x_m, s)$, $F_m(s) = F(x_m, s)$, $\mathcal{L}_{k} u_m(s) = (p{\mathcal{H}_k}-q)u_m(s)$ and the local truncation error $R_{m, k}(t) = \frac{1}{\Gamma(\alpha)} \int_0^t(\frac{\partial^2 u(x_m, s)}{\partial x^2} -{\mathcal{H}_k}u(x_m, s))\; \mathrm{d}s$. Here, $k = 1, 2$.

For the temporal discretization, we divide the interval $[0, T]$ into $0 = t_0 < t_1 < \cdots < t_k < t_{k+1} < \cdots < t_N = T$, with graded meshes $t_j = T(j / N)^{\gamma}(j = 0, 1, 2, \ldots, N,$ and $\gamma \geq 1$). Here, $N$ is a positive integer. Denote $\tau_{j+1} = t_{j+1}-t_{j}(0 \leq j \leq N-1)$. For the sake of discussion, we denote

and

with $s\in[t_j, t_{j+1}]$ and $k = 1, 2$.

Using the piecewise linear interpolation $\Pi^1 \varphi_{m, k}(s)$ for the approximation of $\varphi_{m, k}(s)$ on the each subinterval $[{t_j}, t_{j+1}]$, we can obtain from Eq (2.5) that

where $u_m^n = u(x_m, t_n)$, the truncation error $R_{m, k}^n$ will be given by Lemma 4.1 and the weights $w_{j}^{(n)}$ are described as

and

Remark 2.2. Some strategies on how to reduce the rounding error in the calculation of the weights are proposed in ^[30]. One may refer to these strategies to improve the accuracy of the computation in practice.

Replacing $u_m^n$ with the numerical solution $U_m^n$ in Eq (2.6) and omitting the small term $R_{m, k}^n$, we obtain the following fully discrete high-order compact difference schemes

with $k = 1, 2$. Here, $m = 1, 2, \cdots, M-1$ and $n = 1, 2, \cdots, N$. The corresponding initial and boundary conditions are given by $U_j^0 = \mathrm{e}^{-\frac{1}{2} \rho x_j} h(x_j), j = 0, 1, \cdots, M$, and $U_0^n = \mathrm{e}^{-\frac{1}{2} \rho x_l} \varphi(t_n), U_M^n = \mathrm{e}^{-\frac{1}{2} \rho x_r} \nu(t_n), n = 1, 2, \cdots, N$.

By the definitions in Eqs (2.3) and (2.4), we respectively have the sixth-order compact difference scheme:

and the eighth-order compact difference scheme:

where

For ease of implementation in practice, in the following we present the concrete matrix forms of the two compact difference schemes (2.8) and (2.9).

By the definition of $\mathcal{A}_1$, we have

where

and

Denote by

and

then the sixth-order compact difference scheme (2.8) can be expressed in the matrix form:

where the matrices ${\bf D}$ and ${\bf E}$ are all $(M-1) \times(M-1)$ five-diagonal matrices given by

and

The column vector ${\bf G}^n$ with $(M-1)$-entries is given by

Similarly, the eighth-order compact difference scheme (2.9) has the following matrix from.

where the two $(M-1)\times(M-1)$ seven-diagonal matrices $\tilde{{\bf D}}$ and $\tilde{{\bf E}}$ are given by

with

and

The column vector $\tilde{{\bf G}}^n$ with $(M-1)$-entries is defined by

We remark that the ghost-point values in the sixth-order and eighth-order compact difference schemes (2.8) and (2.9) can be computed by the following extrapolation formulas ^[5].

for the compact difference scheme (2.8), and

for the compact difference scheme (2.9).

3.   Stability analysis

We now present the stability of the two high-order compact difference schemes (2.8) and (2.9) by Fourier method (cf., e.g. ^[3]), respectively.

3.1.   Stability analysis of the scheme (2.8)

Let $\hat{U}_m^n$ be the solution of Eq (2.8). Define $\varepsilon_m^n = U_m^n-\hat{U}_m^n, m = 1, \cdots, M-1, n\geq 0$. It follows from Eq (2.8) that the round off error equation has the following form:

where $m = 1, \cdots, M-1$ and $n \geq 1$.

Define the following grid function:

from which we have the following Fourier series expansion of $\varepsilon^n(x)$:

where the Fourier coefficient $\xi_n(k)$ is given by

Denote $\varepsilon^n = (\varepsilon_1^n, \varepsilon_2^n, \cdots, \varepsilon_{M-1}^n)^T$. Using the definition of the discrete $L^2$-norm:

and Parseval's identity:

we deduce that

So, letting $\varepsilon_m^n = \xi_n \mathrm{e}^{i\ell m h}$ with $\ell = 2\pi k/L$ and substituting it into the Eq (3.1), we have

where we have used the two identifies $\cos 2\ell h = \frac{\mathrm{e}^{2i \ell h}+\mathrm{e}^{-2i \ell h}}{2 }$ and $\cos \ell h = \frac{\mathrm{e}^{i \ell h}+\mathrm{e}^{-i\ell h}}{2}$. It follows that

where the coefficients $\mathcal{B}_j(j = 1, 2)$ are defined by

Lemma 3.1. For the coefficients $\mathcal{B}_j(j = 1, 2)$ given in Eq (3.3), one has the following two inequalities:

where $K_1 = \left|\frac{15(3q+4\ell^2 p)}{37}\right|$.

Proof. The first inequality is obvious in view of $p > 0, q > 0, h > 0$ and $0 < \alpha < 1$. For the second inequality, we have

where the last inequality holds for sufficiently small $h$. Thus, the proof is completed. □

Lemma 3.2. Suppose that $\xi_n(n\geq 1)$ is the solution of Eq (3.2), we have

Proof. First, applying Lemma 3.1 to Eq (3.2), one has

Noting that the property of weights $w_j^{(n)}\leq C \tau_{j+1}(t_n-t_j)^{\alpha-1}$ holds (see [11, Lemma 3.1]). So, the desired result then follows from the modified Gronwall inequality (see [11, Lemma 3.3]). □

We are now in a position to give the stability of the sixth-order compact difference scheme (2.8).

Theorem 3.3. The sixth-order compact difference scheme (2.8) is stable.

Proof. In view of the estimate result (3.4), we obtain

which implies that the scheme (2.8) is stable. □

3.2.   Stability analysis of the scheme (2.9)

Similar to the discussion in the previous subsection for Eq (2.8), we assume that $\tilde{U_m^n}$ is the approximate solution of $U_m^n$ and define $\tilde{\varepsilon_m^n} = U_m^n - \tilde{U_m^n}$, $m = 0, 1, \cdots, M, n = 0, 1, \cdots, N$. So, we can get the roundoff error equation for Eq (2.9):

where $m = 1, \cdots M-1$ and $n \geq 1$.

Here, we similarly define $\tilde{\varepsilon}_m^n = \tilde{\xi}_n\mathrm{e}^{i\ell mh}$ with $\ell = 2\pi k/L$. Substituting it into the above equation, one has

from which we obtain

where

Since $\mathcal{C}_1+\mu^{(\alpha)}\mathcal{C}_2\neq 0$, we divide the Eq (3.5) by this coefficient to get

One can obtain the following lemmas in a similar way of the proof presented in the stability of the sixth-order compact difference scheme (2.8). Thus the proofs are omitted.

Lemma 3.4. The following estimates hold.

where $K_2 = \left| \frac{7(45q+68\ell^2p)}{279}\right|$.

Lemma 3.5. Assume that $\tilde{\xi}_n(n\geq 1)$ is the solution of Eq (3.6), we have

Theorem 3.6. The eighth-order compact difference scheme (2.9) is stable.

Proof. The proof is analogous to the proof of Theorem 3.3, so the details are omitted for the sake of brevity. □

4.   Convergence analysis

Lemma 4.1. Under Assumption (2.1), the truncation error in the high-order compact difference schemes (2.7) has the following form:

where $\dagger = 6$ for $k = 1$ and $\dagger = 8$ for $k = 2$.

Proof. Applying the error estimates presented in [5, Lemma 1] and [13, Lemma 2.1] to Eqs (2.5) and (2.6), respectively, one can readily obtain the desired result. The proof is thereby finished. □

From here on, unless otherwise noting, we always choose the grading parameter $\gamma = 1/\alpha$ in order to obtain the highest possible rate of convergence $O(N^{-2}+h^\dagger)$.

4.1.   Convergence analysis of the scheme (2.8)

Denote the error $e_m^n = u(x_m, t_n)-U_m^n$. Here, $U_m^n$ is the numerical solution of the sixth-order compact difference scheme (2.8). By Eqs (2.6) and (2.8), we have the following error equation:

where $m = 1, 2, \cdots, M-1$, and $n = 1, 2, \cdots, N$.

Similar to the stability analysis in the previous section, we define the grid functions as follow:

and

So, we expand $e^n(x)$ and ${R}_1^n(x)$ into the following Fourier series:

and

The Fourier coefficients $\zeta_n(k)$ and $\eta_n(k)$ are given by

and

In a similar way, we have

and

So, we define $e_m^n = \zeta_n \mathrm{e}^{i \ell mh}$ and $R_{m, 1}^n = \eta_n \mathrm{e}^{i \ell mh}$ with $\ell = 2\pi k/L$. It follows from Eqs (4.1) that

where the definitions of $\mathcal{B}_j(j = 1, 2)$ are given by Eq (3.3).

Next, we show that the solution $\zeta_n$ of Eq (4.3) is bounded.

Lemma 4.2. Suppose that $\zeta_n(n\geq 1)$ is the solution of Eq (4.3), we have

Proof. Using Lemma 4.1, we have

It follows that the term $\sum_{k = -\infty}^{\infty}\vert\eta_n(k) \vert^2$ in Eq (4.2) converges. We conclude that there exists a positive constant $C$ such that

for $n\geq 1$.

So, by Eq (4.3), one has

Utilizing Lemma 3.1 and the property of the cofficient $|\mathcal{B}_1+\mu^{(\alpha)} \mathcal{B}_2|\neq 0$, one gets

Since $w_j^{(n)}\leq C \tau_{j+1}(t_n-t_j)^{\alpha-1}$ (see [11, Lemma 3.1]), applying the modified Gronwall inequality again (see [11, Lemma 3.3]), we can obtain the desired estimate result. So, the proof is completed. □

Now, we present the error estimate for the sixth-order compact difference scheme (2.8).

Theorem 4.3. Let $U_m^n$ and $u(x_m, t_n)$ be solutions of the sixth-order compact difference scheme (2.8) and the model problem (2.1), respectively. Under Assumption 2.1, we have the following error estimate:

Proof. In view of Lemma 4.4, we have

Thus, we complete the proof. □

4.2.   Convergence analysis of the scheme (2.9)

Denote the error $\tilde{e}_m^n = u(x_m, t_n)-U_m^n$. Here, $U_m^n$ is the numerical solution of the eighth-order compact difference scheme (2.9). By Eqs (2.6) and (2.9), we have the following error equation:

where $m = 1, 2, \cdots, M-1$, and $n = 1, 2, \cdots, N$.

Similarly, we define $\tilde{e}_m^n = \tilde{\zeta}_n \mathrm{e}^{i \ell mh}$ and $R_{m, 2}^n = \tilde{\eta}_n \mathrm{e}^{i \ell mh}$ with $\ell = 2\pi k/L$. So, we have

where $\mathcal{C}_j(j = 1, 2)$ is given by Eq (3.6).

Analogous to the proofs of error estimate in the sixth-order compact difference scheme (2.8), we only list here the lemma and theorem and without providing the corresponding proof details for the sake of brevity.

Lemma 4.4. Suppose that $\tilde{\zeta}_n(n\geq 1)$ is the solution of Eq (4.4), we have

Theorem 4.5. Let $U_m^n$ and $u(x_m, t_n)$ be the solutions of eighth-order compact difference scheme (2.9) and the model problem (2.1), respectively. Under Assumption 2.1, we have the following error estimate:

5.   Numerical examples

In this part, we provide numerical examples to verify the accuracy of the two high-order compact difference schemes (2.8) and (2.9). Besides, the numerical simulation for the time-fractional B-S model which describes the European options price are performed to illustrate the practicability of our proposed schemes.

Example 5.1. Consider problem (2.1) with $r = 0.5$, $\sigma = \sqrt{2}$ and $T = 1$. The initial and boundary conditions are given by $u(x, 0) = x^3+x^2+1$ and $u(0, t) = t^{\alpha}+1, u(1, t) = 3(t^{\alpha}+1)$. The source term $f$ is constructed such that the exact solution is $u(x, t) = (t^\alpha + 1)(x^3+x^2+1)$.

Since the fractional order $\alpha\in(0, 1)$, the problem solution is not sufficiently smooth, especially near $t = 0$. The errors of the numerical solution in $L^2$-norm are defined by $E(M, N) = \max_{1\leq j\leq N}\|U^j-u\left(\cdot, t_j\right)\|_2$. The temporal and spatial convergence orders are computed by $\log_{\frac{N_2}{N_1}}\left(\frac{E(M, N_1)}{E(M, N_2)}\right)$ and $\log_{\frac{M_2}{M_1}}\left(\frac{E(M_1, N)}{E(M_2, N)}\right),$ respectively. Here, $N_1 < N_2$ and $M_1 < M_2$. Numerical results are shown in Tables 1 and 2.

Notice that we have set $M = \lceil\sqrt[3]{N}\rceil$ in order to observe the temporal order more clearly when testing the temporal convergence order, see Table 1. So, according to the convergence results, i.e., Theorems 4.3 and 4.5, the theoretical temporal convergence order would be no less than two. From the numerical results, we indeed observe that the convergence orders of proposed schemes (2.8) and (2.9) are consistent with the theoretical accuracy, which are almost $O(N^{-2} + h^6)$ and $O(N^{-2} + h^8)$, respectively. This demonstrates that the proposed schemes (2.8) and (2.9) are robust and accurate when solving the non-smooth solution problem.

Example 5.2. We consider the following European put option problem:

with the parameters: $\sigma = 0.25, r = 0.05, S_l = 0.1, S_r = 100$ and $K = 50$. The final time is $T = 1$(year).

In view of the relationship between the Eqs (1.1) and (2.1), we now test the accuracy of our derived numerical schemes (2.8) and (2.9) for practical use. For sake of discussion, we only focus on the temporal convergence order. The convergence order in time is computed by $\log_2\Big(\frac{E(M, N)}{E(M, N/2)}\Big)$ with the errors $E(M, N) = \|U^N - U^{N/2}\|_2$ and the numerical solution $U^N(x_j) \approx \mathscr{C}(\mathrm{e}^{x_j}, T-t_N)$. In practice, in order for the spatial errors to not contaminate the temporal error, we will make the numbers of the spatial nodes $M$ large enough. The numerical results obtained by the two numerical schemes (2.8) and (2.9), and the two ones (11) and (13) in paper ^[1], are shown in Table 3.

The numerical results presented in Table 3 indicate that the numerical schemes (2.8) and (2.9) can maintain the second-order accuracy in time, while the two ones (11) and (13) in paper ^[1] only reach first-order accuracy. This implies that our schemes give better convergence results in time when numerically solving the time-fractional B-S Eq (1.1).

Finally, we perform the numerical simulation for the time-fractional B-S model (1.1).

Example 5.3. Consider the numerical simulations of the European put option and the European call option. The corresponding initial and boundary conditions are listed as follows.

● European put option: $\mathscr{C}(S, T) = \max \{K-S, 0\}, \mathscr{C}(S_l, \tilde{t}) = K\mathrm{e}^{-r(T-\tilde{t})}$, and $\mathscr{C}(S_r, \tilde{t}) = 0$.

● European call option: $\mathscr{C}(S, T) = \max\{S-K, 0\}, \mathscr{C}(S_l, \tilde{t}) = 0$, and $\mathscr{C}(S_r, \tilde{t}) = S_r-K\mathrm{e}^{-r(T-\tilde{t})}$;

Here, the parameters are chosen as $\sigma = 0.25, r = 0.05, S_l = 0.1, S_r = 100$, and $K = 50$, which is the same as in Eq (5.1). The final times are both chosen as $T = 1$(year). Notice that the numerical schemes (2.8) and (2.9) still hold for solving the classical B-S model, i.e., the case $\alpha = 1$ in time-fractional B-S model (1.1), so by applying the sixth-order compact difference scheme (2.8) for $\alpha = 0.3, 0.8,$ and $1$, we obtain the corresponding numerical simulations and present them in Figures 1 and 2.

From Figures 1 and 2, we can see that, close to the strike price ($K = 50$), the time-fractional derivatives have a significant effect on the option price (see the zoomed parts in Figures 1 and 2) which further verifies the conclusion of the numerical simulation in ^[20]. It seems that, compared with the classical B-S model, the time-fractional B-S model has richer dynamical behaviour in the characteristics of complex movements that appeared in the financial market, and the dynamical mechanism behind it needs to be further studied.

6.   Conclusions

In this paper, we construct two high-order compact finite difference schemes based on temporal graded meshes for solving the time-fractional B-S equation. The stability and convergence of the proposed schemes on the temporal graded meshes are proved by Fourier method. Finally, the numerical examples demonstrate that the two proposed schemes are robust with the high-order accuracy (namely, $O(N^{-2}+h^6)$ and $O(N^{-2}+h^8)$, by choosing the grading parameter $\gamma = 1/\alpha$), even if the problem solution has weak singularity near $t = 0$, which are in a good agreement with the theoretical results.

For the two high-order compact finite difference schemes proposed in this paper, one possible way to further improve their computational efficiency is by employing the idea of sum-of-exponentials (see, e.g., ^[8]) for the integral term appearing in the equivalent form (2.2) of the time-fractional B-S equation. This will be one of our upcoming works.

Use of AI tools declaration

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

Acknowledgments

This work was supported by the Guangxi Natural Science Foundation under grant numbers 2021GXNSFBA196027 and 2023GXNSFAA026315.

Conflict of interest

The authors declare there is no conflict of interest.