Nonuniform Finite Difference Scheme for the Three-Dimensional Time-Fractional Black – Scholes Equation

. In this study, we present an accurate and e ﬃ cient nonuniform ﬁ nite di ﬀ erence method for the three-dimensional (3D) time-fractional Black – Scholes (BS) equation. The operator splitting scheme is used to e ﬃ ciently 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 e ﬀ ects of the fractional-order α in the time-fractional BS model. Numerical experiments demonstrate the e ﬃ ciency and fastness of the proposed scheme.


Introduction
We consider the following 3D version of the time-fractional Black-Scholes (BS) model [1]: where uðx, y, z, tÞ is the option value at time t and u T ðx, y, zÞ is the payoff function at time t = T, where 0 < α < 1 and L BS = 1 2 σ 2 x x 2 ∂ 2 u ∂x 2 + 1 2 σ 2 y y 2 ∂ 2 u ∂y 2 + 1 2 σ 2 z z 2 ∂ 2 u ∂z 2 + ρ xy σ x σ y xy ∂ 2 u ∂x∂y + ρ yz σ y σ z yz ∂ 2 u ∂y∂z Here, x, y, and z, and σ x , σ y , and σ z are the prices and volatilities of the underlying assets x, y, and z, respectively. Additionally, ρ xy , ρ yz , and ρ zx are the correlation values between two subscript asset variables, and r 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 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 [3]. 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][9][10] and numerical [11][12][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 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 [16]. Option derivatives, European vanilla options [17], and double barrier options [18] are analytically priced under the time-fractional 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 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  Journal of Function Spaces 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 ( [27] improved the spatial fourth-order scheme with a temporal accuracy order of 2 − α and performed stability and convergence analysis on their proposed scheme. Golbabai and Nikan [28] numerically solved the time-fractional BS equation using the moving least-squares method. The authors in [29] 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. [30] modified an L1 scheme to solve the time-fractional BS equation. The modified L1 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 time-fractional BS model into a time-fractional subdiffusion model, and then applied L 1 -2 formula for the Caputo time-fractional derivative. This scheme applied a quadratic B-spline collocation scheme for space. By using the compact quadratic spline collocation (QSC) scheme, this scheme yields 3 − α-order and 4-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][27][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 τ = T − t; then, where η = T − ξ is used. Let Uðx, y, z, τÞ = uðx, y, z, T − τÞ; 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 Uðx, y, z, 0Þ = u T ðx, y, zÞ for ðx, y, z , τÞ ∈ Ω × ð0, T. 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-ornothing option.  Here, x 1 = y 1 = z 1 = 0, x N x = L x , y N y = L y , and z N z = L z . Δτ = T/N τ is the time step, and N τ is the number of time steps. Figure 1 illustrates an example of a three-dimensional nonuniform mesh. Let U n ijk be the numerical approximation of Uðx i , y j , z k , nΔτÞ and τ p = pΔτ. The left hand side term in Equation (7) can be approximated by the following numerical quadrature formula:

Numerical Solutions
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 H H  H H  8h  8h  4h  4h  2h 2h h

Journal of Function Spaces
Let δτ = ðΔτÞ α Γð2 − αÞ for simplicity of exposition; then we sequentially solve the following equations [34]: for 1 ≤ i ≤ N x , 1 ≤ j ≤ N y , and 1 ≤ k ≤ N z . 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):   Journal of Function Spaces

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.

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 K = 100 and the cash is C = 100. The parameter values are r = 0:03, σ = 0:3, and Δτ = 1/365. We use a uniform mesh h = 1 with L = 200. Let α U n i be the numerical approximation of the solution, where i = 1, ⋯, N x , n = 0, ⋯, N τ , and 0 < α ≤ 1. A linear boundary condition can be applied. Figure 3 illustrates the numerical solutions of cashor-nothing European call option for different fractionalorders α = 0:4,0:6,0:8, and 1:0. Figures 3(a) and 3(b) show the numerical results with a relatively short maturity T = 0:1 and long maturity T = 2:5, respectively. We can observe the different solution profiles according to α. Figures 4(a) and 4(b) show the differences in the numerical solutions between α = 1:0 and α = 0:4,0:6,0:8,1:0, i.e., 1:0 U Nτ − α U Nτ , for T = 0:1 and T = 2:5, 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 x = 100 with α = 0:4,0:6,0:8, and 1.0 as maturity increases, 0 ≤ T ≤ 2:5. 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.

Cash-or-Nothing Option.
We investigated three-asset cash-or-nothing European call option with the following payoff function:  Figure 6 shows the mesh structure with a mesh size h 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 x 2 < h, then we reset x 2 = 0:5x 3 so that we can apply the linear boundary condition. Similarly, if Table 1 presents the three-asset cash-or-nothing European call option prices with various variable h and time step Δτ. Here, we use α = 0:5. We can confirm that the option prices obtained with each h value converge as the time step Δτ becomes smaller. We adopt the reference solution Uðx i , y j , z k , TÞ which uses h = 1 and time step Δτ = 0:1/80.
We consider the piecewise-uniform mesh Ω Here, D ± i are the upper and lower bounds of each uniform mesh and are defined as follows: where 2 × m i is the number of points in mesh Ω i for i = 1, 2, 3, 4. In particular, m 4 = bD − 3 /ð8hÞc − 1 where bxc is the maximum integer not greater than x. From now on, we use m 1 = 5, m 2 = 5, and m 3 = 4 in our numerical experiments. Figure 7 shows the piecewise-uniform mesh structure defined as Ω x for pricing the three-asset cash-or-nothing option considered in this section. d ± i , which is defined in Figure 7, is an index of the point x with the D ± i values defined above. Figure 12: Nonuniform mesh for the three-asset ELS.

Journal of Function Spaces Journal of Function Spaces
Given the same option, Table 2 lists the option prices with respect to α. Here, h = 1 and Δτ = 0:1/80 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 T = 0:1, 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 h = 1. We use the maturity time T = 0:1, 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 h = 2, 3, ⋯, 8. 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 16h. Likewise, we define D ± 5 = D ± 4 ± ð16hÞm 5 , and m 5 = bD − 4 /ð16hÞc − 1. We compute the CPU time and price of using the nonuniform mesh with m 1 = 1, m 2 = 2, m 3 = 3, and m 4 = 1, 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 h = 8. 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.

Equity-Linked Security.
We consider a three-asset equity-linked security (ELS) option that contains knock-inbarrier (D). 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 (K 1 , K 2 , ⋯, K 6 ) on the early exercise date before maturity, the contract provides the specific returns (c 1 , c 2 , ⋯, c 6 ) 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 (d) and otherwise suffers losses. The payoff structure of ELS is illustrated in Figure 11.
For α = 0:8, 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 piecewise-uniform mesh, as shown in Figure 12, with fixed points (0, D/2, D, K 1 , K 2 , ⋯, K 6 , x0, L) and m 1 = 5, m 2 = 20, where x0 is the current underlying price.
When using a nonuniform grid as shown in Figure 12, the ELS price is 8767 and the elapsed time is 15.3783. When comparing this result to the result obtained using the uniform mesh (h = 2), the relative error of the price is 0.0205 and the elapsed time is 470 times shorter.

Conclusions
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 α = 1:0 and α < 1:0 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 higherdimensional 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.