COMPACT FINITE DIFFERENCE SCHEMES OF THE TIME FRACTIONAL BLACK-SCHOLES MODEL∗

In this paper, three compact difference schemes for the time-fractional Black-Scholes model governing European option pricing are presented. Firstly, in order to obtain the fourth-order accuracy in space by applying the Padé approximation, we eliminate the convection term of the B-S equation by an exponential transformation. Then the time fractional derivative is approximated by L1 formula, L2− 1σ formula and L1− 2 formula respectively, and three compact difference schemes with oders O(∆t2−α +h), O(∆t +h) and O(∆t3−α + h) are constructed. Finally, numerical example is carried out to verify the accuracy and effectiveness of proposed methods, and the comparisons of various schemes are given. The paper also provides numerical studies including the effect of fractional orders and the effect of different parameters on option price in time-fractional B-S model.


Introduction
Black and Scholes [2] proposed the famous option pricing formula in 1973 , the classical Black-Schels model showed that the stock price changes must follow the assumption of geometric Brownian motion. The Black-Scholes model has been increasingly popular because it effectively models the option value and provides a mechanism for extracting implied volatilities. However, one of the biggest drawbacks of the classical B-S model is that it can't capture large movements or jumps over small time steps in the dynamic process of stock price changes. As early as the 1960s, Mandelbrot [15] observed the long-tailed distribution of relative stock price changes, and deduced that the use of α−stable Lévy motion was instead of standard geometric Brownian motion. Because fractional derivatives can describe the characteristics of memory and inheritance and are very close to Lévy processes, fractional differential equations have become powerful tools for studying fractal geometry and fractal dynamics, and have been widely used to model anomalous diffusion or α−stable Lévy processes.
There are mainly two types of B-S models which follow a fractal transmission system: a spatial-fractional B-S model and a time-fractional B-S model. As for a fractional derivative in space, Carr and Wu [3]introduced the FMLS (finite moment log stable) model. Cartea et al. [4] showed that the price of European options satisfies the FPED (fractional partial differential equation) with spatial fractional derivative under FMLS model. Chen et al. [5] derived an explicit closed-form analytical solution for the equation proposed by [4]. In [22] Zhang et al. proposed an implicit discrete scheme for the tempered fractional B-S model governing a European double-knock-out barrier option. Regarding the option pricing with respect to time-fractional derivative, Wyss [19] gave a time-fractional B-S equation for a European vanilla options. Cartea [6] found the value of European-style derivatives satisfies a FPDE with the Caputo time-fractional derivative in modeling stock price by using tick-to-tick data. Applying fractional Taylor formula, Jumarie [11] deduced a time and space fractional B-S equations. Liang et al. [13] established a model for option pricing of two parameters-fractional Black-Scholes-Merton differential equation and obtained explicit option pricing formulas for European options. Chen et al. [7] simplified the model of [13].
In this paper we consider the following time-fractional Black-Scholes model [7] where U (S, τ ) denotes the price of an option with S being the asset price and τ being the current time, r > 0 is the risk-free interest rate, σ > 0 is the volatility of underlying asset, T > 0 is the expiry time. ∂ α U (S,τ ) ∂τ α (0 < α ≤ 1) is a modified right Riemann-Liouville fractional derivative defined as When α = 1 the model (1.1) reduces the classical B-S model. With a growing number of time-fractional B-S models being proposed, different approaches to solve the fractional models have also been developed. Song et al. [16] and Zhang et al. [23] gave different implicit difference schemes for time-fractional B-S equation of European put options, and their convergence rates are respectively O(∆t + h 2 ) and O(∆t 2−α + h 2 ) (Among them, ∆t denotes the temporal step size, h denotes the spatial step size, 0 < α < 1 is the order of fractional derivative, the same as below). Based on the work of [23], Staelen and Hendy [17] constructed a fourth order implicit difference scheme in space and 2 − α order in time. Due to the slow calculation speed of implicit difference method, Yang et al. proposed the Explicit-Implicit scheme and Implicit-Explicit scheme with convergence rate [20] and pure alternative segment explicit-implicit parallel difference scheme with convergence rate O(∆t 2−α + h 2 ) in [24], respectively. Chen et al. [8] proposed the first-order upwind finite difference scheme with a uniform mesh for American puts under the generalized mixed fractional Brownian motion (GMFBM) model. Zhou and Gao [25] developed a Laplace transform method and a boundary-searching finite difference method for a free-boundary time-fractional B-S equation of American option pricing problem with convergence rate O(∆t 2−α + h 2 ). Koleva and Vulkov [12] presented a weighted finite difference method for a timefractional B-S equation with convergenc rate O(∆t + h 2 ). In [9], Cen et al. made an integral discretization scheme in time coordinate direction and employed a central difference scheme for the spatial discretization for a time-fractional B-S equation. Numerical experiments show that their proposed scheme is more accurate and robust when α is close to 0. The convergence rate of the method is O(∆t + h 2 ). By employing the universal difference method, Yang et al. [21] solved the time-space fractional B-S model with the boundary conditions satisfied by standard European call options, and its convergence order is O(∆t + h 2 ).
From the previous work, we can see that there are many low order numerical schemes to solve the time-fractional Black-Scholes model. In this paper, we will propose three compact difference schemes for time-fractional B-S model governing European options to improve numerical accuracy.
The remainder of the paper is organized as follows. In Section 2, eliminating the convection term of the B-S equation by an exponential transformation, the original equation is transformed into an equivalent form. In Section 3, three compact finite difference schemes are presented for solving the time-fractional B-S model. In Section 4, numerical examples are carried out to verify the high accuracy and efficiency of our methods. A conclusion is given in Section 5.

Time fractional Black-Scholes model and its equivalent model
In order to eliminate the variable coefficient S in the model (1.1), we introduce the following transformations: Note that Zhang et al. [23] the modified right R-L fractional derivative can be transformed the following Caputo' form: Model (1.1) can be rewritten as (2.1) For solving the above model numerically, we need to truncate the original unbounded domain into a finite interval, and add a source term f (x, t) to the righthand side of the equation without loss of generality. The model (2.1) can be formulated as the following form: To obtain the forth-order accuracy in space of the scheme by utilizing the Padé approximation, we first multiply 2 σ 2 on the both sides of the equation in model (2.2) as: Then let 1 − 2r σ 2 = β, by introducing the exponential transformation that is similar to Liao [14]: we can eliminate the convection term in Eq. (2.3) and transform it into For convenience, let σ 2 2 = s and 1 2σ 2 (r − σ 2 2 ) 2 + r = w, it is easy to see that s > 0, w > 0. Therefore the model (2.2) can be represented as the following time-fractional diffusion equation with the initial and boundary conditions:

Three new fourth-order compact finite difference schemes
The main purpose of this section is to construct new fourth-order compact difference schemes for problem (2.4)-(2.6). We start to introduce some definitions. Let t n = n · ∆t, n = 0, 1, 2, · · · , N ; x i = a + i · h, i = 0, 1, 2, · · · , M be the uniform time and space mesh, where ∆t = T N and h = b−a M are time step size and spatial step size respectively.
and define the grid function spaces as: x denotes the second-order central difference operator, that is For any ϕ ∈Φ h , define the operator A as

Algorithm 1
Eq. (2.4) at point (x i , t n ) can be formulated as: Applying L 1 formula (see [18]), the Caputo time- As for the second-order spacial derivatives, using Padé scheme, we have 3) into (3.1) and then multiplying the operator A on both sides of (3.1), we get There exists a positive constant C 1 such that In Eq.(3.4), denoting v n i as the approximate solution of v(x i , t n ) and omitting the higher term (r 1 ) n i , the implicit compact finite difference scheme (3.5) for problem (2.4)-(2.6) with initial and boundary discretizations is given as follows:

Algorithm 2
In Eq.(3.7), when n = 1,ĉ Employing linear interpolation between t n−1 and t n for terms and v(x i , t n−1+η ) in Eq.(3.6) respectively, we have Using Padé scheme for the second-order derivatives in space, we obtain and let A = I + h 2 12 δ 2 x . Eq.(3.6) is reformulated as the following form There exists a positive constant C 2 such that In Eq.(3.8), denoting v n i as the approximate solution of v(x i , t n ) and omitting the higher term (r 2 ) n i , we construct the following implicit discretization scheme (3.9) for problem (2.4)-(2.6) equipped with initial and boundary discretizations:

Algorithm 3
Eq.(2.4) at point (x i , t n ) can be represented as the following form (3.10) Using the L1 − 2 formula (see [10]), we discretize the Caputo time-fractional As for the spacial derivatives, using Padé scheme, we have There exists a positive constant C 3 such that In Eq.(3.11), denoting v n i as the approximate solution of v(x i , t n ) and omitting the higher term (r 3 ) n i , we obtain the following implicit discrete scheme (3.12) for problem (2.4)-(2.6) with initial and boundary discretizations as follows: (3.12)

Numerical experiments
In this section, an example with an exact solution is presented to demonstrate the high accuracy of the new schemes proposed in Section 3. Moreover the numerical results on the three schemes are compared. After that, we show the effectiveness of the method by applying Algorithm 3 to several different European option pricing problems.

Example 4.1. Consider the following time-fractional model with nonhomogeneous boundary conditions
is chosen such that the exact solution is V (x, t) = (t 3 + 1)(sin πx + 1).

M ax−error(∆t)
M ax−error(∆t/2) . Table 1 shows that the numerical errors by Algorithm 2 and Algorithm 3 are obviously much smaller than that by Algorithm 1, and the computational errors by Algorithm 3 are smaller than that by Algorithm 2. The corresponding temporal convergence orders of the three schemes are 2−α order in Algorithm 1, second-order in Algorithm 2 and 3 − α order in Algorithm 3 separately.

M ax−error(h)
M ax−error(h/2) . In order to eliminate the influence of the temporal approximation, the time step is used by 1/10000 in Algorithm 2 and Algorithm 3, whereas in Algorithm 1 the temporal stepsize is very very small to obtain its spatial convergence order. Table 2 shows that the spatial convergence orders of the three difference schemes are all indeed 4.
All of these are agree with theoretical analysis in references [23][24][25], which demonstrates the three compact difference schemes obtained in Section 3 are effective and accuracy.
Due to Staelen and Hendy have studied the numerical solution of the time fractional Black-Scholes model of order 0 < α < 1 in [17], where they constructed a numerical scheme of fourth order in space and 2 − α order in time. So it is necessary to compare that scheme with Algorithm 1 in this paper.
with the parameters : r = 0.5 and σ = √ is chosen such that the exact solution is V (x, t) = (t + 1) 2 (x 4 + x 2 + 1). The results are shown in Table 3 and Table 4.
In order to compare the two schemes, we use the same machine equipped with AMD 3.6-GHZ 4 Core processor and the codes are written in Matlab software. From Table 3 and Table 4, we can see that the two methods are both 2 − α order convergence in time and 4 order convergence in space. In addition, the computational errors of the two difference schemes are similar.  Let's take European put option as an example to illustrate the effect of different parameters on option price in time-fractional B-S model.

Example 4.3. Consider the time-fractional B-S model
with α = 0.5, and S a , S b are same as Example 4.3.   Figure  3(a) illustrates the effect of volatility of the stock price movement on option price. It can be seen that when the stock price is near the exercise price, the higher the volatility is, the higher the option price, which confirms a well-known statement in the real financial world: high risk, high return. Figure 3(b) shows the effect of risk-free interest rate on option price, here the parameters are T = 1, K = 20, σ = 0.3 and r = 0.005, 0.05, 0.1, 0.2, separately. From Figure 3(b) we can see that the higher the interest rate is, the lower the option will be.
Take r = 0.05, T = 1, σ = 0.3, K = 15, 20, 25 and 30, Figure 3(c) examines the effect of exercise price on the option price. As we can see, when the exercise price increases, the option price goes up too.
Finally, we analyse the effect of expiration date on the option price. The results are plotted in Figure 3(d), and the parameters are r = 0.05, K = 20, σ = 0.3, and T = 0.5, 1, 2, 3 years respectively. Figure 3(d) shows that when the stock price is low enough, an option with shorter expiration date is more profitable than an option with longer expiration date. While the stock price is high, the curve with longer expiration date is above the curve with shorter expiration date.
The above results match what happens in the real market very well.

Conclusion
In this work, the numerical approximation of the time-fractional Black-Scholes model has been studied. We firstly transformed the time-fractional B-S equation to a time-fractional diffusion equation by an exponential transformation. By using fourth-order Padé approximation to the second-order spatial derivatives, the spatial accuracy has been improved to 4. Then the time-fractional derivative was approximated by the L1 formula, L2 − 1 σ formula and L1 − 2 formula respectively, and we constructed three compact difference schemes with convergence rates O(∆t 2−α +h 4 ), O(∆t 2 + h 4 ) and O(∆t 3−α + h 4 ). Finally, numerical examples showed the accuracy and effectiveness of the proposed methods. The extension of the method to the Black-Scholes model with multiple degrees of freedom will be the future work for us. In addition, it will be also interesting to price other fractional models.