# Numerical Method for American Option Pricing under the Time-Fractional Black–Scholes Model

Authors: Yesen Sun; Wenxiu Gong; Hongliang Dai; Long Yuan
Journal: Mathematical Problems in Engineering
 (2023)

---

## Abstract

The fractional Black–Scholes model has had limited applications in financial markets. Instead, the time-fractional Black–Scholes equation has attracted much research interest. However, it is difficult to obtain the analytic expression for American option pricing under the time-fractional Black–Scholes model. This paper will present an operator-splitting method to price the American options under the time-fractional Black–Scholes model. The fractional partial differential complementarity problem (FPDCP) that the American option price satisfied is split into two subproblems: a linear boundary value problem and an algebraic system. A high-order compact (HOC) scheme and a grid stretching (GS) method are considered for the linear boundary problem. Furthermore, numerical results show that the HOC scheme with a GS method gives an accurate numerical solution for American options under the time-fractional Black–Scholes model.


---

## Body

## 1. Introduction

A financial option is an important hedging instrument, which is one of the most popular and important financial derivatives. It is very important for financial researchers and practitioners to give a fair and reasonable option price. In the field of finance, one of the most classical option pricing models is the Black–Scholes model [1]. The basic assumptions are that the underlying asset price follows geometric Brownian motion and that volatility is a constant. However, the implied volatility depends on strike and maturity. Hence, the assumptions are inconsistent with the behavior of the underlying assets.A lot of models have been introduced as improvements to the Black–Scholes model. Fractional calculus is one possible way to explain longer term memory effects of stock price variations. With the discovery of the fractional structure of financial markets, geometric Brownian motion has been replaced by fractional Brownian motion which provides a powerful tool to explain longer memory effects of underlying asset prices. Some fractional Black–Scholes models have been proposed for options pricing. Wyss [2] priced a European option by the time-fractional Black–Scholes equation and gave a closed-form solution. Cartea and del-Castillo-Negrete [3] derived the FMLS, KOBol, and CGMY models for pricing exotic options with jumps. Jumarie [4] obtained the time-space fractional Black–Scholes equation based on the fractional Taylor formula and Ito lemma. Liang et al. [5] further introduced a single parameter and a biparameter fractional Black–Scholes model. Chen et al. [6] provided a simplified fractional Black–Scholes version of a model by Liang et al. Farhadi et al. [7] also proposed a time fractional Black–Scholes equation by considering the effect of trend memory.With the applications of the time fractional Black–Scholes equation in financial theory, different methods have been proposed for the solution of the time-fractional Black–Scholes equation. There are many researchers who present the analytical solution of the fractional Black–Scholes equation, see [2, 4–6, 8, 9]. The common method is Laplace transform coupled with other methods, such as homotopy perturbation methods and homotopy analysis methods or Laplace decomposition methods, see [9]. However, it is difficult to compute the solutions obtained by these methods. Therefore, it is necessary to develop some numerical methods for the fractional Black–Scholes equation.Song and Wang [10] developed an implicit finite difference scheme for pricing options under Jumarie’s model. Zhang et al. [11] developed a θ finite difference scheme to the time-fractional Black–Scholes equation. Zhang et al. [12] constructed a discrete implicit numerical scheme for the time-fractional Black–Scholes equation governing European options. However, all of the above numerical methods focus on European option pricing. There are only a few studies on the numerical pricing methods of American options under the time-fractional Black–Scholes equation. It is well known that American option pricing problems can be formulated as a linear complementarity problem (LCP) under geometric Brownian motion. The used methods for the LCP are the penalty method [13, 14], projected SOR method [15, 16], projected LU (PLU) method [17], and operator-splitting method [18–21]. Chen et al. [22] used a predictor-corrector method for pricing American options under a finite moment log-stable model. Chen et al. [23] proposed an operator-splitting method for American options under time fractional Black-Scholes equation.In this paper, we consider the operator-splitting method for the time-fractional Black–Scholes equation that the American options satisfied. The fractional partial differential complementarity problem (FPDCP) is split into two sub-problems: a linear boundary value problem and an algebraic system. We consider a high-order compact (HOC) scheme with a grid stretching (GS) method for the linear boundary problem. An HOC scheme is widely used in option pricing because of its high order convergence. Tangman et al. [24] discussed the applications of the HOC scheme for pricing American options and European options with stochastic volatility. An HOC finite difference scheme with a uniform mesh that is fourth order accurate in space and second order accurate in time is introduced for option pricing models with stochastic volatility [25]. Düring et al. [26] derived HOC finite difference schemes for option pricing in stochastic volatility models on nonuniform grids. Moreover, HOC finite difference schemes are used to discretize a nonlinear Black–Scholes equation which models transaction costs arising in the hedging of portfolios [27] and a class of parabolic partial differential equations with time and space-dependent coefficients as well as with mixed-second order derivative terms in n spatial dimensions [28]. The HOC finite difference scheme with a grid stretching method is also considered for European options under the time-fractional Black–Scholes equation (29). At last, we also show that the HOC scheme with a GS method gives an accurate numerical solution for American options under the time-fractional Black–Scholes model.The paper is organized as follows: In Section2, we describe some models for the time-fractional Black–Scholes equation. In Section 3, we first present the GS method for the fractional Black–Scholes equation. Then, we adopt the operator-splitting method for the time-fractional Black–Scholes equation that the American options satisfied. An efficient HOC scheme is constructed at last. In Section 4, the numerical experiments of the convergence and other numerical studies are shown. The conclusions are given in Section 5.


## 2. Model Description

The assumption of the Black–Scholes option pricing is that the stock priceSt follows the geometric Brownian motion.(1)dSt=μStdt+σStdBt,where μ and σ are constants and Bt is the standard Brownian motion. Based on the above assumption, Black and Scholes derived that the price of a European option VS,t satisfies the partial differential equation as(2)∂VS,t∂t+12σ2S2∂2VS,t∂S2+rS∂VS,t∂S−rVS,t=0.With discovery of the fractional structure of financial markets, the fractional Black–Scholes partial differential equations are introduced to extend financial theory.Caputo’s fractional derivative [30] is defined as(3)∂αVS,t∂tα=1Γ1−α∫0tt−ζ−αVζS,ζdζ,0<α<1.The left-modified Riemann–Liouville derivative is defined by(4)∂LαVS,t∂tα=1Γ1−αddt∫0tVS,ζ−VS,0t−ζαdζ,0<α<1.The right-modified Riemann–Liouville derivative is defined as follows:(5)∂RαVS,t∂tα=1Γ1−αddt∫tTVS,ζ−VS,Tζ−tαdζ,0<α<1.Wyss [2] used the right-modified Riemann–Liouville derivative deriving the time-fractional Black–Scholes model.(6)∂RαVS,t∂tα=−12σ2S2∂2VS,t∂S2−rS∂VS,t∂S+rVS,t,0<α≤1,V0,t=pt,V+∞,t=qt,VS,T=gS.Jumarie [4] took account of the stochastic processes of a time-fractional order. The stock value St satisfies the fractional differential equation as(7)dSt=μStdt+σStwtdtα/2,0<α≤1,where wt=t1−α/1−α! is a normalized Gaussian white noise, i.e., the zero mean and the unit variance. Using the equation and combining with the Ito Lemma and the fractional Taylor’s series of the price of the stock option, Jumarie [31] derived the fractional Black–Scholes equation. Jumarie’s model is given by(8)∂LαVS,t∂tα=−Γ1+α2σ2S2∂2VS,t∂S2−t1−αΓ2−αrS∂VS,t∂S+t1−αΓ2−αrVS,t,V0,t=pt,V+∞,t=qt,VS,T=gS.Farhadi et al. [7] used the left-modified Riemann–Liouville derivative deriving another time-fractional Black–Scholes model by considering the fractional trend effect. Farhadi’s model is obtained by the following equation:(9)∂LαVS,t∂tα=−12Γ21+ασ2S2∂2VS,t∂S2−rS∂VS,t∂S+rVt1−αΓ2−α,V0,t=pt,V+∞,t=qt,VS,T=gS.In order to convert these three models into initial value problems, we use variable substitutionτ=T−t.(10)∂RαVS,t∂tα=1Γ1−αddt∫tTVζS,ζ−VζS,Tζ−tαdζ=1Γ1−α−ddτ∫T−τTVζS,ζ−VζS,Tζ−T−ταdζ=−1Γ1−αddτ∫0τVζS,T−t−VζS,0τ−tαdt=−∂LαVS,τ∂τα.According to [11], Caputo’s derivative is consistent with the left-modified Riemann–Liouville derivative. Hence, we replace ∂LαVS,τ/∂tα by ∂αVS,τ/∂tα.Wyss’ model can be converted into(11)∂αV∂τα=12σ2S2∂2V∂S2+rS∂V∂S−rV.Similarly, taking the same substitution and using the following chain rule (Corollary 2.3in [4]),(12)fxαux=Γ2−αuα−1xfuαuxuxαx.We have(13)Vτα=−T−τα−1τ1−αVtα.Now, (8) and (9) can be transformed as follows:(14)T−τ1−ατα−1∂αV∂τα=Γ1+α2σ2S2∂2V∂S2+T−τ1−αΓ2−αrS∂V∂S−T−τ1−αΓ2−αrV,τα−1∂αV∂τα=12Γ21+ασ2S2∂2V∂S2+rS∂V∂S−rV1Γ2−α.The initial and boundary value conditions become(15)V0,τ=pτ,V+∞,τ=qτ,VS,0=gS.These European option price models are linear initial boundary value problems which can be easily solved. In our paper, we consider the problem of pricing American options under Wyss’ model. Jumarie’s model and Farhadi’s model can be found in Appendixes A and B.


## 3. Numerical Method

Using the change of the variablex=lnS, (11) can be rewritten as(16)∂αV∂τα=LV,V−∞,τ=pτ,V+∞,τ=qτ,Vx,0=gx,where(17)LV=12σ2∂2V∂x2+r−12σ2∂V∂x−rV.The solution region isΣ0=−∞≤x≤+∞,0≤τ≤T. To solve (16) numerically, we truncate the unbounded domain Σ0 into a finite interval Σ1=xmin≤x≤xmax,0≤τ≤T. We use the GS method at the strike price which can decrease the error and yield a better convergence. Let y∈0,1 be the transformed coordinate, then(18)x=fy=1ξsinhc2y+c11−y+x0,where c1=sinh−1ξxmin−x0, c2=sinh−1ξxmax−x0, and ξ is a stretching parameter. Jacobian and Hessian of fy are(19)Jy=c2−c1ξcoshc2y+c11−y,Hy=c2−c12ξsinhc2y+c11−y.Then, we have(20)LV=ωy∂2V∂y2+αy∂V∂y−rV,where(21)ωy=σ22Jy2,αy=2r−σ22Jy−σ2Hy2Jy3.This is the time-fractional Black–Scholes model for European options. Unlike the European options, the American option pricing problem is a free boundary problem. We consider the following linear complementarity problem for American put options. Assuming thatUy,τ is the price of American put options, we have(22)∂αU∂τα≥LU,U≥g,∂αU∂τα-LUU−g=0.The initial and boundary value conditions of the FPCDP are(23)Uy,0=gy=maxK−e1/ξsinhc2y+c11−y+x0,0,U0,τ=K,U1,τ=0.In this section, we present the numerical method for solving (8). Since closed-form solutions for (8) do not exist, a numerical method using a high-order compact scheme with operator splitting is developed for its solution.Assuming thatφy,τ is a slack function, FPCDP (8) becomes(24)∂αU∂τα-LU=φ,φ≥0,U≥g,φU−g=0.Letτj=j∆τ, j=0,1,2,⋯,N, yi=ih, and i=0,1,2,⋯,M, where Δτ=T/N and h=1/M are the time and space steps, respectively. According to [12], ∂αU/∂τα can be approximated by the following difference formula:(25)∂αUx,τj∂τα=ρ∑n=0j−1ωnUj−n−Uj−n−1+OΔτ2−α=ρω0Uj−Uj−1+ρ∑n=1j−1ωnUj−n−Uj−n−1+OΔτ2−α,where(26)ρ=1Γ2−αΔτα,ωn=n+11−α−n1−α.Then, we introduce the operator-splitting method for (24). Here, we introduce an intermediate variable U˜, so (24) will be split into two simple equations. First, we solve the following boundary value problem:(27)ρω0U˜j−Uj−1+∑n=1j−1ωnUj−n−Uj−n−1=LU˜j+φj−1,U˜j0=K,U˜j1=0.From (27), we can get the intermediate variable U˜j. Second, we solve the following algebraic equation:(28)ρω0Uj−U˜j=φj−φj−1,φj≥0,Uj≥g,φjUj−g=0.We have the following solution:(29)φj,Uj=0,U˜j−φj−1ρω0,φj−1≤ρω0U˜j−g,g−U˜jρω0+φj−1,g,otherwise.To verify the effectiveness of our method, we add (27) and (28) together and get(30)ρ∑n=0j−1ωnUj−n−Uj−n−1=LU˜j+φj,φj≥0,Uj≥g,φjUj−g=0.For (27), we use the HOC scheme. We rewrite (27) as(31)ωyi∂2U˜ij∂y2+αyi∂U˜ij∂y=fij,where(32)fij=ρω0U˜ij−Uij−1+∑n=1j−1ωnUij−n−Uij−n−1+rU˜ij−φij−1.Let(33)δyU˜ij=U˜i+1j−U˜i−1j2h,δy2U˜ij=U˜i+1j−2U˜ij+U˜i−1jh2.Taking Taylor expansion toU˜i+1j and U˜i−1j, we have(34)U˜i+1j=U˜ij+h∂U˜ij∂y+h22!∂2U˜ij∂y2+h33!∂3U˜ij∂y3+h44!∂4U˜ij∂y4+h55!∂5U˜ij∂y5+Oh6,U˜i−1j=U˜ij−h∂U˜ij∂y+h22!∂2U˜ij∂y2−h33!∂3U˜ij∂y3+h44!∂4U˜ij∂y4−h55!∂5U˜ij∂y5+Oh6.Hence,(35)U˜i+1j−2U˜ij+U˜i−1jh2=∂2U˜ij∂y2+h212∂4U˜ij∂y4+Oh4,U˜i+1j−U˜i−1j2h=∂U˜ij∂y+h26∂3U˜ij∂y3+Oh4,and(36)∂2U˜ij∂y2=δy2U˜ij−h212∂4U˜ij∂y4+Oh4,∂U˜ij∂y=δyU˜ij−h26∂3U˜ij∂y3+Oh4.Substituting (13) and (14) into equation (12), we have(37)ωyiδy2U˜ij+αyiδyU˜ij=fij+Rij+Oh4,where(38)Rij=ωyih212∂4U˜ij∂y4+αyih26∂3U˜ij∂y3.We can derive∂4U˜ij/∂y4 and ∂3U˜ij/∂y3 by differentiating (12) with respect to y(39)∂3U˜ij∂y3=b1yi∂2U˜ij∂y2+b2yi∂U˜ij∂y+b3yi∂fij∂y,∂4U˜ij∂y4=d1yi∂2U˜ij∂y2+d2yi∂U˜ij∂y+d3yi∂2fij∂y2+d4yi∂fij∂y,where(40)b1y=−ωy′+αyωy,b2y=−αy′ωy,b3y=1ωy,d1y=2ωy′2+3ωy′αy+αy2−ωyωy″−2αy′ωyωy2,d2y=2ωy′αy′+αyαy′−αy″ωyωy2,d3y=1ωy,d4y=−2ωy′+αyωy2.Then,Rij can be given as follows:(41)Rij=∂2U˜ij∂y2ωyih212d1yi+αyih26b1yi+∂U˜ij∂yωyih212d2yi+αyih26b2yi+∂2fij∂y2ωyih212d3yi+∂fij∂yωyih212d4yi+αyih26b3yi.Hence, (37) can be written as(42)AyU˜ij=Byfij+Oh4,where(43)Ay=e1yiδy2+e2yiδy,By=I+e3yiδy2+e4yiδy,e1yi=ωyi−ωyih212d1yi−αyih26b1yi,e2yi=αyi−ωyih212d2yi−αyih26b2yi,e3yi=ωyih212d3yi,e4yi=ωyih212d4yi+αyih26b3yi.Using the Crank–Nicolson difference scheme for (42), we obtain the HOC scheme(44)AyU˜ij+U˜ij+12=Byρω0U˜ij+1−Uij+∑n=1jωnUij−n+1−Uij−n+rU˜ij+U˜ij+12−φij,where(45)AyU˜ij=e1yiδy2U˜ij+e2yiδyU˜ij=p1yiU˜i−1j+p2yiU˜ij+p3yiU˜i+1j,p1yi=e1yih2−e2yi2h,p2yi=−2e1yih2,p3yi=e1yih2+e2yi2h,ByU˜ij=U˜ij+e3yiδy2U˜ij+e4yiδyU˜ij=q1yiU˜i−1j+q2yiU˜ij+q3yiU˜i+1j,q1yi=e3yih2−e4yi2h,q2yi=1−2e3yih2,q3yi=e3yih2+e4yi2h.Let(46)P=p2y1p3y1p1y2p2y2p3y2⋱⋱⋱p1yM−2p2yM−2p3yM−2p1yM−1p2yM−1M−1×M−1,Q=q2y1q3y1q1y2q2y2q3y2⋱⋱⋱q1yM−2q2yM−2q3yM−2q1yM−1q2yM−1M−1×M−1.Then, (44) becomes(47)G1U˜j+1=G2U˜j+B,B=2Qρ∑n=1jωn−ωn−1Uj−n+1−ωjU0−2Qφj+C,C=−g˜1U˜Mj+1eM−1+g˜2U˜MjeM−1−g˜3U˜0j+1e1+g˜4U˜0je1,where(48)G1=P−r+2ρω0Q,G2=−P+rQ,g˜1=p3y−2q3yρω0−q3yr,g˜2=−p3y+q3yr,g˜3=p1y−2q1yρω0−q1yr,g˜4=−p1y+q1yr,e1=1,0,…,0M−1′,eM−1=0,…,0,1M−1′.


## 4. Numerical Experiments

The numerical performance of the proposed method is described in this section. We first display the nodes distribution with the GS method and compare our option pricing method with other methods. Then, we present the effect ofα to option prices. Last, we give the order of convergence.The parameters are(49)S0=30,K=50,r=0.01,T=1,σ=0.1,Smax=2K,ξ=10lnK,M=80,N=80.We first display the nodes distribution with the GS method by (18). The result is shown in Figure 1. In Figure 1, we can see that the nodes are not evenly distributed. Figure 2 shows the American put option price with α=0.6 and α=0.4. Table 1 presents the comparison between the operator-splitting (OS) method based on HOC-GS, projected LU method, and operator splitting [23] of American option prices. The results in Table 1 show the effectiveness of our method.Figure 1
Nodes distribution with the grid stretching parameterξ=10/lnK.Figure 2
American put option price withα=0.6 and α=0.4.
(a)(b)Table 1
The comparison between different methods of American option prices.
MethodsGrid (N, M)αS0=30S0=40S0=50OS (HOC-GS)28,280.420.000010.00001.67200.620.000010.00001.7198PLU28,280.420.000010.00001.67130.620.000010.00001.7217OS28,280.420.000010.00001.67470.620.000010.00001.7232Figure 3
American option prices at differentα.Now, we consider the effect of the parameterα. The parameters are set as S0=45,K=50,r=0.01,T=1,σ=0.4,Smax=2K,ξ=10/lnK, and M=200,N=100. Figure 3 shows the American option price at different α. Table 2 shows the computational time for varying α and M,N. Figure 4 shows the American option price at different α and stock price. Figures 4(a) and 4(b) provide that the American put option price increases in α at the at-the-money region and out-of-the-money region. From the two figures, we can see how α affects the American put option prices. In Figure 5, we compare the American put option price under Wyss’, Jumarie’s, and Farhadi’s models by our proposed method. In the figure, we can see the similar results between Wyss’ and Jumarie’s models and the pricing of Farhadi’s model is higher than that of the other two models.Table 2
The computational time for varyingα and M,N.
M,Nα0.10.40.50.71620.01820.01720.01850.01683240.01850.01810.01970.02156480.02390.02090.02320.0229128160.11620.10650.09630.1094256320.66110.67390.66220.6709512649.50179.54549.88409.62471024128195.4311181.0966181.579182.3388Figure 4
American option prices at differentα and stock prices.
(a)(b)Figure 5
American option prices under Wyss’, Jumarie’s, and Farhadi’s models.Next, we check the convergence rates of the HOC-GS scheme. The orders of convergence arelogΔt1/Δt2ek/ek/2 for the time step size ∆t and logh1/h2ek/ek/2 for the spatial step size h, where ek=Vk−Vk/2, Vk is the numerical solution at the final time step, and ∙ is l2 norm or l∞ norm.We check the convergence order for European options. The convergence order in time and space directions for European options is listed in Table3. Tables 4 and 5 show the orders of convergence in time and space for American put options. In Tables 4 and 5, we can see that convergence orders in time and spatial are oscillation. α has an effect on convergence orders. When α is bigger, the convergence orders are more stable and about 2−α which is better than the method in [23]. The reason of the oscillation is that the two subproblems are coupled which influence convergence orders of the HOC-GS scheme. In addition, the orders oscillate but still appear better than those of the method proposed in [23].Table 3
Numerical convergence order for European put options in time and space forα=0.9.
N (M = 50)OrderM (N = 50)Order1001.31791003.95782001.04602003.82654001.04514003.73658001.03968003.7661Table 4
Numerical convergence order in time for varyingα when M=320.
N, α0.10.40.7l2l∞l2l∞l2l∞200.98241.31561.10071.19461.27511.3337400.97230.96140.99140.811.13010.9804802.91.38460.9580.87451.06521.3026160−1.36250.58590.9411.07040.95231.3219Table 5
Numerical convergence order in space for varyingα when N=320.
M, α0.10.40.7l2l∞l2l∞l2l∞400.76780.5851.77471.89812.52822.523580−1.4272−12.05313.45942.19142.29951600.789911.86351.7372.23262.1155


## 5. Conclusion

This paper gives an efficient method to solve the time-fractional Black–Scholes equation for American option pricing. We first employ the operator-splitting method to solve FPDCP that the American option price satisfied. FPDCP is split into two subproblems: a linear boundary problem and an algebraic system. The two subproblems are easy to solve by numerical methods. For the linear boundary problem, we consider an HOC scheme and a GS method which have a higher order convergence. The numerical results demonstrate the effectiveness of our method. Also, we study the effect ofα on the option price.



---

*Source: 4669161-2023-02-20.xml*