# Finite Difference/Fourier Spectral for a Time Fractional Black–Scholes Model with Option Pricing

Authors: Juan He; Aiqing Zhang
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

We study the fractional Black–Scholes model (FBSM) of option pricing in the fractal transmission system. In this work, we develop a full-discrete numerical scheme to investigate the dynamic behavior of FBSM. The proposed scheme implements a knownL1 formula for the α-order fractional derivative and Fourier-spectral method for the discretization of spatial direction. Energy analysis indicates that the constructed discrete method is unconditionally stable. Error estimate indicates that the 2−α-order formula in time and the spectral approximation in space is convergent with order OΔt2−α+N1−m, where m is the regularity of u and Δt and N are step size of time and degree, respectively. Several numerical results are proposed to confirm the accuracy and stability of the numerical scheme. At last, the present method is used to investigate the dynamic behavior of FBSM as well as the impact of different parameters.


---

## Body

## 1. Introduction

The classical option pricing model is proposed by Black and Scholes [1], which is based on the assumption that stocks and options are in an “ideal state” in the market and Samuelson’s model [2]:(1)dS=μSdt+σSdBt,where S be the the stock value, Bt is the Brownian motion with the unit variance, and μ and σ are two constants. But in many cases, fractional Brownian motion is more accurate than integer order [3]. On the other hand, more and more diffusion processes were found to be non-Fickian [4, 5], and the fractional order stochastic differential equation is considered as an extension of the stochastic differential equation. One view is that fractional order option trading equation is regarded as nonrandom growth process caused by Brownian motion. Therefore, Jumarie [6] and Liang et al. [7] considered fractional Brownian motion in Samuelson’s model equation:(2)dS=μSdt+σSdBt,α,dS=μSdt+σSwtdtα.Combining Itô lemma and fractional Taylor expansion of the option priceV, Jumarie [6] obtained the following FBSM:(3)∂αVS,t∂tα=rVS,t1−α!−rSα∂αVS,t∂Sαt1−α−α!31−α!22α!σ2S2∂2αVS,t∂S2α=0,∂αVS,t∂tα=rVS,t−rS∂V∂St1−α1−α!−α!2σ2S2∂2VS,t∂S2=0.Jumarie [8] derived new families of the exact solution of the above equations. Moreover, traditional pricing models for double barrier options are often biased when price changes are considered as fractal transmission systems. Chen et al. [9] revealed that it would be better to use fractional order Black–Scholes equation to explain the pricing in fractal transmission systems.There is a lot of work in the modeling and calculation of fractional equations. Yang et al. [10, 11] developed a new definition of fractional derivative. The advantage of this definition is that it does not contain singular kernel. Inc et al. [12, 13] studied the isolated solutions of a class of fractional equations with Kerr law nonlinearity by Riccati–Bernoulli method. The exact dark optical and periodic singular soliton solution is obtained. Singhet al. [14–16] solved a series of fractional equations by using homotopy analysis technique and Laplace transform algorithm.As we all know, as an effective method,L1 formula [17] has been widely used in the calculation of fractional differential equations. Langlands and Henry [18], Sun and Wu [19], and Lin and Xu [20] discussed the error estimate for L1 scheme. De Staelena and Hendybc [21] constructed a numerical method of fourth-order finite difference in space and 2−α in time. Stability, uniqueness, and error estimates are analyzed. Zhang et al. [22] presented a discrete implicit finite scheme to solve time fractional Black–Scholes model. They discussed the stability and error estimation of numerical schemes by Fourier analysis. Unfortunately, their analysis methods are local approach. Due to the importance of FBSM, it is necessary to reconstruct an efficient numerical method and analyze global stability and error estimates.In this work, we will develop an efficient full-discrete scheme to approximate the Black–Scholes model withα-order fractional derivative. We apply L1 method to discretize the direction of time and Fourier-spectral method to discretize the direction of space. Using the energy analysis method, we discuss the stability and error estimate of the fully discrete numerical method. The detailed analysis shows that the scheme is unconditionally energy stable, and the error estimates indicate that our full-discrete scheme can achieve 2−α-order accuracy in time and exponential accuracy in space direction. Finally, some numerical examples are conducted to support the theoretical claims. At the same time, the dynamic behavior of FBSM is studied by the proposed method.We organize the rest of the paper as follows. Section2 will briefly introduce the FBSM. In Section 3, we develop a time-discrete method for FBSM and then present its discrete energy law. In Section 4, we will study the error estimate of the full-discrete scheme. In Section 5, we present accuracy/stability tests and numerous numerical examples to demonstrate the validity of the full-discrete method. In addition, we will discuss the properties of the solution of the FBSM. Some concluding remarks are given in Section 6.


## 2. Black–Scholes Model

In this work, we will consider the following time fractional Black–Scholes model:(4)∂αVS,t∂tα+12σ2S2∂2VS,t∂S2+rS∂VS,t∂S−rVS,t=0,S,t∈0,+∞×0,T,(5)V0,t=pt,V+∞,t=qt,(6)VS,T=wS,where V is the price of the option, S the price of the underlying asset, r the interest rate, and σ the volatility of the stock price, 0<α≤1; the time fractional derivative ∂αVS,t/∂tα is defined by(7)∂αV⋅,t∂tα:=∂tαV⋅,t=1Γ1−α∫0tt−μ−α∂V⋅,μ∂μdμ.This is a linear parabolic partial differential equation which has been studied extensively.We transform the problem to an initial value problem by using the time to maturet:=T−τ, and we then set x=lnS,VS,τ=uex,T−τ; we can rewrite (4) as(8)∂tαux,τ−ξ∂x2ux,τ−ω∂xux,τ+rux,τ=0,x,τ∈0,+∞×0,T,where ξ=1/2σ2,ω=r−ξ, and with the following boundary (barrier) and initial conditions,(9)u−∞,τ=pτ,u−∞,τ=qτ,ux,0=u0x,a<x<b.In order to solve the above model by numerical method, it is necessary to truncate the original unbounded region into a finite interval. Therefore, we will consider problem (8) in bounded interval 0,2π. Then, we will study the following problem:(10)∂ταux,τ−ξ∂x2ux,τ−ω∂xux,τ+rux,τ=0,x,τ∈0,2π×0,T,(11)u0,τ=u2π,τ,(12)ux,0=u0x.Remark 1.
In fact, one can choose homogeneous or inhomogeneous boundary conditions. It all depends on the actual option price. We have tested it, and it does not make any difference in actual numerical examples.


## 3.2−αOrder Numerical Method

Here, we will develop the time-discrete method for equation (10). First, given a positive integer K, set Δt=T/K be the time step size, and denote τn=nΔt,0≤n≤K as the mesh point. Then, we introduce an L1 method to discrete the Caputo fractional derivative of order α:(13)∂tαux,τ=1Γ2−α∑j=0nbju⋅,τn+1−j−u⋅,τn−jΔtα+OΔt2−α,where bj=j+11−α−j1−α.Lemma 1 (see [20, 23, 24]).
The coefficients of formula (13) satisfy the following properties:(14)1=b0>b1>b2>⋯bj⟶0,asj⟶+∞.Then, we can obtain the following time-discrete scheme:(15)1a0un+1−∑j=0n−1bj−bj+1un−j−bnu0−ξ∂x2un+1−ω∂xun+1+run+1=0,n≥0,where a0=ΔtαΓ2−α. It should be noted that if n=0, we can rewrite the above equation as(16)1a0u1−u0−ξ∂x2u1−ω∂xu1+ru1=0.First of all, we have the following energy stability results for time-discrete (15).Theorem 1.
The time-discrete scheme (15) is unconditionally stable. It satisfies the following energy dissipation law:(17)uk+1≤u0,k=0,1,…,K.Proof.
Whenn=0, computing the L2 inner product of (16) with 2u1, we obtain(18)2u1−u0,u1−2a0ξ∂x2u1,u1−a0ω∂xu1,u1+a0ru1,u1=0.
It is easy to verify that the following formula is correct:(19)2A−B,A=A2−B2+A−B2.
Thus,(20)u12−u02+u1−u02+2a0a∂xu12+2a0ru12=0.
Giving up some positive terms, we have(21)u12≤u02.
Assume the following inequality holds:(22)uj2≤u02,j=2,3,…,n.
Next, we will showun+12≤u02 is still valid. If j=n+1, taking the L2 inner product of (15) with 2un+1, we derive(23)2un+12+2a0ξ∂xun+12+2a0run+12=2∑j=0n−1bj−bj+1un−j+bnu0,un+1≤∑j=0n−1bj−bj+1un−j2+un+12+bnu02+un+12.
Note the fact that(24)∑j=0n−1bj−bj+1+bn=1.
Thus, we get(25)un+12≤∑j=0n−1bj−bj+1un−j2+bnu02≤∑j=0n−1bj−bj+1+bnu02=u02.
This yields (17).


## 4. Error Estimate for Full Discretization

In this part, we will study the Fourier-spectral method for the time-discrete method (15). First, we define SN as the polynomial space. Define πN:L2Ω⟶SN be the L2-projection operator which satisfies(26)πNϕ−ϕ,ψ=0,∀ψ∈SN.We have the following estimate [25]:(27)ϕ−πNϕl≤CNl−mϕm,∀ϕ∈HmΩ,m>l≥0.Then, we can develop the following full-discrete scheme:(28)1a0uNn+1−∑j=0n−1bj−bj+1uNn−j−bnuN0,ϕN−ξ∂x2uNn+1,ϕN−ω∂xun+1,ϕN+run+1,ϕN=0,ϕN∈SN.We now present the stability results of the fully discrete scheme (28).Theorem 2.
LetuNnn=1M−1 be the solution of (28), then we derive(29)uNn+12≤uN02.Next, we begin to analyze the error estimates of the full-discrete scheme (28). Define the following error function:(30)Rn+1=1Γ2−α∑j=0nbju⋅,tn+1−j−u⋅,tn−jΔtα−1Γ1−α∫0tn+1∂u⋅,μ∂μdμtn+1−μα.From [20, 19, 24, 23], we know that Rn+1 satisfies(31)Rn+1≤CΔt2−α.We also define the following error functions:(32)e˜un=πNutn−uNn,e^un=utn−πNutn,eun=e^un+e˜un=utn−uNn.Lemma 2.
Fora0 and bn, we have the following results:(33)a0≤2bnΓ1−αTα.Proof.
(34)bn=n+11−α−n1−α=n1−α1+1n1−α−1=n1−α1−αn+1−α−α2!1n2+⋯≥n1−α1−αn+1−α−α2!1n2=n1−α1−α2n+1−α2n+1−α−α2!1n2.
Note that(35)1−α2n+1−α−α2!1n2≥0.
Therefore, we obtain(36)a0bn≤2nαΔtαΓ2−α1−α≤2Γ1−αTα.Remark 2.
In [23, 24], readers will also find similar parameter estimation. This result is very useful for us to analyze error estimate.Theorem 3.
For the constructed numerical scheme (28), we have the following error estimate:(37)uτk−uk≤CΔt2−α+N1−m,k=0,1,…,K=TΔt.Proof.
Forn=0, we can write equation (28) as(38)1a0uN1−uN0,ϕN−ξ∂x2uN1,ϕN−ω∂xuN1,ϕN+ruN1,ϕN=0.
Subtracting (38) from (10) at τ1, we note that(39)∂ταu⋅,τ1−uN1−uN0a0=∂ταu⋅,τ1−u⋅,τ1−u⋅,τ0a0+u⋅,τ1−u⋅,τ0a0−πNu⋅,τ1−πNu⋅,τ0a0+πNu⋅,τ1−πNu⋅,τ0a0−uN1−uN0a0=R1+1a0I−πNu⋅,τ1−u⋅,τ0+1a0e˜u1−e˜u0,∂x2u⋅,τ1−∂x2uN1=∂x2u⋅,τ1−πN∂x2u⋅,τ1+πN∂x2u⋅,τ1−∂x2uN1=∂x2I−πNu⋅,τ1+∂x2e˜u1,∂xu⋅,τ1−∂xuN1=∂xu⋅,τ1−πN∂xu⋅,τ1+πN∂xu⋅,τ1−∂xuN1=∂xI−πNu⋅,τ1+∂xe˜u1,u⋅,τ1−uN1=u⋅,τ1−πNu⋅,τ1+πNu⋅,τ1−uN1=I−πNu⋅,τ1+e˜u1.
Then, we have(40)e˜u1−e˜u0,ϕN+a0ξ∂xe˜u1,∂xϕN−a0ω∂xe˜u1,ϕN+a0re˜u1,ϕN=−a0R1,ϕN+πN−Iu⋅,τ1−u⋅,τ0,ϕN+a0ξπN−I∂xu⋅,τ1,∂xϕN−a0ωI−πNu⋅,τ1,∂xϕN+a0rπN−Iu⋅,τ1,ϕN,ϕN∈SN.
SetϕN=2e˜u1, we have(41)2e˜u12+2a0ξ∂xe˜u12+2a0re˜u12≤a01rb0R12+rb0e˜u12+a0ξ1b0πN−I∂xu⋅,τ12+b0∂xe˜u12+a0ωωξb0I−πNu⋅,τ12+ξb0ω∂xe˜u12.
Dropping some positive terms, we find(42)e˜u12≤C1a0b0Δt4−2α+C2a0b0N1−m+C3a0b0N−m.
Assume(43)e˜uk2≤C1a0bk−1Δt4+C2a0bk−1N1−m+C3a0bj−1N−m,k=2,3,…,n.
Next, we will prove that it holds also fork=n+1. Subtracting (28) from a reformulation of (10) at tn+1, we find(44)e˜un+1−∑j=0n−1bj−bj+1e˜un−j−bne˜u0,ϕN+a0ξ∂xe˜un+1,∂xϕN−a0ω∂xe˜un+1,ϕN+a0re˜un+1,ϕN=−a0Rn+1,ϕN+πN−Iu⋅,τn+1−∑j=0n−1bj−bj+1u⋅,τn−j−u⋅,τ0,ϕN+a0ξπN−I∂xu⋅,τn+1,∂xϕN−a0ωI−πNu⋅,τn+1,∂xϕN+a0rπN−Iu⋅,τn+1,ϕN.
LetϕN=2e˜un+1, we have(45)2e˜un+12+2a0ξ∂xe˜un+12+2a0re˜un+12≤2∑j=0n−1bj−bj+1e˜un−j−bne˜u0,e˜un+1−2a0Rn+1,e˜un+1+2a0ξπN−I∂xu⋅,τn+1,∂xe˜un+1−2a0ωI−πNu⋅,τn+1,∂xe˜un+1≤∑j=0n−1bj−bj+1e˜un−j2+e˜un+12+a01rRn+12+re˜un+12+a0ξπN−I∂xu⋅,τn+12+∂xe˜un+12+a0ωωξI−πNu⋅,τn+12+ξω∂xe˜un+12.
Thus, we have(46)e˜un+12≤∑j=0n−1bj−bj+1e˜un−j2+a0rRn+12+a0aπN−I∂xu⋅,τn+12+a0b2aI−πNu⋅,τn+12≤∑j=0n−1bj−bj+1C1a0bj−1Δt4+C2a0bj−1N1−m+C3a0bj−1N−m+bnC1a0bnΔt4+C2a0bnN1−m+C3a0bnN−m.
Note thatbn−j≥bn, thus(47)e˜un+12≤C1a0bnΔt4+C2a0bnN1−m+C3a0bnN−m∑j=0n−1bj−bj+1+bn=C1a0bnΔt4+C2a0bnN1−m+C3a0bnN−m.
Note that(48)uτk−uk≤e^uk+e˜uk.
This ends the proof.


## 5. Numerical Examples

In this section, several numerical examples will be present to confirm the accuracy and applicability of the full-discrete scheme (28). We consider a rectangular computed domain of 0,2π. In order to better simulate the periodic boundary conditions, the Fourier-spectral method will be used to discretize space direction.


### 5.1. Verification of Convergence of Numerical Method

First, in order to conduct a time accuracy test, an exact solution will be constructed to evaluate the convergence of the full-discrete scheme (28).Example 1.
We consider the following FBSM withα-order Caputo derivative:(49)∂ταux,τ−ξ∂x2ux,τ−ω∂xux,τ+rux,τ=fx,τ,where(50)fx,τ=2T2−αΓ3−αsinx+ξτ2sinx−ωτ2cosx+rτ2sinx.
It is easy to verify that the exact solution will beu=τ2sinx.
We setN=128. The default values for the parameters are set as T=1,ξ=1,ω=0,andr=1. In Table 1, we show the temporal convergence orders of various time steps. As can be seen from Table 1, our full-discrete scheme is close to 2−α-order accuracy in time, which is confirm with the result in Theorem 3.Table 1
Temporal convergence orders of various time steps for Example1.
α/ΔtΔt = 1/200Δt = 1/400Δt = 1/800Δt = 1/1600Δt = 1/3200α=0.11.81661.82631.83451.84151.8475α=0.31.66801.67461.67971.68371.6869α=0.51.48921.49251.49471.49631.4974α=0.61.39391.39601.39741.39831.3989α=0.71.29641.29791.29871.29921.2995α=0.91.09801.09901.09951.09971.0999FixN=128, in Figure 1, we give L2 error for different α. It is obvious that our numerical scheme has good convergence in time direction. Let T=0.1,Δt=10−6,ξ=r=0.5,ω=0, and u0=cos8x. Figure 2 shows that the full-discrete scheme (28) has excellent convergence behavior in space direction.Figure 1
Error in time direction for differentα.Figure 2
TheL2 and L∞ errors in space direction with α=0.1.


### 5.2. Effect of Various Parameters

This section is devoted to investigate the dynamic behavior of FBSM equation with differentα. In the following numerical experiments, we fix ξ=0.5,w=0.1,r=0.6,N=64,Δt=0.01,andu0=sinx. From Figures 3–5, we know that α has certain influence on the solution, as α increases, and the solution becomes smoother. In order to test the influence of ξ on the option price, we set r=0.6,u0=sinx/2 and let ξ change at the same time. Figure 6 shows that the parameter ξ has a significant effect on the price of options, and there will be an inflection point around x=2.8. Finally, we investigated the influence of r on the option price, and the results are shown in Figure 7, it can be seen that when r increases, the option price also increases.Figure 3
The dynamic behavior of solution to FBSM equation atα=0.2.Figure 4
The dynamic behavior of solution to FBSM equation atα=0.5.Figure 5
The dynamic behavior of solution to FBSM equation atα=0.9.Figure 6
The influence of variousξξ=0.1,0.2,0.3 on the option price.Figure 7
The influence of variousrξ=0.3,0.4,0.5 on the option price.


## 5.1. Verification of Convergence of Numerical Method

First, in order to conduct a time accuracy test, an exact solution will be constructed to evaluate the convergence of the full-discrete scheme (28).Example 1.
We consider the following FBSM withα-order Caputo derivative:(49)∂ταux,τ−ξ∂x2ux,τ−ω∂xux,τ+rux,τ=fx,τ,where(50)fx,τ=2T2−αΓ3−αsinx+ξτ2sinx−ωτ2cosx+rτ2sinx.
It is easy to verify that the exact solution will beu=τ2sinx.
We setN=128. The default values for the parameters are set as T=1,ξ=1,ω=0,andr=1. In Table 1, we show the temporal convergence orders of various time steps. As can be seen from Table 1, our full-discrete scheme is close to 2−α-order accuracy in time, which is confirm with the result in Theorem 3.Table 1
Temporal convergence orders of various time steps for Example1.
α/ΔtΔt = 1/200Δt = 1/400Δt = 1/800Δt = 1/1600Δt = 1/3200α=0.11.81661.82631.83451.84151.8475α=0.31.66801.67461.67971.68371.6869α=0.51.48921.49251.49471.49631.4974α=0.61.39391.39601.39741.39831.3989α=0.71.29641.29791.29871.29921.2995α=0.91.09801.09901.09951.09971.0999FixN=128, in Figure 1, we give L2 error for different α. It is obvious that our numerical scheme has good convergence in time direction. Let T=0.1,Δt=10−6,ξ=r=0.5,ω=0, and u0=cos8x. Figure 2 shows that the full-discrete scheme (28) has excellent convergence behavior in space direction.Figure 1
Error in time direction for differentα.Figure 2
TheL2 and L∞ errors in space direction with α=0.1.


## 5.2. Effect of Various Parameters

This section is devoted to investigate the dynamic behavior of FBSM equation with differentα. In the following numerical experiments, we fix ξ=0.5,w=0.1,r=0.6,N=64,Δt=0.01,andu0=sinx. From Figures 3–5, we know that α has certain influence on the solution, as α increases, and the solution becomes smoother. In order to test the influence of ξ on the option price, we set r=0.6,u0=sinx/2 and let ξ change at the same time. Figure 6 shows that the parameter ξ has a significant effect on the price of options, and there will be an inflection point around x=2.8. Finally, we investigated the influence of r on the option price, and the results are shown in Figure 7, it can be seen that when r increases, the option price also increases.Figure 3
The dynamic behavior of solution to FBSM equation atα=0.2.Figure 4
The dynamic behavior of solution to FBSM equation atα=0.5.Figure 5
The dynamic behavior of solution to FBSM equation atα=0.9.Figure 6
The influence of variousξξ=0.1,0.2,0.3 on the option price.Figure 7
The influence of variousrξ=0.3,0.4,0.5 on the option price.


## 6. Conclusion

In this paper, a new full-discrete numerical method is developed to solve the FBSM. An efficient2−α-order and unconditionally energy stable method is constructed by combining the L1 approach in time and Fourier method in space direction. It is proved that the full-discrete converges to the order OΔt2−α+N−s+N−m globally. Numerical examples demonstrate the robustness and accuracy of the developed full-discrete method, numerically. Finally, we also study the properties of the solution of the FBSM.



---

*Source: 1393456-2020-09-04.xml*