#### Abstract

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.

#### 1. Introduction

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  and numerical  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  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.

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.

#### 4. 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 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.

#### Appendix

The following Listing 1 is a MATLAB code for pricing the three-asset cash-or-nothing European call option.

 1 c l e a r ; c l c ; 2 L=200; x v o l=0. 3 ; y v o l=0. 3 ; z v o l=0. 3 ; r=0. 0 3 ; rho xy=0. 5 ; rho yz=0. 5 ; 3 rho zx=0. 5 ;K1=100;K2=100;K3=100;T=0.1 ; dt=0.1 /80;Nt=c e i l (T/dt ) ; 4 h=1;m1=5;m2=5;m3=4; xr=K1+0. 5 h : h :K1+0. 5 h+m1; 5 xr=[ xr (1:end -1) xr(end): 2 h : xr(end)+2hm2 ] ; 6 xr=[ xr (1:end -1) xr(end): 4 h : xr(end)+4hm3 ] ; 7 m4=f l o o r ( (L- xr(end)) /8) ; 8 xr=[ xr (1:end -1) xr(end): 8 h : xr(end)+8hm4 ] ; 9 i f xr(end)=K1, y>=K2, z>=K3, 1 ) =100;V=U; 16 alp=0. 5 ; s=1. 0 /( dtˆ alp gamma(2 - alp ) ) ; 17 ax=z e r o s (1 ,Nx-2 ) ; dx=ax ; cx=ax ; 18 f o r i =2:Nx-1 19 ax ( i -1 )=r x ( i ) hx( i ) /(hx( i -1 ) (hx ( i -1 )+hx ( i ) ) ) . . . 20 - ( x v o l x ( i ) ) ˆ2/(hx ( i -1 ) ( hx( i -1 )+hx ( i ) ) ) ; 21 dx ( i -1 )=s+( x v o l x ( i ) ) ˆ2/(hx( i -1 ) hx ( i ) ) . . . 22 - r x ( i ) ( hx( i ) -hx( i -1 ) ) /(hx ( i -1 ) hx( i ) )+r / 3 ; 23 cx ( i -1 )=- r x ( i ) hx ( i -1 ) /(hx ( i ) (hx ( i -1 )+hx( i ) ) ) . . . 24 - ( x v o l x ( i ) ) ˆ2/(hx ( i ) (hx ( i -1 )+hx ( i ) ) ) ; 25 end 26 dx ( 1 )=dx ( 1 )++2ax ( 1 ) ; cx ( 1 )=cx ( 1 ) - ax ( 1 ) ; 27 ax (Nx-2 )=ax (Nx-2 ) - cx (Nx-2 ) ; dx (Nx-2 )=dx(Nx-2 )++2cx (Nx-2 ) ; 28 bx=ax ; by=ax ; bz=ax ; 29 f o r n=1:Nt 30 F=z e r o s (Nx-2 ,Ny-2 ,Nz -2 ) ; 31 i f n>1 32 f o r j =1:n -1 33 F=F+((n - j+1) ˆ(1 - alp ) - (n - j ) ˆ(1 - alp ) ) (U( 2 :Nx-1 , 2 :Ny-1 , 2 :Nz -1 , j+1) . . . 34 -U( 2 :Nx-1 , 2 :Ny-1 , 2 :Nz -1 , j ) ) ; 35 end 36 end 37 V( : , : , : , n)=U( : , : , : , n) ; 38 f o r j =2:Ny-1 39 f o r k=2:Nz -1 40 f o r i =2:Nx-1 41 bx ( i -1 )=s V( i , j , k , n) - s F( i -1 , j -1 , k -1 ) . . . 42 +( rho xy x v o l y v o l x ( i ) y ( j ) (V( i +1, j +1,k , n)+V( i -1 , j -1 , k , n) . . . 43 -V( i -1 , j+1,k , n) -V( i +1, j -1 , k , n) ) / ( ( hx ( i ) hy( j ) )+(hx ( i -1 ) hy ( j ) ) . . . 44 +(hx( i ) hy( j -1 ) )+(hx ( i -1 ) hy( j -1 ) ) ) . . . 45 +rho yz y v o l z v o l y ( j ) z ( k ) (V( i , j+1,k+1,n)+V( i , j -1 , k -1 , n) . . . 46 -V( i , j -1 , k+1,n) -V( i , j +1,k -1 , n) ) / ( ( hy ( j ) hz ( k ) )+(hy ( j -1 ) hz ( k ) ) . . . 47 +(hy( j ) hz ( k -1 ) )+(hy ( j -1 ) hz ( k -1 ) ) ) . . . 48 +rho zx x v o l z v o l x ( i ) z ( k ) (V( i +1, j , k+1,n)+V( i -1 , j , k -1 , n) . . . 49 -V( i -1 , j , k+1,n) -V( i +1, j , k -1 , n) ) / ( ( hx ( i ) hz ( k ) )+(hx ( i -1 ) hz ( k ) ) . . . 50 +(hx( i ) hz ( k -1 ) )+(hx ( i -1 ) hz ( k -1 ) ) ) ) ; 51 end 52 U( 2 :Nx-1 , j , k , n+1)=thomas3 ( ax , dx , cx , bx ) ; 53 end 54 end 55 U( 1 , 2 :Ny-1 , 2 :Nz -1 , n+1)=2U( 2 , 2 :Ny-1 , 2 :Nz -1 , n+1) . . . 56 -U( 3 , 2 :Ny-1 , 2 :Nz -1 , n+1) ; 57 U( : , 1 , 2 :Nz -1 , n+1)=2U( : , 2 , 2 :Nz -1 , n+1) -U( : , 3 , 2 :Nz -1 , n+1) ; 58 U( : , : , 1 , n+1)=2U( : , : , 2 , n+1) -U( : , : , 3 , n+1) ; 59 U(Nx , 2 :Ny-1 , 2 :Nz -1 , n+1)=2U(Nx-1 , 2 :Ny-1 , 2 :Nz -1 , n+1) . . . 60 -U(Nx-2 , 2 :Ny-1 , 2 :Nz -1 , n+1) ; 61 U( : ,Ny , 2 :Nz -1 , n+1)=2U( : ,Ny-1 , 2 :Nz -1 , n+1) -U( : ,Ny-2 , 2 :Nz -1 , n+1) ; 62 U( : , : ,Nz , n+1)=2U( : , : ,Nz -1 , n+1) -U( : , : ,Nz -2 , n+1) ; 63 f o r k=2:Nz -1 64 f o r i =2:Nx-1 65 f o r j =2:Ny-1 66 by ( j -1 )=s U( i , j , k , n+1) ; 67 end 68 V( i , 2 :Ny-1 , k , n+1)=thomas3 ( ax , dx , cx , by ) ; 69 end 70 end 71 V( 1 , 2 :Ny-1 , 2 :Nz -1 , n+1)=2V( 2 , 2 :Ny-1 , 2 :Nz -1 , n+1) . . . 72 -V( 3 , 2 :Ny-1 , 2 :Nz -1 , n+1) ; 73 V( : , 1 , 2 :Nz -1 , n+1)=2V( : , 2 , 2 :Nz -1 , n+1) -V( : , 3 , 2 :Nz -1 , n+1) ; 74 V( : , : , 1 , n+1)=2V( : , : , 2 , n+1) -V( : , : , 3 , n+1) ; 75 V(Nx , 2 :Ny-1 , 2 :Nz -1 , n+1)=2V(Nx-1 , 2 :Ny-1 , 2 :Nz -1 , n+1) . . . 76 -V(Nx-2 , 2 :Ny-1 , 2 :Nz -1 , n+1) ; 77 V( : ,Ny , 2 :Nz -1 , n+1)=2V( : ,Ny-1 , 2 :Nz -1 , n+1) -V( : ,Ny-2 , 2 :Nz -1 , n+1) ; 78 V( : , : ,Nz , n+1)=2V( : , : ,Nz -1 , n+1) -V( : , : ,Nz -2 , n+1) ; 79 f o r j =2:Ny-1 80 f o r i =2:Nx-1 81 f o r k=2:Nz -1 82 bz (k -1 )=s V( i , j , k , n+1) ; 83 end 84 U( i , j , 2 :Nz -1 , n+1)=thomas3 ( ax , dx , cx , bz ) ; 85 end 86 end 87 U( 1 , 2 :Ny-1 , 2 :Nz -1 , n+1)=2U( 2 , 2 :Ny-1 , 2 :Nz -1 , n+1) . . . 88 -U( 3 , 2 :Ny-1 , 2 :Nz -1 , n+1) ; 89 U( : , 1 , 2 :Nz -1 , n+1)=2U( : , 2 , 2 :Nz -1 , n+1) -U( : , 3 , 2 :Nz -1 , n+1) ; 90 U( : , : , 1 , n+1)=2U( : , : , 2 , n+1) -U( : , : , 3 , n+1) ; 91 U(Nx , 2 :Ny-1 , 2 :Nz -1 , n+1)=2U(Nx-1 , 2 :Ny-1 , 2 :Nz -1 , n+1) . . . 92 -U(Nx-2 , 2 :Ny-1 , 2 :Nz -1 , n+1) ; 93 U( : ,Ny , 2 :Nz -1 , n+1)=2U( : ,Ny-1 , 2 :Nz -1 , n+1) -U( : ,Ny-2 , 2 :Nz -1 , n+1) ; 94 U( : , : ,Nz , n+1)=2U( : , : ,Nz -1 , n+1) -U( : , : ,Nz -2 , n+1) ; 95 end 96 f i g u r e ( 1 ) ; c l f ; colormap ( [ 0 0 0 ] ) ; 97 mesh ( x , y ,U( : , : , min ( f i n d ( z>100) )++1,n+1) ) ; 98 Pr i c e=i n t e rp 3 (y , x , z ,U( : , : , : , Nt+1) ,K1,K2,K3) 99 f u n c t i on x=thomas3 ( alpha , beta ,gamma, f ) 100 n=length ( f ) ; 101 f o r i =2:n 102 mult=alpha ( i ) / be ta ( i -1 ) ; 103 be ta ( i )=be ta ( i ) -multgamma( i -1 ) ; 104 f ( i )=f ( i ) -mult f ( i -1 ) ; 105 end 106 x (n)=f (n) / be ta (n) ; 107 f o r i=n -1 : -1 : 1 108 x ( i )=( f ( i ) -gamma( i ) x ( i+1) ) / be ta ( i ) ; 109 end 110 end

#### 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.