A High-Order Numerical Method Based on a Spatial Compact Exponential Scheme for Solving the Time-Fractional Black–Scholes Model

: This paper investigates a high-order numerical method based on a spatial compact exponential scheme for solving the time-fractional Black–Scholes model. Firstly, the original time-fractional Black–Scholes model is converted into an equivalent time-fractional advection–diffusion reaction model by means of a variable transformation technique. Secondly, a novel high-order numerical method is constructed with ( 2 − α ) accuracy in time and fourth-order accuracy in space based on a spatial compact exponential scheme, where α is the fractional derivative. The uniqueness of solvability of the derived numerical method is rigorously discussed. Thirdly, the unconditional stability and convergence of the derived numerical method are rigorously analyzed using the Fourier analysis technique. Finally, numerical examples are presented to test the effectiveness of the derived numerical method. The proposed numerical method is also applied to solve the time-fractional Black–Scholes model, whose exact analytical solution is unknown; numerical results are illustrated graphically.


Introduction
Due to the nonlocal nature of the fractional differential and integral operators, fractional calculus has gained more and more attention from researchers.Fractional calculus has provided a powerful tool for modeling anomalous behaviours from various fields of science and engineering, such as material science, biology, physics, control science, fluid mechanics and finance [1][2][3].
In the financial market, the famous Black-Scholes model [4,5] was derived to determine the option price under some strict assumptions.However, the model still has some drawbacks; for example, it cannot capture some significant movements or jumps within a small time step [6].Hence, the Black-Scholes model does not often show the actual performance of stock prices in a realistic financial market.Since the fractional differential and integral operators make it easy to characterize the behaviors with memory and hereditary properties, the fractional Black-Scholes model was generalized [7,8].Liang et al. [9] established a novel two-parameter time-fractional Black-Scholes-Merton model for option pricing.Chen et al. [10] considered the variation in the option price as a fractal transmission system and derived a time-fractional Black-Scholes model.Zhang et al. [11] gave a comprehensive review for various fractional Black-Scholes equations.
Since the exact solutions for various kinds of fractional partial differential equations are not often available, numerical methods for solving these equations must be considered.Various efficient numerical methods have already been investigated [12][13][14][15][16][17], such as finite difference methods, finite volume methods, finite element methods and spectral methods.For the time-fractional Black-Scholes model, Zhang et al. [6] constructed a numerical scheme with (2 − α) accuracy in time and second-order spatial accuracy, and analyzed the stability and convergence.Soon, De Staelen and Hendy [18] constructed a high-order difference method with fourth-order spatial accuracy, and analyzed the stability and convergence.Similar results can be found in [19].Roul [20] constructed a numerical method based on a B-spline collocation technique.An et al. [21] presented a space-time spectral method and analyzed the stability and convergence of the numerical method.Song and Lyu [22] investigated a high-order and fast numerical method for the time-fractional Black-Scholes equation, with a weak initial singularity of the solution.Abdi et al. [23] introduced two compact finite difference methods.Taghipour and Aminikhah [24] proposed an efficient spectral collocation method.Kaur and Natesan [25] proposed a numerical method based on the cubic spline method.Zhang and Zheng [26] proposed the finite element method for solving the variable-order time-fractional Black-Scholes model.Kazmi [27] developed an efficient numerical method with second-order temporal and spatial accuracy.
In this paper, we investigate the following time-fractional Black-Scholes model [10]: where V(S, τ) represents the price of an option, S is the asset price, τ is the current time, σ is the volatility of the underlying asset, r is the risk-free interest rate, T is the expiry time, the function ψ(S) is the terminal condition which denotes the terminal payoff function of the option, and the functions p(τ) and q(τ) are boundary conditions which denote the rebates paid when the corresponding barrier is hit [10]; more details can be found in [10].The time-fractional derivative in model ( 1) is the modified right Riemann-Liouville time-fractional derivative When the fractional derivative α = 1, the above mentioned model (1) becomes the classical Black-Scholes model.While 0 < α < 1, the above mentioned model (1) represents a time-fractional Black-Scholes model [10].An example describing the European double barrier knock-out call option can be found in Section 4. Letting S = e x and τ = T − t, and denoting u(x, t) = V(e x , T − t), the above mentioned time-fractional Black-Scholes model (1) can be rewritten as and the time-fractional derivative is transformed into the following form: Actually, the fractional derivative 0 D α t u(x, t) in model (3) represents the Caputo timefractional derivative C 0 D α t u(x, t) which is already derived in reference [6,10] 0 In order to obtain the numerical solution of the above time-fractional Black-Scholes model (3), it is necessary to truncate the original unbounded domain (0, ∞) × (0, T) into a finite domain (x L , x R ) × (0, T).Hence, the aim of this paper is to investigate the following time-fractional advection-diffusion reaction model: as a source term is added for numerical experimental validation.
The outline of the paper is arranged as follows.In Section 2, a high-order numerical method is derived with (2 − α) temporal accuracy and fourth-order spatial accuracy, and the uniqueness of solvability is rigorously discussed.In Section 3, the stability and the convergence of the derived high-order numerical method is strictly analyzed using the Fourier analysis technique.In Section 4, numerical examples are given to test the effectiveness of the proposed numerical method.The proposed numerical method is also applied to solve the time-fractional Black-Scholes model, whose exact analytical solution is unknown; the call option prices for different fractional derivative α are illustrated graphically.Finally, some conclusions are summarized.

Construction of the High-Order Numerical Method
In this section, we will derive a high-order numerical method based on a spatial compact exponential scheme for the derived fractional advection-diffusion reaction model (6).
Firstly, let M and N be two positive integers, h = x R −x L M the space step, and ∆t = T N the time step.Denote Ω h = {x j |x j = x L + jh, j = 0, 1, 2, . . ., M} and Ω τ = {t n |t n = n∆t, n = 0, 1, 2, . . ., N}. Denote the grid function u n j = u(x j , t n ) and f n j = f (x j , t n ).Introduce the following notations: In order to derive a high-order numerical discrete scheme in space, we first investigate the corresponding convection-diffusion equation where g(x, t) is a suitable smooth function.A fourth-order compact exponential scheme has already been derived [28,29] Bu j = Ag j + O(h 4 ), where with Substituting g(x, t) = C 0 D α t u(x, t) + cu(x, t) − f (x, t) in the Equation ( 9) at the mesh point (x j , t n ), we can obtain the following semi-discrete high-order compact exponential approximation scheme: Next, in order to construct a fully-discrete numerical approximation method, we introduce the following two lemmas: where α ∈ (0, 1) and Lemma 2. The coefficients a n (n = 0, 1, 2, 3, . ..) in ( 13) satisfy From Lemma 1, it is known that the Caputo time-fractional derivative Substituting ( 15) into ( 12) and denoting µ = ∆t −α Γ(2−α) , we can obtain that, for n = 1, and, for n ≥ 2, where Substituting u n j with its numerical approximations U n j and omitting R n j , we have derived a high-order numerical scheme as follows: and, for n ≥ 2, Finally, we discretize the initial condition and the boundary conditions by Next, we will discuss the unique solvability of the derived numerical method ( 19)-( 22).
Proof.The matrix form of the derived numerical method ( 19)-( 22) can be written as where with Noting that and letting z = bh 2a , we can obtain Therefore, Using Lemma 3 in [29], we can obtain Therefore, which means the coefficient matrix P is a strictly diagonal dominant matrix, hence, matrix P is non-singular.Therefore, the solution is unique.

Stability and Convergence Analysis
For the derived high-order numerical method ( 19)-( 22), the stability and convergence analysis will be rigorously discussed by the Fourier analysis method.

Stability Analysis
For simplicity, we first denote and r 1 , r 2 , r 3 are already denoted in ( 25)- (27).Suppose that Ũn j is the approximation solution of the above numerical scheme ( 19)- (22).Let with ρ n 0 = ρ n M = 0, and T .Therefore, we can obtain that, for n = 1, and, for n ≥ 2, we have Define the grid function and periodically expand ρ n (x) in the Fourier series [6] with the period where Using the Parseval equality and the following discrete norm we have Based on the above analysis, we assume that the solution of the above error Equations (31) and (32) has the expression where σ = 2πl/L.
Lemma 5.The following inequality is established: Proof.First, it is known from Lemma 4 that, for n = 1, Assuming that then hence, the proof is completed by means of mathematical induction.
Theorem 2. The derived high-order numerical method ( 19)-( 22) is unconditionally stable with respect to the initial value, that is, Proof.From Lemma 5, we have hence, the proof is completed.19) and (20) from Equations ( 16) and( 17), we can obtain the error equations

Convergence Analysis
and, for n ≥ 2, where Define the following grid functions: and Then e n (x) and R n (x) can be expressed in the following form: where Similarly as in stability analysis, we have Due to the convergence of the series on the right hand side of the Equation ( 54), there exists a positive constant C 2 > 0, such that Based on the above analysis, we assume that the solutions of the above error Equations (49) and (50) have the expressions where σ = 2πl/L.Substituting the expressions (56) into ( 49) and (50), we have and, for n ≥ 2, Lemma 6.For 0 < ∆t < 1, we have Proof.Since µ = ∆t −α Γ(2−α) , ∆t ∈ (0, 1), obviously we can obtain that is, Hence, the proof is complete.

Lemma 7.
The following inequality is established: where C 2 is denoted in (55).
Proof.It is known from Lemma 6 that, for n = 1, hence, the proof is finished by means of mathematical induction.
Theorem 3. The derived high-order numerical method ( 19)-( 22) is unconditionally convergent with convergence order O(∆t 2−α + h 4 ), that is, there exists a positive constant C > 0, such that Proof.By means of ( 18), we have From Lemma 7, we have where therefore, the proof is finished.

Numerical Examples
In this section, we present some numerical examples to verify the effectiveness and accuracy of the derived high-order numerical method.
Define Error The convergence order of the derived high-order numerical method in time is calculated by and the convergence order of the derived high-order numerical method in space is calculated by Order = log 2 e(∆t, 2h) e(∆t, h) .
We solve the problem using the derived high-order numerical method ( 19)- (22).Table 1 shows the calculated errors and the corresponding convergence orders in time for different fractional derivative α = 0.3, 0.5, 0.7 and a fixed h = 1/1000.From Table 1, we can obviously see that the convergence order in time is O ∆t 2−α , which is in agreement with the theoretical analysis.Table 2 shows the calculated errors and the corresponding convergence orders in space for different fractional derivative α = 0.3, 0.5, 0.7 and a fixed ∆t = 1/100,000.From Table 2, one can see that the convergence order in space is O h 4 , which is in agreement with the theoretical analysis.Table 3 shows the comparison of error and convergence order in space, calculated using the derived high-order numerical method ( 19)-( 22) and previous research [6,18,19] when the fractional derivative α = 0.7 and ∆t = 1/100,000.From Table 3, one can see that the derived high-order numerical method ( 19)-( 22) gains the same convergence order in space as [18,19], which is higher than [6].Furthermore, all the above analyses demonstrate that our derived high-order numerical method ( 19)-( 22) is effective.Figure 1 illustrates the comparison between the exact solutions and the numerical solutions at t = 1/3, 1/2, 2/3 when α = 0.2. Figure 2 illustrates the comparison between the exact solutions and the numerical solutions at x = 1/4, 1/2, 3/4 when α = 0.9.From Figures 1 and 2, one can obviously observe that the numerical solutions are remarkably consistent with the exact solutions, which indicates the efficiency of the derived high-order numerical method.Example 2 (unknown solution).Consider the following time-fractional Black-Scholes model, which describes the European double barrier knock-out call option [6,10]: Here, the parameter values are taken to be the same as in reference [6,10] : σ = 0.45, r = 0.03, D = 0.01, T = 1, K = 10 .The double barrier knock-out call is obtained if the price of a knock-out call whose strike is within the two barriers is worked out; more details can be found in [10].Example 2 (unknown solution).Consider the following time-fractional Black-Scholes model, which describes the European double barrier knock-out call option [6,10]: Here, the parameter values are taken to be the same as in reference [6,10] : σ = 0.45, r = 0.03, D = 0.01, T = 1, K = 10 .The double barrier knock-out call is obtained if the price of a knock-out call whose strike is within the two barriers is worked out; more details can be found in [10].Example 2 (unknown solution).Consider the following time-fractional Black-Scholes model, which describes the European double barrier knock-out call option [6,10]: Here, the parameter values are taken to be the same as in reference [6,10]: σ = 0.45, r = 0.03, D = 0.01, T = 1, K = 10.The double barrier knock-out call is obtained if the price of a knock-out call whose strike is within the two barriers is worked out; more details can be found in [10].
Since the exact analytical solution is not known, we solve this problem by the derived high-order numerical method ( 19)-( 22) for different fractional derivative value.Numerical results are presented graphically.
Figure 3 illustrates the European double barrier knock-out call option prices for different time-fractional derivative α = 0.1, 0.3, 0.5, 0.7, 0.9, 1.From Figure 3, one can see that how the call option prices are affected by the strike price K and the time-fractional derivative α.Furthermore, we can conclude that the time-fractional Black-Scholes model can display more characteristics of the call option price changes than the classical Black-Scholes model (α = 1).Since the exact analytical solution is not known, we solve this problem by the derived high-order numerical method ( 19)-( 22) for different fractional derivative value.Numerical results are presented graphically.
Figure 3 illustrates the European double barrier knock-out call option prices for different time-fractional derivative α = 0.1, 0.3, 0.5, 0.7, 0.9, 1.From Figure 3, one can see that how the call option prices are affected by the strike price K and the time-fractional derivative α.Furthermore, we can conclude that the time-fractional Black-Scholes model can display more characteristics of the call option price changes than the classical Black-Scholes model (α = 1).

Conclusions
Due to the fact that the exact solution of the fractional partial differential equation is not easy to obtain, and can even be impossible, we have derived a new high-order numerical method based on a spatial compact exponential scheme for solving the time-fractional Black-Scholes model.Firstly, by means of the variable transformation technique, we convert the original time-fractional Black-Scholes equation into an equivalent time-fractional advectiondiffusion reaction equation, then we construct a new high-order numerical method with (2 − α) accuracy in time and fourth-order accuracy in space based on spatial compact exponential scheme, which is different from the previous research [6,18,19].We have rigorously discussed the uniqueness of solvability, and the unconditional stability and convergence of the derived high-order numerical method.Finally, two examples, whose exact solutions are known and unknown, are given to illustrate the efficiency of the derived high-order numerical method.

Conclusions
Due to the fact that the exact solution of the fractional partial differential equation is not easy to obtain, and can even be impossible, we have derived a new high-order numerical method based on a spatial compact exponential scheme for solving the time-fractional Black-Scholes model.Firstly, by means of the variable transformation technique, we convert the original time-fractional Black-Scholes equation into an equivalent time-fractional advectiondiffusion reaction equation, then we construct a new high-order numerical method with (2 − α) accuracy in time and fourth-order accuracy in space based on spatial compact exponential scheme, which is different from the previous research [6,18,19].We have rigorously discussed the uniqueness of solvability, and the unconditional stability and convergence of the derived high-order numerical method.Finally, two examples, whose exact solutions are known and unknown, are given to illustrate the efficiency of the derived high-order numerical method.

Figure 3 .
Figure 3.European double barrier knock-out call option prices for different fractional derivative α.

Author Contributions:
Conceptualization, B.Y.; methodology, X.H. and B.Y.; software, X.H. and B.Y.; validation, X.H. and B.Y.; formal analysis, X.H. and B.Y.; investigation, X.H. and B.Y.; writingoriginal draft preparation, X.H. and B.Y.; writing-review and editing, X.H. and B.Y.; supervision, B.Y.; funding acquisition, B.Y.All authors have read and agreed to the published version of the manuscript.Funding: This research was funded by Shandong Provincial Natural Science Foundation, China, grant number ZR2022MA021 and China Postdoctoral Science Foundation, Grant number 2018T110679.Data Availability Statement: Data are contained within the article.Acknowledgments:The authors thank the editors and the referees for their valuable comments and helpful suggestions.

Figure 3 .
Figure 3.European double barrier knock-out call option prices for different fractional derivative α.

Table 1 .
Error and convergence order in time for different fractional derivative α.

Table 2 .
Error and convergence order in space for different fractional derivative α.