Numerical Methods for Differential and Integral EquationsView this Special Issue
Nonuniform Finite Difference Scheme for the Three-Dimensional Time-Fractional Black–Scholes Equation
In this study, we present an accurate and efficient nonuniform finite difference method for the three-dimensional (3D) time-fractional Black–Scholes (BS) equation. The operator splitting scheme is used to efficiently solve the 3D time-fractional BS equation. We use a nonuniform grid for pricing 3D options. We compute the three-asset cash-or-nothing European call option and investigate the effects of the fractional-order in the time-fractional BS model. Numerical experiments demonstrate the efficiency and fastness of the proposed scheme.
We consider the following 3D version of the time-fractional Black–Scholes (BS) model : where is the option value at time and is the payoff function at time , where and
Here, , , and , and , , and are the prices and volatilities of the underlying assets , , and , respectively. Additionally, , , and are the correlation values between two subscript asset variables, and is the interest rate. Black and Scholes published in 1973 their paper which described the BS model and option pricing formula . This has become an important fundamental topic for studying financial engineering and financial theory. However, the option pricing formula is based on the assumption that the returns of asset prices follow a Gaussian distribution. This means that the volatility of the underlying asset price is constant until the time to maturity of the option contracts. This has become a weakness of this formula. Many researchers and traders have found that rare events such as drastic drops in financial markets are much more frequent than would be anticipated based on Gaussian distributions and that the distribution of the returns of asset prices has a fat tail. Therefore, in real financial markets, many researchers have begun to develop models that more accurately reflect real market. The study of stable distribution has arisen naturally during the study of heavy-tailed distributions and has been applied in finance to develop models of extreme events that occur rarely. Because a stable probability distribution captures unpredictable events well, it is now more suitable for financial markets than the BS model. Time-fractional analysis is closely connected to stable probability distributions . Many researchers in the financial field have attempted to generalize the BS model in the fractional-order based on the fact that fractional derivatives and integrals provide powerful tools for explaining the memory and hereditary traits of different substances. The use of the fractional BS model for the high volatility of the stock market is one such generalization. There are two types of fractional derivatives as space-fractional [4, 5] and time-fractional derivatives [6, 7]. Regarding the time-fractional model, researchers have focused on the analytical [8–10] and numerical [11–13] methods. The finite difference method (FDM) is known as the most famous evaluation tool in quantitative finance and is more stable than Monte Carlo simulation (MCS). FDM has been applied in various studies [14, 15]. One researcher who has solved a two-dimensional time-fractional BS model using an implicit FDM proposed a fast biconjugate gradient stabilized scheme to solve the linear system to speed up computation and save storage space . Option derivatives, European vanilla options , and double barrier options  are analytically priced under the time-fractional BS equation. In , the authors developed a homotopy perturbation method to obtain the analytical solutions for the fractional BS equation. Khajehnasiri and Safavi presented the Boubaker operation matrix for the time fractional derivative which approximates the solution of the fractional BS . The authors in  proposed a novel operator splitting scheme for pricing American options using the time-fractional BS equation. They provided the effects of the fractional orders and the comparison of fractional equations through the numerical analysis. The paper also used the FDM of the Crank–Nicolson scheme for pricing European options based on the fractional BS equation . They demonstrated that the proposed scheme has unconditional stability and convergent property through the numerical results. In , the authors represented that the fractional partial differential equation (PDE) has been successfully applied in option pricing problems and it is more suitable for empirical financial markets. They used a fast preconditioned iterative method for pricing rainbow options based on a two-dimensional fractional PDE. This method demonstrated the accuracy and efficiency of numerical studies. The author in  proposed the application of homotopy analysis method (HAM) for pricing European call option based on time-fractional BS equation. He demonstrated the accuracy, effectiveness, and suitability of HAM through comparative tests. The pricing equation based on a space-time fractional PDE is presented in . The author calculated European call and put options based on space-time fractional BS equation using the technique of Adomian decomposition method under the FDM. In , option derivatives were numerically priced using the -method for the time-fractional BS equation. These schemes are both first-order and second-order accurate in time and space, respectively. De Staelen and Hendy  improved the spatial fourth-order scheme with a temporal accuracy order of and performed stability and convergence analysis on their proposed scheme. Golbabai and Nikan  numerically solved the time-fractional BS equation using the moving least-squares method. The authors in  solved the fractional three-dimensional (3D) chaotic process using the Adams–Bashforth–Moulton (ABM) method. They implemented an alternative numerical method based on the ABM method to reduce the computational cost and demonstrated that the proposed method is efficient and effective. She et al.  modified an scheme to solve the time-fractional BS equation. The modified time method is based on a change of variable and then obtains optimal error estimates. In , the authors removed the convection term with exponential transformation, transforming the time-fractional BS model into a time-fractional subdiffusion model, and then applied - formula for the Caputo time-fractional derivative. This scheme applied a quadratic -spline collocation scheme for space. By using the compact quadratic spline collocation (QSC) scheme, this scheme yields -order and -order convergence in time and space, respectively. The complexity of calculations and CPU time are very important when applying numerical methods to solve high-dimensional problems. Although numerical studies have been conducted on the one-asset [26–28] and two-asset [16, 32] options, there is a lack of research on higher-dimensional numerical methods of more than two assets. Therefore, in this paper, we present the 3D time-fractional BS equation for pricing three-asset cash-or-nothing European call option. Let us consider the following change of the variable ; then, where is used. Let ; then, Equation (5) becomes where we have used the integration by parts and the Leibniz integral rule. Therefore, after the change of variables, Equation (1) becomes with the initial condition for . When we solve the 3D time-fractional BS equation, there are difficulties in terms of memory shortage and computational cost because of the nonlocal property of the temporal derivative, which is the left hand side term in Equation (7). Therefore, we need efficient numerical schemes for this type of time-fractional PDE. First, the numerical scheme should be stable so that relatively large time steps can be used; otherwise, the computational cost will increase exponentially. Second, at each time step, the numerical solution scheme should be fast. To satisfy these conditions, in this study, we present an accurate and efficient nonuniform finite difference method for the 3D time-fractional BS model.
This paper is organized as follows. In Section 2, the proposed numerical scheme is described. In Section 3, numerical results are presented. In Section 4, conclusions are drawn. In the appendix, we provide the MATLAB code for the numerical implementation for the three-asset cash-or-nothing option.
2. Numerical Solutions
Let be the computational domain discretized in nonuniform intervals , and for , , and . Here, , and . is the time step, and is the number of time steps. Figure 1 illustrates an example of a three-dimensional nonuniform mesh.
Let be the numerical approximation of and . The left hand side term in Equation (7) can be approximated by the following numerical quadrature formula:
Therefore, we propose the following discretization of Equation (7) using Equation (8). where
The numerical derivatives are defined as and the other terms are similarly defined. Additional details can be found in [33, 34]. We solve the discrete Equation (9) using the operator splitting method. First, let where
Let for simplicity of exposition; then we sequentially solve the following equations : for , and . Note that if we sum up these three equations (14)–(16), then we obtain Equation (9). For the detailed numerical solution, algorithm with source program code of Equations (14)–(16) can be found in . We use the linear boundary condition, specifically, for example, in the case of Equation (14) (see Figure 2):
3. Numerical Experiments
Numerical experiments were conducted using MATLAB R2020b software on an Intel(R) Core(TM) i7-7700 CPU @3.60 GHz machine with 8 GB of memory.
3.1. Effect of Fractional-Order
In this subsection, we investigate the effects of the fractional-order by considering one-asset cash-or-nothing European call option. The payoff function of cash-or-nothing European call option is defined as where the strike price is and the cash is . The parameter values are , , and . We use a uniform mesh with . Let be the numerical approximation of the solution, where , , and . A linear boundary condition can be applied. Figure 3 illustrates the numerical solutions of cash-or-nothing European call option for different fractional-orders , and . Figures 3(a) and 3(b) show the numerical results with a relatively short maturity and long maturity , respectively. We can observe the different solution profiles according to .
Figures 4(a) and 4(b) show the differences in the numerical solutions between and , i.e., , for and , respectively. The lower the for a short maturity option, the higher the price of the option is in in the money (ITM) and undervalued in out of the money (OTM), see Figure 4(a). However, for long maturity option, the result is contrary to the result of short maturity option (see Figure 4(b)).
Figure 5 shows numerical solutions at with , and 1.0 as maturity increases, . For short maturity times, the solutions with lower values diffuse rapidly, but as the maturity time increases, the solutions with lower values diffuse slowly. The results of the two subaxes in Figure 5 can be interpreted similarly to Figure 4.
3.2. Three-Asset Options with Nonuniform Mesh
3.2.1. Cash-or-Nothing Option
We investigated three-asset cash-or-nothing European call option with the following payoff function: where the strike prices are and the cash is . The parameter values are , , , , , , and , , and . Figure 6 shows the mesh structure with a mesh size for the convergence test. Note that we straddle the strike point such that the strike point is in the middle of two neighboring points. If , then we reset so that we can apply the linear boundary condition. Similarly, if , then we reset .
Table 1 presents the three-asset cash-or-nothing European call option prices with various variable and time step . Here, we use . We can confirm that the option prices obtained with each value converge as the time step becomes smaller. We adopt the reference solution which uses and time step .
We consider the piecewise-uniform mesh with
Here, are the upper and lower bounds of each uniform mesh and are defined as follows: where is the number of points in mesh for . In particular, where is the maximum integer not greater than . From now on, we use , , and in our numerical experiments. Figure 7 shows the piecewise-uniform mesh structure defined as for pricing the three-asset cash-or-nothing option considered in this section. , which is defined in Figure 7, is an index of the point with the values defined above.
Given the same option, Table 2 lists the option prices with respect to . Here, and are taken and the other parameters are the same as in the test above. The code for the numerical implementation for this test is provided in the appendix.
Figure 8 shows the option prices according to the value of . For , the option prices obtained tend to be undervalued as decreases, as is the case with one underlying asset.
Figure 9 shows the CPU time and prices of three-asset cash-or-nothing option. In Figure 9, the dotted curve is the reference price of using the uniform mesh with mesh size . We use the maturity time , and the other parameters are the same as in the tests previously. In Figure 9, the solid and dashed curves are the CPU time and prices, respectively, with respect to uniform mesh with mesh size . Here, the uniform mesh is constructed as shown in Figure 6. Figure 10 was constructed in a similar manner to Figure 7. Here, we add a piecewise-uniform mesh with a step size . Likewise, we define , and . We compute the CPU time and price of using the nonuniform mesh with , and , which is constructed in Figure 10. We can confirm that the difference between the reference and numerical solutions obtained with each mesh is greater when using the uniform mesh, despite using the number of same grid points when the uniform mesh size is . Additionally, the elapsed time is similar to using the uniform mesh. In other words, nonuniform meshes are faster and more accurate compared to uniform mesh.
3.2.2. Equity-Linked Security
We consider a three-asset equity-linked security (ELS) option that contains knock-in-barrier (). The complex profit structure of ELS complicates pricing. To briefly explain return of ELS on one asset, if the underlying asset price is higher than the predetermined exercise prices () on the early exercise date before maturity, the contract provides the specific returns () and is terminated. Otherwise, the contract will continue on the next early exercise date. If the contract is not terminated by maturity, it depends on whether the underlying asset touched the knock-in-barrier. If the underlying asset did not touch the knock-in-barrier, it provides dummy return () and otherwise suffers losses. The payoff structure of ELS is illustrated in Figure 11.
For , we performed the comparison test with the uniform mesh under the same conditions considered in Section 3.2  on the nonuniform mesh. For additional information on numerical testing, please refer to the thesis . The nonuniform mesh was constructed using the piecewise-uniform mesh, as shown in Figure 12, with fixed points () and , where is the current underlying price.
When using a nonuniform grid as shown in Figure 12, the ELS price is and the elapsed time is 15.3783. When comparing this result to the result obtained using the uniform mesh (), the relative error of the price is 0.0205 and the elapsed time is 47 times shorter.
In this study, we presented an efficient and accurate nonuniform FDM for the 3D time-fractional BS equation. In numerical experiments, we investigated the effects of the fractional-order by considering one-asset cash-or-nothing European call option. The lower the value of , and the shorter the maturity of the option, and the larger the difference in option prices between the and except for at the money (ATM). Because of the complexity of calculations and CPU time for computation on high-dimensional options takes longer, there is a lack of research on higher-dimensional numerical methods with more than three assets in the time-fractional BS equation. We used the nonuniform implicit FDM with operator splitting scheme for pricing three-asset cash-or-nothing options and ELS. Here, we use the operator splitting method to solve the discrete system of equations and linear boundary conditions efficiently. Based on the use of the nonuniform implicit FDM, the numerical solution computation could be fast, and the numerical scheme could be stable even if relatively large time steps are used. Our results suggest that the proposed method is faster and more accurate than the uniform mesh. We demonstrated the efficiency and fastness of the proposed method through numerical experiments. Although there have been theoretical analyses (stability analysis, truncation error, and convergence analysis) of the one-dimensional time-fractional BS equation [1, 27], there is a lack of research on multidimensional theoretical analysis of the time-fractional BS equation. Therefore, for future work, we will perform the theoretical analysis of the multidimensional time-fractional BS equation and compute various financial assets based on the multidimensional time-fractional BS equation using the proposed method and continue to improve our method.
The following Listing 1 is a MATLAB code for pricing the three-asset cash-or-nothing European call option.
No data were used to support this study.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
The corresponding author (J.S. Kim) was supported by the Brain Korea 21 FOUR from the Ministry of Education of the Republic of Korea.
J. Huang, Z. Cen, and J. Zhao, “An adaptive moving mesh method for a time-fractional Black–Scholes equation,” Advances in Difference Equations, vol. 2019, no. 1, 14 pages, 2019.View at: Google Scholar
F. Black and M. Scholes, “The pricing of options and corporate liabilities,” Journal of Political Economy, vol. 81, no. 3, pp. 637–654, 1973.View at: Publisher Site | Google Scholar
A. A. Tateishi, H. V. Ribeiro, and E. K. Lenzi, “The role of fractional time-derivative operators on anomalous diffusion,” Frontiers in Physics, vol. 5, 2017.View at: Publisher Site | Google Scholar
Y. Chen, H. Gu, and L. Ma, “Variational method to -Laplacian fractional Dirichlet problem with instantaneous and noninstantaneous impulses,” Journal of Function Spaces, vol. 2020, Article ID 8598323, 8 pages, 2020.View at: Publisher Site | Google Scholar
D. Prathumwan and K. Trachoo, “On the solution of two-dimensional fractional Black–Scholes equation for European put option,” Advances in Difference Equations, vol. 2020, no. 1, 9 pages, 2020.View at: Google Scholar
Y. Kumar and V. K. Singh, “Computational approach based on wavelets for financial mathematical model governed by distributed order fractional differential equation,” Mathematics and Computers in Simulation, vol. 190, pp. 531–569, 2021.View at: Publisher Site | Google Scholar
M. X. Zhou, A. S. V. Kanth, K. Aruna et al., “Numerical solutions of time fractional Zakharov-Kuznetsov equation via natural transform decomposition method with nonsingular kernel derivatives,” Journal of Function Spaces, vol. 2021, Article ID 9884027, 17 pages, 2021.View at: Publisher Site | Google Scholar
S. Kumar, A. Yildirim, Y. Khan, H. Jafari, K. Sayevand, and L. Wei, “Analytical solution of fractional Black–Scholes European option pricing equation by using Laplace transform,” Journal of Fractional Calculus and Applications, vol. 2, no. 8, pp. 1–9, 2012.View at: Google Scholar
J. S. Duan, L. Lu, L. Chen, and Y. L. An, “Fractional model and solution for the Black–Scholes equation,” Mathematical Methods in the Applied Sciences, vol. 41, no. 2, pp. 697–704, 2018.View at: Google Scholar
K. Zhang, “Existence and uniqueness of analytical solution of time-fractional Black–Scholes type equation involving hyper-Bessel operator,” Mathematical Methods in the Applied Sciences, vol. 44, no. 7, pp. 6164–6177, 2021.View at: Publisher Site | Google Scholar
M. M. Khader, “The numerical solution for BVP of the liquid film flow over an unsteady stretching sheet with thermal radiation and magnetic field using the finite element method,” International Journal of Modern Physics C, vol. 30, no. 11, pp. 1950080–1950088, 2019.View at: Publisher Site | Google Scholar
Y. Gu and H. Sun, “A meshless method for solving three-dimensional time fractional diffusion equation with variable-order derivatives,” Applied Mathematical Modelling, vol. 78, pp. 539–549, 2020.View at: Publisher Site | Google Scholar
P. Roul, “A high accuracy numerical method and its convergence for time-fractional Black-Scholes equation governing European options,” Applied Numerical Mathematics, vol. 151, pp. 472–493, 2020.View at: Publisher Site | Google Scholar
M. M. Khader, “Fourth-order predictor-corrector FDM for the effect of viscous dissipation and Joule heating on the Newtonian fluid flow,” Computers & Fluids, vol. 182, pp. 9–14, 2019.View at: Publisher Site | Google Scholar
M. M. Khader and R. P. Sharma, “Evaluating the unsteady MHD micropolar fluid flow past stretching/shirking sheet with heat source and thermal radiation: implementing fourth order predictor-corrector FDM,” Mathematics and Computers in Simulation, vol. 181, pp. 333–350, 2021.View at: Publisher Site | Google Scholar
H. Zhang, F. Liu, S. Chen, and M. Shen, “A fast and high accuracy numerical simulation for a fractional Black–Scholes model on two assets,” Annals of Applied Mathematics, vol. 36, no. 1, pp. 91–110, 2020.View at: Google Scholar
W. Wyss, “The fractional Black–Scholes equation,” Fractional Calculus and Applied Analysis, vol. 3, no. 1, pp. 51–61, 2000.View at: Google Scholar
W. Chen, X. Xu, and S. P. Zhu, “Analytically pricing double barrier options based on a time-fractional Black- Scholes equation,” Computers & Mathematics with Applications, vol. 69, no. 12, pp. 1407–1419, 2015.View at: Publisher Site | Google Scholar
A. N. Fall, S. N. Ndiaye, and N. Sene, “Black-Scholes option pricing equations described by the Caputo generalized fractional derivative,” Chaos, Solitons & Fractals, vol. 125, pp. 108–118, 2019.View at: Publisher Site | Google Scholar
A. A. Khajehnasiri and M. Safavi, “Solving fractional Black-Scholes equation by using Boubaker functions,” Mathematical Methods in the Applied Sciences, vol. 4, no. 11, pp. 8505–8515, 2021.View at: Google Scholar
C. Chen, Z. Wang, and Y. Yang, “A new operator splitting method for American options under fractional Black- Scholes models,” Computers & Mathematics with Applications, vol. 77, no. 8, pp. 2130–2144, 2019.View at: Publisher Site | Google Scholar
G. Krzyżanowski, M. Magdziarz, and Ł. Płociniczak, “A weighted finite difference method for subdiffusive Black-Scholes model,” Computers & Mathematics with Applications, vol. 80, no. 5, pp. 653–670, 2020.View at: Publisher Site | Google Scholar
X. Chen, D. Ding, S. L. Lei, and W. Wang, “A fast preconditioned iterative method for two-dimensional options pricing under fractional differential models,” Computers & Mathematics with Applications, vol. 79, no. 2, pp. 440–456, 2020.View at: Publisher Site | Google Scholar
S. E. Fadugba, “Homotopy analysis method and its applications in the valuation of European call options with time-fractional Black-Scholes equation,” Chaos, Solitons & Fractals, vol. 141, article 110351, 2020.View at: Publisher Site | Google Scholar
L. Song, “A space-time fractional derivative model for European option pricing with transaction costs in fractal market,” Chaos, Solitons & Fractals, vol. 103, pp. 123–130, 2017.View at: Publisher Site | Google Scholar
X. Zhang, S. Shuzhen, W. Lifei, and Y. Xiaozhong, “-Difference numerical method for solving time-fractional Black–Scholes equation. Highlights of Science paper online,” China Science and Technology Papers, vol. 7, no. 13, pp. 1287–1295, 2014.View at: Google Scholar
R. H. De Staelen and A. S. Hendy, “Numerically pricing double barrier options in a time-fractional Black-Scholes model,” Computers & Mathematics with Applications, vol. 74, no. 6, pp. 1166–1175, 2017.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
M. S. Hashemi, M. Inc, and A. Yusuf, “On three-dimensional variable order time fractional chaotic system with nonsingular kernel,” Chaos, Solitons & Fractals, vol. 133, article 109628, 2020.View at: Publisher Site | Google Scholar
M. She, L. Li, R. Tang, and D. Li, “A novel numerical scheme for a time fractional Black-Scholes equation,” Journal of Applied Mathematics and Computing, vol. 66, no. 1-2, pp. 853–870, 2021.View at: Publisher Site | Google Scholar
Z. Tian, S. Zhai, H. Ji, and Z. Weng, “A compact quadratic spline collocation method for the time-fractional Black-Scholes model,” Journal of Applied Mathematics and Computing, vol. 66, no. 1-2, pp. 327–350, 2021.View at: Publisher Site | Google Scholar
W. Chen and S. Wang, “A 2nd-order ADI finite difference method for a 2D fractional Black-Scholes equation governing European two asset option pricing,” Mathematics and Computers in Simulation, vol. 171, pp. 279–293, 2020.View at: Publisher Site | Google Scholar
J. Bodeau, G. Riboulet, and T. Roncalli, “Non-uniform grids for PDE in finance,” 2000, https://papers.ssrn.com/sol3/papers.cfm?abstract_id=1031941.View at: Google Scholar
S. Kim, D. Jeong, C. Lee, and J. Kim, “Finite difference method for the multi-asset Black–Scholes equations,” Mathematics, vol. 8, no. 3, p. 391, 2020.View at: Publisher Site | Google Scholar
W. Lee, An efficient finite difference method for the three-dimensional time-fractional Black–Scholes equation, Graduate School, Korea University, 2021.