Abstract
In this study, we present an accurate and efficient nonuniform finite difference method for the threedimensional (3D) timefractional Black–Scholes (BS) equation. The operator splitting scheme is used to efficiently solve the 3D timefractional BS equation. We use a nonuniform grid for pricing 3D options. We compute the threeasset cashornothing European call option and investigate the effects of the fractionalorder in the timefractional BS model. Numerical experiments demonstrate the efficiency and fastness of the proposed scheme.
1. Introduction
We consider the following 3D version of the timefractional Black–Scholes (BS) model [1]: 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 [2]. 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 heavytailed 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. Timefractional analysis is closely connected to stable probability distributions [3]. Many researchers in the financial field have attempted to generalize the BS model in the fractionalorder 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 spacefractional [4, 5] and timefractional derivatives [6, 7]. Regarding the timefractional 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 twodimensional timefractional 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 [16]. Option derivatives, European vanilla options [17], and double barrier options [18] are analytically priced under the timefractional BS equation. In [19], 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 [20]. The authors in [21] proposed a novel operator splitting scheme for pricing American options using the timefractional 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 [22]. They demonstrated that the proposed scheme has unconditional stability and convergent property through the numerical results. In [23], 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 twodimensional fractional PDE. This method demonstrated the accuracy and efficiency of numerical studies. The author in [24] proposed the application of homotopy analysis method (HAM) for pricing European call option based on timefractional BS equation. He demonstrated the accuracy, effectiveness, and suitability of HAM through comparative tests. The pricing equation based on a spacetime fractional PDE is presented in [25]. The author calculated European call and put options based on spacetime fractional BS equation using the technique of Adomian decomposition method under the FDM. In [26], option derivatives were numerically priced using the method for the timefractional BS equation. These schemes are both firstorder and secondorder accurate in time and space, respectively. De Staelen and Hendy [27] improved the spatial fourthorder scheme with a temporal accuracy order of and performed stability and convergence analysis on their proposed scheme. Golbabai and Nikan [28] numerically solved the timefractional BS equation using the moving leastsquares method. The authors in [29] solved the fractional threedimensional (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. [30] modified an scheme to solve the timefractional BS equation. The modified time method is based on a change of variable and then obtains optimal error estimates. In [31], the authors removed the convection term with exponential transformation, transforming the timefractional BS model into a timefractional subdiffusion model, and then applied  formula for the Caputo timefractional 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 highdimensional problems. Although numerical studies have been conducted on the oneasset [26–28] and twoasset [16, 32] options, there is a lack of research on higherdimensional numerical methods of more than two assets. Therefore, in this paper, we present the 3D timefractional BS equation for pricing threeasset cashornothing 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 timefractional 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 timefractional 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 timefractional 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 threeasset cashornothing 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 threedimensional 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 [34]: 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 [34]. 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) i77700 CPU @3.60 GHz machine with 8 GB of memory.
3.1. Effect of FractionalOrder
In this subsection, we investigate the effects of the fractionalorder by considering oneasset cashornothing European call option. The payoff function of cashornothing 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 cashornothing European call option for different fractionalorders , 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 .
(a)
(b)
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)).
(a)
(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. ThreeAsset Options with Nonuniform Mesh
3.2.1. CashorNothing Option
We investigated threeasset cashornothing 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 threeasset cashornothing 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 piecewiseuniform 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 piecewiseuniform mesh structure defined as for pricing the threeasset cashornothing 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 threeasset cashornothing 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 piecewiseuniform 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. EquityLinked Security
We consider a threeasset equitylinked security (ELS) option that contains knockinbarrier (). 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 knockinbarrier. If the underlying asset did not touch the knockinbarrier, 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 [35] on the nonuniform mesh. For additional information on numerical testing, please refer to the thesis [35]. The nonuniform mesh was constructed using the piecewiseuniform 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.
4. Conclusions
In this study, we presented an efficient and accurate nonuniform FDM for the 3D timefractional BS equation. In numerical experiments, we investigated the effects of the fractionalorder by considering oneasset cashornothing 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 highdimensional options takes longer, there is a lack of research on higherdimensional numerical methods with more than three assets in the timefractional BS equation. We used the nonuniform implicit FDM with operator splitting scheme for pricing threeasset cashornothing 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 onedimensional timefractional BS equation [1, 27], there is a lack of research on multidimensional theoretical analysis of the timefractional BS equation. Therefore, for future work, we will perform the theoretical analysis of the multidimensional timefractional BS equation and compute various financial assets based on the multidimensional timefractional BS equation using the proposed method and continue to improve our method.
Appendix
The following Listing 1 is a MATLAB code for pricing the threeasset cashornothing European call option.

Data Availability
No data were used to support this study.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
The corresponding author (J.S. Kim) was supported by the Brain Korea 21 FOUR from the Ministry of Education of the Republic of Korea.