# Lattice Boltzmann Method for the Generalized Black-Scholes Equation

Authors: Fangfang Wu; Duoduo Xu; Ming Tian; Yingying Wang
Journal: Advances in Mathematical Physics
 (2023)

---

## Abstract

In this paper, an efficient lattice Boltzmann model for the generalized Black-Scholes equation governing option pricing is proposed. The Black-Scholes equation is firstly equivalently transformed into an initial value problem for a partial differential equation with a source term using the variable substitution and the derivative rules, respectively. Then, applying the multiscale Chapman-Enskog expansion, the amending function is expanded to second order to recover the convective and source terms of the macroscopic equation. The D1Q3 lattice Boltzmann model with spatial second-order accuracy is constructed, and the accuracy of the established D1Q5 model is greater than second order. The numerical simulation results demonstrate the effectiveness and numerical accuracy of the proposed models and indicate that our proposed models are suitable for solving the Black-Scholes equation. The proposed lattice Boltzmann model can also be applied to a class of partial differential equations with variable coefficients and source terms.


---

## Body

## 1. Introduction

It is widely acknowledged among the researchers that the lattice Boltzmann method (LBM) has been developed into an alternative and promising numerical scheme for simulating fluid flows and modeling physics in fluids [1–4]. Compared with traditional numerical methods, LBM is a mesoscopic numerical method which is based on microscopic models and mesoscopic kinetic equations, and it has many outstanding advantages, such as its algorithmic simplicity, parallel computation, and easy handling of complicated boundary conditions. Recently, LBM has been further adopted to solve some partial differential equations, including convection-diffusion equation [5–7], Burgers’ equation [8–10], wave equation [11, 12], and partial differential equations with variable coefficients [13–15].With the rapid development of the global financial market, the risks of international commodity trading have gradually increased. Among the many financial derivatives, options, futures, and forward contracts are the most basic. Compared with other financial derivatives, options have greater flexibility and can create more investment opportunities for investors, a popular hedging tool for investors [16]. A traded option in a financial market is a contract which gives its owner the right to buy (call option) or sell (pull option) a fixed quantity of assets of a specified stock at a fixed price (exercise or strike price) on (European option) or before (American option) a given date (expiry date) [17]. However, the complex and stochastic nature of financial market makes it challenging to rationally determine the values of options. A huge improvement has been made in the option pricing theory since the development of the Black-Scholes option pricing model by Fischer Black and Myron Scholes in 1973 and previously by Robert Merton in Ref. [18, 19]. Subsequently, numerous studies proposed to use the generalized Black-Scholes model to describe the dynamic process of the underlying asset and solve the generalised Black-Scholes equation numerically. The generalized Black-Scholes equation (BSE) has the form [20, 21]
(1)∂V∂t~+rS,t~−dS,t~S∂V∂S+12σ2S,t~S2∂2V∂S2−rS,t~V=0,S,t~∈0,Smax×0,T,with the final condition
(2)VS,T=φS,S∈0,Smax,where φS is a smooth function and boundary conditions
(3)V0,t~=0,t~∈0,T,(4)VSmax,t~=Smaxexp−∫t~TdSmax,sds−Eexp−∫t~TrSmax,sds,where VS,t~ is the European option price at underlying asset price S and at time t~, rS,t~>0 is the risk-free interest rate, dS,t~ is the dividend, and σS,t~ is the volatility. rS,t~,dS,t~ and σS,t~ are all the functions of S and t~. E is the exercise price, T is the expiry date, and Smax is the suitably chosen positive number. The BSE can be solved exactly when the coefficients are constant or space-independent [20–22]. However, finding the closed form solution of the BSE is difficult when the coefficients are discontinuous ordinary functions or generalized functions. Therefore, it is essential to choose an efficient and accurate numerical algorithm to solve this type of problems.Ross et al. [23] were the first to apply the lattice approach to the BSE, and it was improved by Hull and White [24]. This method is equivalent to an explicit time-stepping scheme. The finite difference method is one of the most common numerical approaches to solve the BSE in order to avoid probabilistic and stochastic methods. For instance, Cen and Le [25] presented a robust and accurate finite difference method with spatial second-order accuracy for a generalized BSE based on a central difference spatial discretization on a piecewise uniform mesh and an implicit time-stepping technique. The method proposed by Rao [26] uses the high-order difference approximation with identity expansion (HODIE) scheme in the spatial direction and the two-step backward differentiation formula in the temporal direction, which has second-order convergence in space as well as in time. Roul and Goura [27] proposed a high-order compact finite difference method based on a uniform mesh to obtain a highly accurate result for a generalized BSE. In addition, Mohammadi [28] developed a collocation method with the quintic B-spline functions for the generalized BSE. Nyoumbi and Tambue [29] presented and analysed the finite volume method with a fitted technique for the two-dimensional BSE with stochastic volatility. In recent years, several methods have been proposed to solve the time-fractional BSE numerically such as radial basis function method by Golbabai et al. [30], localized kernel-based meshless method by Nikan et al. [31], and the moving least-squares (MLS) method by Golbabai and Nikan [32].Since most numerical methods in the literature are either too complicated to implement or with very low computational efficiency, it is desirable to have an alternative easily performed method for the option pricing. Up to now, rarely have researchers conducted numerical investigation on BSE with LBM. Therefore, this paper aims to propose an efficient lattice Boltzmann model for the generalized BSE. A second-order multiscale Chapman-Enskog expansion of the amending function is used to recover the convection and source terms in the macroscopic equations. The lattice Boltzmann models with second-order and higher second-order accuracy are established using the D1Q3 model and D1Q5 model, respectively. European call option, binary option, and butterfly spread option are simulated efficiently to test the performance of the proposed method. Naturally, the proposed lattice Boltzmann model can be adopted to solve American options if it is used together with a technique for free boundary problems.The remainder of this paper is as follows. In Section2, Equations (1) and (2) are transferred into the equivalent initial value problem through variable substitution and derivative rule in order to resolve the backward in time problem of the option pricing. In Section 3, the lattice Boltzmann model for the generalized BSE is described. The proposed model is verified by several numerical examples in Section 4. Finally, conclusions are summarized in Section 5.


## 2. Preliminaries

The BSE (1) is a degenerate parabolic partial differential equation with variable coefficients. Further, the partial differential equations given in Equations (1)–(4) are backward in time. In order to construct lattice Boltzmann model for this type of problems, using variable substitutions t~=T−t,S=ex and derivative rule, Equations (1)–(4) are equivalent to the following form:
(5)∂u∂t+∂Ax,tu∂x−∂2Bx,tu∂x2=Fx,t,u,with initial condition(6)ux,0=φexpx,x∈0,xmax,and boundary conditions(7)u0,t=0,t∈0,T,uxmax,t=expxmax−∫0td^xmax,sds−Eexp−∫0tr^xmax,sds,where
(8)ux,t=Vexpx,T−t,Ax,t=d^x,t−r^x,t+12σ∧2x,t+2σ^x,tσ^xx,t,Bx,t=12σ∧2x,t,Fx,t,u=d^xx,t−r^xx,t+σ^x2x,t+σ^x,tσ^xxx,t+σ^x,tσ^xx,t−r^x,tu,σ^x,t=σexpx,T−t,r^x,t=rexpx,T−t,d^x,t=dexpx,T−t.σ^xx,t, d^xx,t, and r^xx,t denote the first-order partial derivatives of σ^x,t, d^x,t, and r^x,t with respect to x, respectively. σ^xxx,t denotes the second-order partial derivative of σ^x,t with respect to x. Next, a lattice Boltzmann model is built for Equation (5).


## 3. Lattice Boltzmann Method

A lattice Boltzmann model applied in this paper is the Bhatnagar-Gross-Krook model with the D1Qb [3] (one dimension in space and b discrete directions) lattice. The evolution law of the particle distribution function in the model is written as
(9)fαx+eαΔt,t+Δt−fαx,t=−1τfαx,t−fαeqx,t+Δthαx,t,where fαx,t and fαeqx,t are the local particle distribution function and the local equilibrium distribution function at time t and position x in direction α, respectively, hαx,t is used for amending the convection and source terms in the recovered macroscopic equation, eα is a superficial phase-space variable, Δt is time step, and τ is the dimensionless relaxation time. The stability of the equation requires that τ>0.5 [33].To recover the macroscopic BSE correctly by using Equation (9), we assume that the equilibrium distribution function fαeqx,t and particle distribution function fαx,t satisfy local mass conservation law, i.e.,
(10)ux,t=∑αfαx,t=∑αfαeqx,t.Now, using the Taylor series expansion to Equation (9), we have
(11)Δt∂∂t+eα∂∂xfα+Δt22∂∂t+eα∂∂x2fα+Δt36∂∂t+eα∂∂x3fα=−1τfα−fαeq+Δthα.The multiscale Chapman-Enskog expansion up to first order in space and second order in time is applied; we have(12)∂∂x=ε∂∂x1,(13)∂∂t=ε∂∂t1+ε2∂∂t2+ε3∂∂t3,where ε is a small expansion parameter. Furthermore, fαx,t and hαx,t can be expressed as
(14)fα=fα0+εfα1+ε2fα2+⋯,(15)hα=εhα1+ε2hα2+⋯.Substituting Equations (12)–(15) into Equation (11), we obtain a series of partial differential equations in different orders of ε:
(16)Oε0:−1τfα0−fαeq=0,i.e.,fα0=fαeq,(17)Oε1:∂∂t1+eα∂∂x1fα0=−1τΔtfα1+hα1,(18)Oε2:∂∂t2fα0+Δt12−τ∂∂t1+eα∂∂x12fα0+Δtτ∂∂t1+eα∂∂x1hα1=−1τΔtfα2+hα2,(19)Oε3:∂∂t3fα0+Δt1−2τ∂∂t2∂∂t1+eα∂∂x1fα0+Δt2τ2−τ+16∂∂t1+eα∂∂x13fα0+Δt2τ2−τ2∂∂t1+eα∂∂x12hα1+Δtτ∂∂t2hα1+τΔt∂∂t1+eα∂∂x1hα2=−1τΔtfα3.From Equations (10), (14), and (16), we have
(20)∑αfα0=ux,t,∑αfαn=0,n>0.In order to recover Equation (5), some constraints on the higher moments of the local equilibrium distribution function are imposed as
(21)∑αeαfα0=0,(22)∑αeα2fα0=λ2Bx,tu,(23)∑αeα3fα0=0.Meanwhile, the amending function satisfies(24)∑αhα1=0,∑αhα2=1ε2Fx,t,u,(25)∑αeαhα1=λ1εAx,tu,∑αeαhα2=0,(26)∑αeα2hα1=0,∑αeα2hα2=0.Here,λ1 and λ2 are the parameters to be determined in the following. According to Equations (15), (24), (25), and (26), we have
(27)∑αhα=ε∑αhα1+ε2∑αhα2=Fx,t,u,(28)∑αeαhα=ε∑αeαhα1+ε2∑αeαhα2=ε∑αeαhα1=λ1Ax,tu,(29)∑αeα2hα=ε∑αeα2hα1+ε2∑αeα2hα2=0.Summing Equation (17) over α and using Equations (20), (21), and (24), we obtain
(30)∂u∂t1=∑αhα1.Similarly, summing Equation (18) over α and using Equation (20), we obtain
(31)∂u∂t2+Δt12−τ∑α∂∂t1+eα∂∂x12fα0+Δtτ∑α∂∂t1+eα∂∂x1hα1=∑αhα2.According to Equations (20)–(22), (24), and (30), we have
(32)∑α∂∂t1+eα∂∂x12fα0=∂2∂t12∑αfα0+2∂2∂t1∂x1∑αeαfα0+∂2∂x12∑αeα2fα0=λ2∂2∂x12Bx,tu.According to Equations (24) and (25), we have
(33)∑α∂∂t1+eα∂∂x1hα1=∂∂t1∑αhα1+∂∂x1∑αeαhα1=λ1ε∂∂x1Ax,tu.Substituting Equations (32) and (33) into Equation (31) yields
(34)∂u∂t2+Δtτλ1ε∂∂x1Ax,tu−Δtτ−12λ2∂2∂x12Bx,tu=∑αhα2.Summing Equation (19) over α and using Equation (20), we obtain
(35)∂u∂t3+Δt2τ2−τ+16∑α∂∂t1+eα∂∂x13fα0+Δt2τ2−τ2∑α∂∂t1+eα∂∂x12hα1+τΔt∑α∂∂t1+eα∂∂x1hα2=0.Using Equations (20)–(26) and (30), we have
(36)∑α∂∂t1+eα∂∂x13fα0=∂3∂t13∑αfα0+3∂3∂x1∂t12∑αeαfα0+3∂3∂x12∂t1∑αeα2fα0+∂3∂x13∑αeα3fα0=3λ2∂3∂x12∂t1Bx,tu,(37)∑α∂∂t1+eα∂∂x12hα1=∂2∂t12∑αhα1+2∂2∂x1∂t1∑αeαhα1+∂2∂x12∑αeα2hα1=2λ1ε∂2∂x1∂t1Ax,tu.Substituting Equations (36) and (37) into Equation (35) yields
(38)∂u∂t3+Δt22λ1ετ2−τ2∂2∂x1∂t1Ax,tu+3Δt2λ2τ2−τ+16∂3∂x12∂t1Bx,tu+τΔt1ε2∂∂t1Fx,t,u=0.


### 3.1. D1Q3 Model

For the D1Q3 model,e0,e1,e2=0,c,−c, where c=Δx/Δt with Δx representing grid spacing. The macroscopic equation with the second order of truncation error is derived by taking Equation (30)×ε+Equation (34)×ε2 and using Equation (27):
(39)∂u∂t+Δtτλ1∂∂xAx,tu−Δtτ−12λ2∂2∂x2Bx,tu+Oε3=Fx,t,u.To recover Equation (5), we let
(40)λ1=1Δtτ,λ2=1Δtτ−0.5.Then, Equation (39) becomes
(41)∂u∂t+∂∂xAx,tu−∂2∂x2Bx,tu+Oε3=Fx,t,u.Coupling Equations (20)–(22), we can derive the local equilibrium distribution function fα0α=0,1,2 as follows:
(42)f00=u−Bx,tuc2Δtτ−0.5,f10=f20=Bx,tu2c2Δtτ−0.5.And according to Equations (27) and (28), a simple case for the amending function hαα=0,1,2 is given as
(43)h0=Fx,t,u,h1=−h2=Ax,tu2cΔtτ.Here,τ is a free parameter to be determined. The numerical accuracy and stability can be improved by choosing appropriate τ.


### 3.2. D1Q5 Model

For the D1Q5 model,e0,e1,e2,e3,e4=0,c,−c,2c,−2c. The macroscopic equation is derived by taking Equation (30)×ε+Equation (34)×ε2+Equation (38)×ε3 and using Equations (27) and (40):
(44)∂u∂t+∂∂xAx,tu−∂2∂x2Bx,tu+Re+Oε4=0,where the truncation error is
(45)Re=2Δt2ε2λ1τ2−τ2∂2∂x∂t1Ax,tu+3Δt2ε3λ2τ2−τ+16∂3∂x2∂t1Bx,tu+τΔtε∂∂t1Fx,t,u.AsΔt and ε are approaching to zero, Equation (44) coincides with Equation (5). Although the D1Q5 model did not achieve third-order convergence in space due to the existence of Re, a convergence higher than second order was constructed in this model.Remark 1.
WhenAx,t=Ax, Bx,t=Bx, and Fx,t,u=Fx,u, using Equation (30), we have
(46)∂2∂x∂t1Axu=∂∂x∂∂t1Axu=∂∂xAx∂u∂t1=0,∂3∂x2∂t1Bxu=∂2∂x2∂∂t1Bxu=∂2∂x2Bx∂u∂t1=0.∂∂t1Fx,u=∂Fx,u∂u∂u∂t1=0
Then, the truncation errorRe=0.Coupling Equations (20)–(23), we can derive the equilibrium distribution function fα0α=0,1,2,3,4 as follows:
(47)f00=u−5Bx,tu4c2Δtτ−0.5,f10=f20=2Bx,tu3c2Δtτ−0.5,f30=f40=−Bx,tu24c2Δtτ−0.5.According to Equations (27)–(29), a simple case for the amending function hαα=0,1,2,3,4 is given as
(48)h0=Fx,t,u,h1=−h2=2Ax,tu3cΔtτ,h3=−h4=−Ax,tu12cΔtτ.


## 3.1. D1Q3 Model

For the D1Q3 model,e0,e1,e2=0,c,−c, where c=Δx/Δt with Δx representing grid spacing. The macroscopic equation with the second order of truncation error is derived by taking Equation (30)×ε+Equation (34)×ε2 and using Equation (27):
(39)∂u∂t+Δtτλ1∂∂xAx,tu−Δtτ−12λ2∂2∂x2Bx,tu+Oε3=Fx,t,u.To recover Equation (5), we let
(40)λ1=1Δtτ,λ2=1Δtτ−0.5.Then, Equation (39) becomes
(41)∂u∂t+∂∂xAx,tu−∂2∂x2Bx,tu+Oε3=Fx,t,u.Coupling Equations (20)–(22), we can derive the local equilibrium distribution function fα0α=0,1,2 as follows:
(42)f00=u−Bx,tuc2Δtτ−0.5,f10=f20=Bx,tu2c2Δtτ−0.5.And according to Equations (27) and (28), a simple case for the amending function hαα=0,1,2 is given as
(43)h0=Fx,t,u,h1=−h2=Ax,tu2cΔtτ.Here,τ is a free parameter to be determined. The numerical accuracy and stability can be improved by choosing appropriate τ.


## 3.2. D1Q5 Model

For the D1Q5 model,e0,e1,e2,e3,e4=0,c,−c,2c,−2c. The macroscopic equation is derived by taking Equation (30)×ε+Equation (34)×ε2+Equation (38)×ε3 and using Equations (27) and (40):
(44)∂u∂t+∂∂xAx,tu−∂2∂x2Bx,tu+Re+Oε4=0,where the truncation error is
(45)Re=2Δt2ε2λ1τ2−τ2∂2∂x∂t1Ax,tu+3Δt2ε3λ2τ2−τ+16∂3∂x2∂t1Bx,tu+τΔtε∂∂t1Fx,t,u.AsΔt and ε are approaching to zero, Equation (44) coincides with Equation (5). Although the D1Q5 model did not achieve third-order convergence in space due to the existence of Re, a convergence higher than second order was constructed in this model.Remark 1.
WhenAx,t=Ax, Bx,t=Bx, and Fx,t,u=Fx,u, using Equation (30), we have
(46)∂2∂x∂t1Axu=∂∂x∂∂t1Axu=∂∂xAx∂u∂t1=0,∂3∂x2∂t1Bxu=∂2∂x2∂∂t1Bxu=∂2∂x2Bx∂u∂t1=0.∂∂t1Fx,u=∂Fx,u∂u∂u∂t1=0
Then, the truncation errorRe=0.Coupling Equations (20)–(23), we can derive the equilibrium distribution function fα0α=0,1,2,3,4 as follows:
(47)f00=u−5Bx,tu4c2Δtτ−0.5,f10=f20=2Bx,tu3c2Δtτ−0.5,f30=f40=−Bx,tu24c2Δtτ−0.5.According to Equations (27)–(29), a simple case for the amending function hαα=0,1,2,3,4 is given as
(48)h0=Fx,t,u,h1=−h2=2Ax,tu3cΔtτ,h3=−h4=−Ax,tu12cΔtτ.


## 4. Numerical Simulations

To show the efficiency of the models proposed above, we simulate the value of various European options with different final and boundary conditions and choices of parameters. The distribution functionsfαx,t are initialized by setting to be fαeqx,t. The nonequilibrium extrapolation scheme proposed by Guo et al. [34] is employed to treat the boundary conditions. In addition, to test the precision of the present models, the maximum absolute error (MAE), the root mean squar error (RMS), and corresponding orders of convergence PMAE and PRMS are defined as follows:
(49)EMAEM,N=max0≤i≤Muxi,t−u∗xi,t,ERMSM,N=∑0≤i≤Muxi,tN−u∗xi,tN2M+1,PMAEM,N=log2EMAEM,NEMAE2M,2N,PRMSM,N=log2ERMSM,NERMS2M,2N,where uxi,t and u∗xi,t are the numerical solution and analytical solution on the M grids in space and N grids in time, respectively. Computations for the problems included in this paper were done by C-Free 5.0 software on an 11th Intel Core CPU 2.40 GHz PC equipped with 16GB RAM.Example 1.
Consider the following general Black-Scholes type equation [35]:
(50)∂u∂t=αx,tuxx+βx,tux+γx,tu+fx,t,x,t∈−2,2×0,1,with
(51)αx,t=0.082+T−tsinexpx2,βx,t=0.061+texp−expx−0.02exp−t−expx−αx,t,γx,t=−0.061+texp−expx,where fx,t is chosen such a way that ux,t=expx−t.In the computation, we takeΔt=10−5N=105, Δx=0.0125M=320, τ=0.9 for the D1Q3 model, and τ=1.01 for the D1Q5 model. The computational domain is fixed on −2,2. Table 1 compares the numerical errors of D1Q3 and D1Q5 models in different times. Figure 1 plots the time evolution of the numerical solution for D1Q5 model and analytical solution.Table 1
Errors for Example1 at different times.
t0.20.40.60.81.0MAED1Q36.9826e-044.4809e-042.6458e-041.3020e-044.2636e-05RMSD1Q31.4863e-041.0714e-047.1613e-054.1477e-051.9956e-05MAED1Q54.2793e-052.9890e-051.6116e-057.1896e-062.9742e-06RMSD1Q59.3560e-067.5258e-064.8232e-062.7756e-061.7843e-06Figure 1
Time evolution of the value of the option for Example1. (a) Numerical solution versus analytical solution at t=0.0,0.2,0.4,0.8,1.0. (b) Time evolution of numerical solution.
(a)(b)To test the accuracy of the D1Q3 model, simulations are performed for the fixed time stepΔt=10−5N=105 and Δt=10−6N=106, respectively. The grid resolutions range from M=50 to M=800 (Δx=0.005 to Δx=0.04). The dimensionless relaxation time τ=0.9 is fixed for the computations. Figure 2 shows the log-log plots of the RMS and MAE versus grid spacing at t=1. It can be seen that the proposed D1Q3 model for BSE is of the second-order accuracy in space (the slopes of the lines are all about 2.0), which is consistent with the theoretical accuracy.Example 2.
Consider Equation (1) for European call option with the final condition VS,T=S−E+ and boundary conditions V0,t~=0 and VS,t~=Sexp−dT−t~−Eexp−rT−t~ as S⟶∞.(a)
In this case, the analytical solution of Equation (1) with σ=0.4, r=0.06, d=0.02, E=1.0, and T=1 is(52)VS,t~=SNd1exp−dT−t~−ENd2exp−rT−t~,Nx=12π∫−∞xexp−12y2dy,d1S,t~=logS/E+r−d+1/2σ2T−t~σT−t~,d2S,t~=d1S,t~−σT−t~
In the computation, we takeτ=0.98, Δt=10−4N=104, and Δx=0.05M=80 for the D1Q5 model. The computational domain is fixed on x∈−2,2. Figure 3 shows the time evolution of the numerical solution and analytical solution. The option pricing is calculated backward in time, which means, the solution at the moment t=0 in Figure 3 is the desired European option price.
To identify the accuracy of the proposed models, simulations are performed for the different grid sizes, time steps, and dimensionless relaxation times, respectively. The grid resolutions fromM=40 to M=640 (Δx=0.0125 to Δx=0.2) are repeatedly calculated. Table 2 includes the numerical errors, orders of convergence, and CPU times (in seconds) obtained from the D1Q3 and D1Q5 models, respectively, when the fixed time step Δt=10−4N=104. From Table 2, it can be seen that the order of convergence of the D1Q3 model is close to 2.0, while that of the D1Q5 model is between 2.0 and 3.0, which is consistent with the theory, and the computational results and computational efficiency of the D1Q5 model are better than those of the D1Q3 model. Figure 4 depicts the log-log plots of the MAE versus the lattice spacing at t=0 for different dimensionless relaxation times and time steps. From Figure 4, it can be found that the slopes of six lines are all close to 2.0, which also confirms that the D1Q3 model is of second-order accuracy in space.
(b)
Consider Equation (1) with σ=0.4, r=0.04, d=0.02, E=1.0, and T=1, and the analytical solution for this case is the same as in Example 2(a)
The computational domain is fixed onx∈−2,2. Tables 3 and 4 give the errors, orders of convergence, and CPU times for the D1Q5 model and HODIE scheme [26] for the BSE, respectively. By comparing the above results, it is shown that our proposed algorithm is more stable and efficient than the HODIE scheme.
(c)
Consider Equation (1) with σS,t~=0.4, rS,t~=0.1+0.02sin10Tt~, dS,t~=0.06S/X, E=400, X=700, and T=1. In this case, the analytical solution is not known
In the computation, we takeτ=1.0, Δt=10−5N=105, Smax=700, and M=100. Figure 5 plots the time evolution of numerical solution, which agrees well with the result in Ref. [36].Example 3.
Consider Equation (1) for the binary European call option with the final condition
(53)VS,T=Q,S≥E,0,S<E,and the boundary conditions
(54)V0,t~=0,VS,t~=Qexp−rT−t~,S⟶∞.(a)
In this case, the analytical solution of Equation (1) with σS,t~=0.4, rS,t~=0.06, dS,t~=0.02, E=1, Q=1, and T=1 is(55)VS,t~=QNd2exp−rT−t~,Nx=12π∫−∞xexp−12y2dy,d1S,t~=logS/E+r−d+1/2σ2T−t~σT−t~,d2S,t~=d1S,t~−σT−t~
In the computation, we takeτ=1.0, Δt=10−5N=105, and Δx=0.025M=160 for the D1Q5 model. The computational domain is fixed on x∈−2,2. Figure 6(a) presents the comparisons between the numerical solutions and analytical solutions at different times, showing that the numerical solutions are in good agreement with the analytical solutions. Figure 6(b) plots the time evolution of numerical solutions.
(b)
Consider Equation (1) with σS,t~=0.42+t~sinS, rS,t~=0.061+T−t~exp−S, dS,t~=0.02expt~−T−S, E=4, Q=2, and T=1. In this case the analytical solution is not known
In the computation, we takeτ=1.0, Δt=10−5N=105, and Δx=0.025M=160. The computational domain is fixed on x∈−2,2. Figure 7 shows the time evolution of the option price. It can be found that the numerical solutions from our method are smooth for all t∈0,1 and consistent with Ref. [34].Example 4.
Here, we consider butterfly spread option which is a combination of three options with three strike prices, in which one contract is purchased with two outside strike prices and two contracts are sold at middle strike price [35].
(a)
Consider Equation (1) for butterfly spread option with σS,t~=0.42+t~exp−S, rS,t~=0.061+T−t~exp−S, dS,t~=0.02expt~−T−S, E1=1, E2=2, E3=3, and T=1. The final condition and boundary conditions for the butterfly spread option are as follows:(56)VS,T=maxS−E1,0−2maxS−E2,0+maxS−E3,0,V0,t~=0,VSmax,t~=0
In this case, the analytical solution is not known.
In the computation, we takeτ=1.0, Δt=10−5N=105, and Δx=0.025M=160. The computational domain is fixed on x∈−2,2. Figure 8 shows the time evolution of the value of the options, which is consistent with the results in Ref. [35].
(b)
Consider Equation (1) with σS,t~=0.4, rS,t~=0.1+0.02sin10Tt~, dS,t~=0.06S/E, E1=40, E2=50, and E3=60. The final condition for the butterfly spread option is as follows:
(57)VS,T=1,S∈E1,E2,−1,S∈E2,E3,0,otherwise
In this case, the analytical solution is not known also.
In the computation, we takeτ=1.0, Δt=10−5N=105, Smax=100, and M=100. Figure 9 plots the time evolution of the value of the option, which is consistent with the results in Ref. [36].Figure 2
Accuracy test for Example1 at t=1 with different time steps.Figure 3
Time evolution of the value of the option for Example2(a). (a) Numerical solution versus analytical solution at t=1.0,0.0. (b) Time evolution of numerical solutions.
(a)(b)Table 2
Comparisons of the D1Q3 and D1Q5 models for Example2(a) with the fixed time step ∆t=10−4.
MD1Q3D1Q5RMSD1Q3PRMSMAED1Q3PMAECPU time (s)RMSD1Q5PRMSMAED1Q5PMAECPU time (s)404.4440e-04—1.0994e-03—0.1013.4386e-04—7.7946e-04—0.015801.1268e-041.982.7161e-042.010.2198.6980e-051.981.9321e-042.010.0471602.7967e-052.016.5878e-052.040.4232.1467e-052.014.6360e-052.060.1173206.0860e-062.201.5139e-052.120.8634.6181e-062.229.7959e-062.240.2356401.2328e-062.303.7704e-062.011.6958.3944e-072.461.5778e-062.630.455Figure 4
Accuracy test for Example2(a) at t=0.Table 3
Errors and corresponding orders of convergence of the D1Q5 model for Example2(b).
M242526272829210211N2−7×1052−6×1052−5×1052−4×1052−3×1052−2×1052−1×105105MAE4.8024e-031.2180e-033.0455e-047.4115e-051.7551e-053.9023e-067.6816e-071.5426e-07PMAE—1.97921.99982.03892.07822.16922.34492.3160RMS2.0846e-035.1770e-041.2885e-043.1874e-057.7768e-061.8462e-064.1112e-078.5616e-08PRMS—2.00962.00642.01532.03512.07462.16702.2636CPU time (s)0.0000.0000.0150.0620.23500.9173.71315.271Table 4
Errors and corresponding orders of convergence of the HODIE scheme [26] for Example 2(b).
M242526272829210N1010×210×2210×2310×2410×2510×26MAE3.4226e-027.4769e-031.7759e-034.4895e-041.1219e-042.8068e-057.0223e-06PMAE—2.19462.07391.98392.00061.99891.9989RMS8.4476e-031.8556e-034.5900e-041.1499e-042.8825e-057.2183e-061.8061e-06PRMS—2.23042.01531.99701.99621.99751.9987Figure 5
Time evolution of the numerical solution for Example2(c).Figure 6
Time evolution of the value of the option for Example3(a). (a) Numerical solution versus analytical solution at t=1.0,0.4,0.2,0.0. (b) Time evolution of numerical solutions.
(a)(b)Figure 7
Time evolution of the value of the option for Example3(b). (a) Numerical solution at t=1.0,0.8,0.4,0.2,0.0. (b) Time evolution of numerical solution.
(a)(b)Figure 8
Time evolution of the value of the option for Example4(a). (a) Numerical solution at t=1.0,0.8,0.4,0.2,0.0. (b) Time evolution of numerical solutions.
(a)(b)Figure 9
Time evolution of the value of the option for Example4(b). (a) Numerical solutions at different times. (b) Time evolution of numerical solutions.
(a)(b)


## 5. Conclusions

In this paper, we presented a lattice Boltzmann model with an amending function for the numerical solution of the generalized BSE governing European options. In order to construct lattice Boltzmann model, the BSE is equivalently transformed into an initial value problem for a partial differential equation with a source term using the variable substitution and derivative rules, respectively. Then, the amending functionhαx,t is expanded to second order using the multiscale Chapman-Enskog expansion to recover the convective and source terms of the macroscopic equation. Numerical experiments with European call options with different parameters are performed to verify the accuracy and efficiency of the proposed models. The numerical results of the D1Q5 model are better than those of the D1Q3 model, and our results are better than the existing numerical results as well. The pricing of European call option, binary option, and butterfly spread option is efficiently simulated in the numerical experiments. Furthermore, our proposed lattice Boltzmann model can be used to solve a class of partial differential equations with variable coefficients and a source term.



---

*Source: 1812518-2023-01-18.xml*