Lattice Boltzmann Method for the Generalized Black-Scholes Equation

Wu, Fangfang; Xu, Duoduo; Tian, Ming; Wang, Yingying

doi:https://doi.org/10.1155/2023/1812518

Advances in Mathematical Physics

On this page

Abstract Introduction Preliminaries Conclusions Data Availability Conflicts of Interest Acknowledgments References Copyright Related Articles

Research Article | Open Access

Volume 2023 | Article ID 1812518 | https://doi.org/10.1155/2023/1812518

Lattice Boltzmann Method for the Generalized Black-Scholes Equation

Fangfang Wu,¹Duoduo Xu,¹Ming Tian,¹and Yingying Wang¹

Academic Editor: Zengtao Chen

Received29 Aug 2022

Revised24 Dec 2022

Accepted05 Jan 2023

Published18 Jan 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.

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] with the final condition where is a smooth function and boundary conditions where is the European option price at underlying asset price and at time , is the risk-free interest rate, is the dividend, and is the volatility. and are all the functions of and . is the exercise price, is the expiry date, and 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 Section 2, 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 and derivative rule, Equations (1)–(4) are equivalent to the following form:

with initial condition

and boundary conditions where

, , and denote the first-order partial derivatives of , , and with respect to , respectively. denotes the second-order partial derivative of with respect to . 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 where and are the local particle distribution function and the local equilibrium distribution function at time and position in direction , respectively, is used for amending the convection and source terms in the recovered macroscopic equation, is a superficial phase-space variable, is time step, and is the dimensionless relaxation time. The stability of the equation requires that [33].

To recover the macroscopic BSE correctly by using Equation (9), we assume that the equilibrium distribution function and particle distribution function satisfy local mass conservation law, i.e.,

Now, using the Taylor series expansion to Equation (9), we have

The multiscale Chapman-Enskog expansion up to first order in space and second order in time is applied; we have where is a small expansion parameter. Furthermore, and can be expressed as

Substituting Equations (12)–(15) into Equation (11), we obtain a series of partial differential equations in different orders of :

From Equations (10), (14), and (16), we have

In order to recover Equation (5), some constraints on the higher moments of the local equilibrium distribution function are imposed as

Meanwhile, the amending function satisfies

Here, and are the parameters to be determined in the following. According to Equations (15), (24), (25), and (26), we have

Summing Equation (17) over and using Equations (20), (21), and (24), we obtain

Similarly, summing Equation (18) over and using Equation (20), we obtain

According to Equations (20)–(22), (24), and (30), we have

According to Equations (24) and (25), we have

Substituting Equations (32) and (33) into Equation (31) yields

Summing Equation (19) over and using Equation (20), we obtain

Using Equations (20)–(26) and (30), we have

Substituting Equations (36) and (37) into Equation (35) yields

3.1. D1Q3 Model

For the D1Q3 model, , where with representing grid spacing. The macroscopic equation with the second order of truncation error is derived by taking Equation (30)Equation (34) and using Equation (27):

To recover Equation (5), we let

Then, Equation (39) becomes

Coupling Equations (20)–(22), we can derive the local equilibrium distribution function as follows:

And according to Equations (27) and (28), a simple case for the amending function is given as

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, . The macroscopic equation is derived by taking Equation (30)Equation (34)Equation (38) and using Equations (27) and (40): where the truncation error is

As 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 , a convergence higher than second order was constructed in this model.

Remark 1. When , , and , using Equation (30), we have Then, the truncation error .

Coupling Equations (20)–(23), we can derive the equilibrium distribution function as follows:

According to Equations (27)–(29), a simple case for the amending function is given as

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 functions are initialized by setting to be . 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 (), the root mean squar error (), and corresponding orders of convergence and are defined as follows: where and are the numerical solution and analytical solution on the grids in space and 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]: with where is chosen such a way that .

In the computation, we take , , for the D1Q3 model, and for the D1Q5 model. The computational domain is fixed on . 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.

(a)

(b)

To test the accuracy of the D1Q3 model, simulations are performed for the fixed time step and , respectively. The grid resolutions range from to ( to ). The dimensionless relaxation time is fixed for the computations. Figure 2 shows the log-log plots of the RMS and MAE versus grid spacing at . 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 and boundary conditions and as (a)In this case, the analytical solution of Equation (1) with , , , , and isIn the computation, we take , , and for the D1Q5 model. The computational domain is fixed on . 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 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 from to ( to ) 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 . 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 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 , , , , and , and the analytical solution for this case is the same as in Example 2(a)The computational domain is fixed on . 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 , , , , , and . In this case, the analytical solution is not knownIn the computation, we take , , , and . 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 and the boundary conditions (a)In this case, the analytical solution of Equation (1) with , , , , , and isIn the computation, we take , , and for the D1Q5 model. The computational domain is fixed on . 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 , , , , , and . In this case the analytical solution is not knownIn the computation, we take , , and . The computational domain is fixed on . 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 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 , , , , , , and . The final condition and boundary conditions for the butterfly spread option are as follows:In this case, the analytical solution is not known.
In the computation, we take , , and . The computational domain is fixed on . 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 , , , , , and . The final condition for the butterfly spread option is as follows: In this case, the analytical solution is not known also.
In the computation, we take , , , and . Figure 9 plots the time evolution of the value of the option, which is consistent with the results in Ref. [36].

(a)

(b)

(a)

(b)

(a)

(b)

(a)

(b)

(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 function 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.

Data Availability

The data used to support the findings of this study are included within the article.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This work was supported by the National Natural Science Foundation of China under Grant 62103289.

References

S. Chen and G. D. Doolen, “Lattice Boltzmann method for fluid flows,” Annual Review of Fluid Mechanics, vol. 30, no. 1, pp. 329–364, 1998.
View at: Publisher Site | Google Scholar
S. Succi, R. Benzi, A. Cali, and M. Vergassola, “The lattice Boltzmann equation: theory and application,” in Microscopic Simulations of Complex Hydrodynamic Phenomena, pp. 187–203, Springer, 1992.
View at: Publisher Site | Google Scholar
Y. H. Qian, S. Succi, and S. A. Orszag, “Recent advances in lattice Boltzmann computing,” Annual Reviews of Computational Physics, vol. III, pp. 195–242, 1995.
View at: Publisher Site | Google Scholar
S. Succi, The Lattice Boltzmann Equation: For Fluid Dynamics and Beyond, Oxford university press, 2001.
B. Shi and Z. Guo, “Lattice Boltzmann simulation of some nonlinear convection-diffusion equations,” Computers & Mathematics with Applications, vol. 61, no. 12, pp. 3443–3452, 2011.
View at: Publisher Site | Google Scholar
L. Wang, B. Shi, and Z. Chai, “Regularized lattice Boltzmann model for a class of convection-diffusion equations,” Physical Review E, vol. 92, no. 4, article 043311, 2015.
View at: Publisher Site | Google Scholar
Z. Chai, B. Shi, and Z. Guo, “A multiple-relaxation-time lattice Boltzmann model for general nonlinear anisotropic convection–diffusion equations,” Journal of Scientific Computing, vol. 69, no. 1, pp. 355–390, 2016.
View at: Publisher Site | Google Scholar
J. Zhang and G. Yan, “Lattice Boltzmann method for one and two-dimensional Burgers equation,” Physica A: Statistical Mechanics and its Applications, vol. 387, no. 19-20, pp. 4771–4786, 2008.
View at: Publisher Site | Google Scholar
F. Liu and S. Weiping, “Numerical solutions of two-dimensional Burgers' equations by lattice Boltzmann method,” Communications in Nonlinear Science and Numerical Simulation, vol. 16, no. 1, pp. 150–157, 2011.
View at: Publisher Site | Google Scholar
H. Lai and C. Ma, “A new lattice Boltzmann model for solving the coupled viscous Burgers' equation,” Physica A: Statistical Mechanics and its Applications, vol. 395, pp. 445–457, 2014.
View at: Publisher Site | Google Scholar
H. Lai and C. Ma, “Lattice Boltzmann model for generalized nonlinear wave equations,” Physical Review E, vol. 84, no. 4, article 046708, 2011.
View at: Publisher Site | Google Scholar
Y. Wang, D. Rui, and Z. Chai, “Lattice Boltzmann model for time-fractional nonlinear wave equations,” Advances in Applied Mathematics and Mechanics, vol. 14, no. 4, pp. 914–935, 2022.
View at: Publisher Site | Google Scholar
W. Fangfang, W. Shi, and F. Liu, “A lattice Boltzmann model for the Fokker-Planck equation,” Communications in Nonlinear Science and Numerical Simulation, vol. 17, no. 7, pp. 2776–2790, 2012.
View at: Publisher Site | Google Scholar
Q. Li, Z. Chai, and B. Shi, “Lattice Boltzmann model for a class of convection-diffusion equations with variable coefficients,” Computers & Mathematics with Applications, vol. 70, no. 4, pp. 548–561, 2015.
View at: Publisher Site | Google Scholar
Z.-Z. Lan, H. Wen-Qiang, and B.-L. Guo, “General propagation lattice Boltzmann model for a variable-coefficient compound KdV-Burgers equation,” Applied Mathematical Modelling, vol. 73, pp. 695–714, 2019.
View at: Publisher Site | Google Scholar
P. Wilmott, S. Howison, and J. Dewynne, The Mathematics of Financial Derivatives, Cambridge U, 2012.
J. M. Manale and F. M. Mahomed, “A simple formula for valuing American and European call and put options,” in Proceeding of the Hanno Rund Workshop on the Differential Equations, University of Natal.
View at: Google Scholar
R. C. Merton, “Theory of rational option pricing,” The Bell Journal of economics and management science, vol. 4, no. 1, pp. 141–183, 1973.
View at: Publisher Site | Google Scholar
F. Black and M. Scholes, “The pricing of options and corporate liabilities,” in World Scientific Reference on Contingent Claims Analysis in Corporate Finance: Volume 1: Foundations of CCA and Equity Valuation, pp. 3–21, World Scientific, 2019.
View at: Publisher Site | Google Scholar
P. Wilmott, J. Dewynne, and S. Howison, Option Pricing: Mathematical Models and Computation 1993, Oxford Financial Press: Oxford, 1993.
R. Kangro and R. Nicolaides, “Far field boundary conditions for Black–Scholes equations,” SIAM Journal on Numerical Analysis, vol. 38, no. 4, pp. 1357–1368, 2000.
View at: Publisher Site | Google Scholar
P. Lucas Jódar, J. C. C. Sevilla-Peris, and R. Sala, “A new direct method for solving the Black-Scholes equation,” Applied Mathematics Letters, vol. 18, no. 1, pp. 29–32, 2005.
View at: Publisher Site | Google Scholar
S. Ross, J. C. Cox, and M. Rubinstein, “Option pricing: a simplified approach,” A Simplified Approach, vol. 7, no. 3, pp. 229–263, 1979.
View at: Publisher Site | Google Scholar
J. C. Hull and A. White, “The use of the control variate technique in option pricing,” Journal of Financial and Quantitative Analysis, vol. 23, no. 3, pp. 237–251, 1988.
View at: Publisher Site | Google Scholar
Z. Cen and A. Le, “A robust and accurate finite difference method for a generalized Black-Scholes equation,” Journal of Computational and Applied Mathematics, vol. 235, no. 13, pp. 3728–3733, 2011.
View at: Publisher Site | Google Scholar
S. Chandra Sekhara Rao, “Numerical solution of generalized Black-Scholes model,” Applied Mathematics and Computation, vol. 321, pp. 401–421, 2018.
View at: Publisher Site | Google Scholar
P. Roul and V. P. Goura, “A new higher order compact finite difference method for generalised Black-Scholes partial differential equation: European call option,” Journal of Computational and Applied Mathematics, vol. 363, pp. 464–484, 2020.
View at: Publisher Site | Google Scholar
R. Mohammadi, “Quintic B-spline collocation approach for solving generalized Black-Scholes equation governing option pricing,” Computers & Mathematics with Applications, vol. 69, no. 8, pp. 777–797, 2015.
View at: Publisher Site | Google Scholar
C. D. Nyoumbi and A. Tambue, “Convergence of a fitted finite volume method for pricing two dimensional assets with stochastic volatilities,” in Computational and Applied Math, Research Gate, 2021.
View at: Google Scholar
A. Golbabai, O. Nikan, and T. Nikazad, “Numerical analysis of time fractional Black–Scholes European option pricing model arising in financial market,” Computational and Applied Mathematics, vol. 38, no. 4, pp. 1–24, 2019.
View at: Publisher Site | Google Scholar
O. Nikan, Z. Avazzadeh, and J. A. Tenreiro Machado, “Localized kernel-based meshless method for pricing financial options underlying fractal transmission system,” in Mathematical Methods in the Applied Sciences, Wiley, 2021.
View at: Publisher Site | Google Scholar
A. Golbabai and O. Nikan, “A computational method based on the moving least-squares approach for pricing double barrier options in a time-fractional Black–Scholes model,” Computational Economics, vol. 55, no. 1, pp. 119–141, 2020.
View at: Publisher Site | Google Scholar
J. D. Sterling and S. Chen, “Stability analysis of lattice Boltzmann methods,” Journal of Computational Physics, vol. 123, no. 1, pp. 196–206, 1996.
View at: Publisher Site | Google Scholar
Z.-L. Guo, C.-G. Zheng, and B.-C. Shi, “Non-equilibrium extrapolation method for velocity and pressure boundary conditions in the lattice Boltzmann method,” Chinese Physics, vol. 11, no. 4, pp. 366–374, 2002.
View at: Publisher Site | Google Scholar
M. K. Kadalbajoo, L. P. Tripathi, and A. Kumar, “A cubic B-spline collocation method for a numerical solution of the generalized Black-Scholes equation,” Mathematical and Computer Modelling, vol. 55, no. 3-4, pp. 1483–1505, 2012.
View at: Publisher Site | Google Scholar
S. Wang, “A novel fitted finite volume method for the Black-Scholes equation governing option pricing,” IMA Journal of Numerical Analysis, vol. 24, no. 4, pp. 699–720, 2004.
View at: Publisher Site | Google Scholar

Copyright

Copyright © 2023 Fangfang Wu et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

PDF Download Citation

Download other formats

Order printed copies

Views

368

Downloads

437

Citations