# Spectral Solutions for Fractional Black–Scholes Equations

Authors: M. A. Abdelkawy; António M. Lopes
Journal: Mathematical Problems in Engineering
 (2022)

---

## Abstract

This paper presents a numerical method to solve accurately the fractional Black–Scholes model of pricing evolution. A fully spectral collocation technique for the two independent variables is derived. The shifted fractional Jacobi–Gauss–Radau and shifted fractional Jacobi–Gauss–Lobatto collocation techniques are utilized. Firstly, the independent variables are interpolated at the shifted fractional Jacobi nodes, and the solution of the model is approximated by means of a sequence of shifted fractional Jacobi orthogonal functions. Then, the residuals at the shifted fractional Jacobi quadrature locations are estimated. As a result, an algebraic system of equations is obtained that can be solved using any appropriate approach. The accuracy of the proposed method is demonstrated using two numerical examples. It is observed that the new technique is more accurate, efficient, and feasible than other approaches reported in the literature. Indeed, the results show the exponential convergence of the method, both for smooth and nonsmooth solutions.


---

## Body

## 1. Introduction

Fractional calculus [1–3] is the branch of calculus that generalizes the derivative and integral operators to noninteger orders. During the last few decades, the fractional calculus gained increasing importance in several fields, such as mechanics [4, 5], biology [6, 7], physics [8], chemistry [9, 10], viscoelasticity [11], engineering [12], and finance/economics [13, 14]. Fractional differential equations can be used to represent accurately a variety of phenomena. However, many equations cannot be solved analytically, and, thus, numerical approaches are required to find approximate solutions. The finite difference [15, 16] and finite element [17–20] methods have been widely used, as well as other numerical schemes, such as the Galerkin [21], Petrov-Galerkin [22, 23], pseudospectral [24–26], and fast predictor-corrector [27] methods.The Black–Scholes model (BSM) is a mathematical description of pricing evolution [28, 29]. Specifically, the model estimates financial instruments via time variation. These instruments (such as stocks or futures) track a lognormal distribution of prices. Based on this hypothesis and taking into account other variables, the BSM derives the price of call options. Despite there being several models that describe pricing evolution, the BSM has been one of the most significant and prevalent in the last decades, and its generalization to fractional order was proposed [30, 31].In [32], the Laplace homotopy analysis approach was employed to solve the fractional BSM (FBSM) in the sense of the Caputo–Fabrizio derivative. In [33], the authors employed the Sumudu and Laplace transforms to address various economic models with different fractional operators. In [34, 35], the Adomian and fractional Adomian decomposition were adopted to solve the FBSM. In [36], the time-fractional BSM was solved using the generalized differential transform technique. Moreover, in [37], the existence and uniqueness solution of the European-type option pricing model was discussed, while in [38], the multivariate Padé approximation was utilized to solve the Caputo fractional European vanilla call option pricing problem.In most circumstances, it is impossible to provide explicit analytical solutions to space and/or time-fractional differential equations. Hence, developing effective numerical techniques is a critical necessity. Many high-accuracy numerical approaches have been devised to tackle diverse issues in various applications. Because of their high-order precision, fractional differential and integral equations [39] have seen tremendous development in recent years. But, contrasting with the effort put into evaluating finite difference and element schemes, little research has been dedicated to designing and assessing global spectral schemes [40]. In cases with periodic boundary conditions, spectral techniques based on the Fourier expansion have been used. When compared to other schemes, spectral techniques take the main place due to their robustness and exponential rates of convergence [41–44]. Despite some drawbacks, such as the inability to represent physical processes in spectral space and the difficulty in parallelizing on distributed memory computers, they reveal excellent accuracy, high speed of convergence, and simplicity in solving many types of differential equations [45]. The spectral collocation method is particularly valuable because it can estimate the solution of a wide range of equations. Moreover, its exponential rate of convergence is extremely useful in delivering very exact solutions. The collocation approach has grown in prominence in recent decades for dealing with specific difficulties posed by fractional derivatives.In this paper, a new spectral collocation approach is proposed to solve the FBSM [31, 46, 47], given by(1)DμZx,τ+12α2x2∂2∂x2Zx,τ+rx∂∂xZx,τ−rYx,τ+Gx,τ=0,x,τ∈℧×℧,with DμZx,τ denoting the right Riemann–Liouville fractional derivative:(2)DμZx,τ=1Γ1−μ∂∂τ∫τTZx,s−Zx,Ts−τμds,where 0<μ<1.Herein, the shifted fractional Jacobi–Gauss–Radau collocation (SFJ-GR-C) and shifted fractional Jacobi–Gauss–Lobatto collocation (SFJ-GL-C) techniques are utilized to solve the FBSM. We interpolate the independent variables using the shifted fractional-order Jacobi nodes, and we approximate the solution of the model by a sequence of shifted fractional-order Jacobi orthogonal functions. After that, the residuals at the shifted fractional-order Jacobi quadrature locations are estimated. This yields an algebraic system of equations that can be solved using a suitable approach. The accuracy of the new method is demonstrated using numerical examples.The paper is organized into four sections: Section2 introduces the new spectral collocation technique. Section 3 applies the method to solve numerical examples. Section 4 summarizes the main conclusions.


## 2. Fully Spectral Collocation Technique

We start by mapping the variableτ as τ=t−T. Additionally, the Riemann–Liouville is transformed to the Caputo fractional derivative DμYx,τ (see [31, 46]):(3)∂μ∂tμYx,t−12α2x2∂2∂x2Yx,t−rx∂∂xYx,t+rYx,t=ℍx,t,x,t∈℧×℧,where ℧≡0,ℒ and ℧≡0,T. Conditions are(4)Yx,0=Θ1x,x∈℧,Y0,t=Θ2t,Yℒ,t=Θ3t,t∈℧,and DμZx,τ is a Caputo fractional derivative.We express the truncated solution as(5)YN,ℳx,t=∑r1=0,…,Nr2=0,…,ℳςr1,r2Gℒ,r1,λ1α1,β1xGT,r2,λ2α2,β2t,where Gℒ,r,λα,βx denotes the shifted fractional Jacobi functions on 0,ℒ (see [48, 49] for extra information). ∂/∂xYN,ℳx,t is determined as [48, 49](6)∂∂xYN,ℳx,t=∑r1=0,…,Nr2=0,…,ℳςr1,r2Gℒ,r1,λ1α1,β1xGT,r2,λ2α2,β2t.Similarly, we have [48, 49](7)∂2∂x2YN,ℳx,t=∑r1=0,…,Nr2=0,…,ℳςr1,r2Gℒ,r1,λ1′′α1,β1xGT,r2,λ2α2,β2t.The time Caputo fractional derivative, on the contrary, is obtained as [48, 49](8)∂μ∂tμYN,ℳx,t=∑r1=0,…,Nr2=0,…,ℳςr1,r2Gℒ,r1,λ1α1,β1xGT,r2,λ2μα2,β2t.Previous computations are indicated at certain nodes as(9)∂μ∂tμYx,tt=tT,ℳ,mα2,β2x=xℒ,N,nα1,β1=∑r1=0,…,Nr2=0,…,ℳςr1,r2Gℒ,r1,λ1α1,β1xℒ,N,nα1,β1GT,r2,λ2μα2,β2tT,ℳ,mα2,β2,∂∂xYN,ℳx,tt=tT,ℳ,mα2,β2x=xℒ,N,nα1,β1=∑r1=0,…,Nr2=0,…,ℳςr1,r2Gℒ,r1,λ1′α1,β1xℒ,N,nα1,β1GT,r2,λ2α2,β2tT,ℳ,mα2,β2,∂2∂x2YN,ℳx,tt=tT,ℳ,mα2,β2x=xℒ,N,nα1,β1=∑r1=0,…,Nr2=0,…,ℳςr1,r2Gℒ,r1,λ1′′α1,β1xℒ,N,nα1,β1)GT,r2,λ2α2,β2tT,ℳ,mα2,β2.Otherwise, the initial boundary can be obtained by(10)∑r1=0,…,Nr2=0,…,ℳςr1,r2Gℒ,r1,λ1α1,β1xGT,r2,λ2α2,β20=Θ1x,∑r1=0,…,Nr2=0,…,ℳςr1,r2Gℒ,r1,λ1α1,β10GT,r2,λ2α2,β2t=Θ2t,∑r1=0,…,Nr2=0,…,ℳςr1,r2Gℒ,r1,λ1α1,β1ℒGT,r2,λ2α2,β2t=Θ3t.Eq. (3) is constrained to be zero at N−1×ℳ points:(11)Δxℒ,N,nα1,β1,tT,ℳ,mα2,β2=ℍxℒ,N,nα1,β1,tT,ℳ,mα2,β2,where(13)Δxℒ,N,nα1,β1,tT,ℳ,mα2,β2=∑r1=0,…,Nr2=0,…,ℳςr1,r2Gℒ,r1,λ1α1,β1xℒ,N,nα1,β1GT,r2,λ2μα2,β2tT,ℳ,mα2,β2−12α2xℒ,N,nα1,β12∑r1=0,…,Nr2=0,…,ℳςr1,r2Gℒ,r1,λ1′′α1,β1xℒ,N,nα1,β1GT,r2,λ2α2,β2tT,ℳ,mα2,β2−rxℒ,N,nα1,β1∑r1=0,…,Nr2=0,…,ℳςr1,r2G′ℒ,r1,λ1α1,β1xℒ,N,nα1,β1GT,r2,λ2α2,β2tT,ℳ,mα2,β2,where n=1,…,N−1,andm=1,…,ℳ, and(14)∑r1=0,…,Nr2=0,…,ℳςr1,r2Gℒ,r1,λ1α1,β1xℒ,N,nα1,β1GT,r2,λ2α2,β20=Θ1xℒ,N,nα1,β1,∑r1=0,…,Nr2=0,…,ℳςr1,r2Gℒ,r1,λ1α1,β10GT,r2,λ2α2,β2tT,ℳ,mα2,β2=Θ2tT,ℳ,mα2,β2,∑r1=0,…,Nr2=0,…,ℳςr1,r2Gℒ,r1,λ1α1,β1ℒGT,r2,λ2α2,β2tT,ℳ,mα2,β2=Θ3tT,ℳ,mα2,β2.Finally, a solvable linear algebraic system is provided.It should be mentioned that the convergence spectral rate of Jacobi polynomials has been thoroughly studied in the literature [50, 51], while the convergence spectral rate of fractional Jacobi functions has been addressed in [48, 49, 52], to cite a few.


## 3. Numerical Results

In this section, we demonstrate the resilience and accuracy of our strategy to solve two problems. The algorithm was implemented in MATHEMATICA version 12.2, running in an Intel(R) Core(TM) i7-10510U CPU @ 2.30 GHz, on a PC with 8.00 GB of RAM. In what follows, we take℧≡0,1 and ℧≡0,1.Example 1.
We solve the FBSM [47].(15)∂μ∂tμYx,t−12α2∂2∂x2Yx,t−r1∂∂xYx,t+r2Yx,t=ℍx,t,with(16)Yx,0=cosπx,x∈0,1,Y0,t=1+t6,Y1,t=−t6−1,t∈0,1,choosing ℍx,t such that(17)Yx,t=t6+1sinπx+π2.
Table1 summarizes the L∞ errors between the exact and approximate solutions achieved by our method and those reported in [47], taking α=0.25,r1=0.01875,r2=0.05, and μ=0.6. We verify that the new strategy is superior and that good estimates are obtained with a small number of collocation points. Figures 1 and 2 present the 3D charts of the numerical solution and the absolute error, respectively, when α1=β1=α2=β2=0,λ1=1,λ2=1, and N=ℳ=20. Figure 3 depicts the corresponding exact and numerical solutions along x. Figures 4 and 5 represent the absolute errors along x and t, respectively, while Figure 6 shows the maximum absolute error ME convergence for the following cases:Case I:α1=1/2,α1=1/3,α2=2/3,β2=1/2,λ1=0.1,λ2=1.Case II:α1=α1=α2=β2=1/2,λ1=0.5,λ2=1.Case III:α1=α1=−α2=−β2=1/2,λ1=0.9,λ2=1, λ2=1.Table 1
TheL∞ errors between the exact and approximate solutions achieved by our method and those reported in [47] for Example 1, taking α=0.25,r1=0.01875,r2=0.05, and μ=0.6.
N,ℳα1=β1=0α1=1/2,β1=0α1=β1=1/2Method in [47] for μ=0.6α2=β2=0α2=−1/2,β2=0α2=β2=−1/2N,ℳL∞(4, 4)2.55732×10−21.09933×10−21.26953×10−2(7,20)8.6814e−07(8, 8)1.60814×10−61.8991×10−61.446×10−6(9,20)4.3312e−09(12, 12)3.43343×10−113.95182×10−112.89728×10−11(11,20)1.8476e−11(16, 16)4.21885×10−155.32907×10−153.58257×10−15(13,20)5.7662e−14(20, 20)4.08849×10−154.66294×10−156.66134×10−16Figure 1
Numerical solution of Example1, when α1=β1=α2=β2=0,λ1=1,λ2=1, and N=ℳ=20.Figure 2
Absolute error of Example1, when α1=β1=α2=β2=0,λ1=1,λ2=1, and N=ℳ=20.Figure 3
Exact and numerical solutions along with Example1, where α1=β1=α2=β2=0,λ1=1,λ2=1, and N=ℳ=20.Figure 4
Absolute error alongt with Example 1, when α1=1/2β1=0=α2=−1/20,β2=0,λ1=1,λ2=1, and N=ℳ=20.Figure 5
Absolute error alongx with Example 1, when α1=1/2,β1=0,α2=−1/2,β2=0,λ1=1,λ2=1, and N=ℳ=20.Figure 6
TheME convergence of Example 1 for three cases.Example 2.
We solve the FBSM [46].(18)∂μ∂tμYx,t−12α2x2∂2∂x2Yx,t−rx∂∂xYx,t+rYx,t=ℍx,t,with(19)Yx,0=ex+x+1,x∈0,1,Y0,t=tμ+1,Y1,t=tμ+e+2,t∈0,1,choosing ℍx,t such that Yx,t=tμ+ex+x+1.
Table2 summarizes the L∞ errors between the exact and approximate solutions achieved by our method and those reported in [46], taking α=0.1,r=0.06, and 0<μ<1. We verify that the new strategy is superior and that good estimates are obtained with a small number of collocation points. Figures 7 and 8 present the 3D charts of the numerical solution and the absolute error, respectively, when α1=−α2=−0.5,β1=β2=0,λ1=0.9,λ2=1, and N=ℳ=14. Figure 9 depicts the corresponding exact and numerical solutions along x. Figures 10 and 11 represent the absolute errors along x and t, respectively, while Figure 12 shows the ME convergence for the following cases:Case I:α1=β1=α2=β2=0,μ=0.6,λ1=1,λ2=1.Case II:α1=0.5,β1=0,α2=−1/2,β2=0,μ=0.6,λ1=1,λ2=1.Case III:α1=β1=−α2=−β2=1/2,μ=0.6,λ1=1,λ2=1.
Figures8, 10, and 11 show that the results are fairly precise, as the absolute error approaches zero. We can also see in Figures 7 and 9 that the approximate solution matches precisely the exact solution. In addition, Figure 12 clearly reveals that the method yields exponential convergence of the error.
To sum up, we emphasize that the majority of numerical techniques for the problem at hand are based on orthogonal polynomials. The use of classical polynomials in problems with nonsmooth solutions leads to low accuracy or even failure to converge. Using the fractional, rather than the classical, Jacobi functions mitigates the problem. With the proposed method, all numerical computations could be completed with good precision and a low number of degrees of freedom. Moreover, our technique outperformed other existing approaches. Finally, we can conclude that the fully spectral collocation approach is a useful, efficient, and acceptable strategy for dealing with problems that have singular solutions.Table 2
TheL∞ errors between the exact and approximate solutions achieved by our method and those reported in [46] for Example 2, taking α=0.1,r=0.06, and 0<μ<1.
N,ℳμ=0.1μ=0.5μ=0.9α1=1/2,β1=0α1=1/3,β1α1=β1=0α1=β1=1/2α1=1/2,β1=0α1=β1=1/2α2=β2=1/2α2=1/2,β2=1/2α2=β2=0α2=β2=1/2α2=−1/2,β2=0α1=β2=−1/2(2,2)7.02×10−31.20449×10−32.49091×10−32.48772×10−37.29645×10−32.48377×10−3(6,6)5.66938×10−84.69012×10−85.81708×10−84.75405×10−85.66938×10−84.75405×10−8(10,10)4.04121×10−142.30926×10−143.06422×10−142.57572×10−144.13003×10−142.75335×10−14(14,14)5.77316×10−151.33227×10−151.33227×10−152.66454×10−153.10862×10−153.00862×10−15Standard central difference method [46]N,ℳ641282565121024μ0.13.3676e−41.3450e−46.0422e−52.9532e−51.5110e−50.56.8645e−42.9990e−41.3884e−46.6733e−53.2753e−50.96.7548e−42.9157e−41.3344e−46.3528e−53.0953e−5Figure 7
Numerical solution of Example2, when α1=−α2=−0.5,β1=β2=0,λ1=0.9,λ2=1, and N=ℳ=14.Figure 8
Absolute error of Example2, when α1=−α2=−0.5,β1=β2=0,λ1=0.9,λ2=1, and N=ℳ=14.Figure 9
Exact and numerical solutions alongx of Example 2, when α1=−α2=−0.5,β1=β2=0,λ1=0.9,λ2=1, and N=ℳ=14.Figure 10
Absolute error alongt of Example 2, when α1=−α2=−0.5,β1=β2=0,λ1=0.9,λ2=1, and N=ℳ=14.Figure 11
Absolute error alongx of Example 2, when α1=−α2=−0.5,β1=β2=0,λ1=0.9,λ2=1, and N=ℳ=14.Figure 12
TheME convergence of Example 2 for the three cases.


## 4. Conclusion

The FBSM was treated using a fully spectral collocation technique for the two independent variablesx and t. It was verified that the new technique is superior in terms of accuracy and efficiency to other methods, for both smooth and nonsmooth solutions. To deal with the FBSM, we devised an approach that yields an algebraic system from which an approximated solution can be computed. The simulation results revealed that the proposed approach is effective for the goal at hand. Furthermore, because of its ease of use, our technique is relevant to a wide range of fractional problems. In the future, we can concentrate on the usage of the spectral Galerkin and the tau approaches for solving more complicated pricing models, such as the tempered FBSM. Finally, we should mention that the maximum absolute error for a given boundary value problem with a smooth solution is exponentially convergent. For nonsmoothness in time (or in space), the method’s order of convergence degrades. This, however, can be mitigated by employing the fractional-order Jacobi functions described herein.



---

*Source: 9365292-2022-07-20.xml*