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

Authors: Sangkwon Kim; Chaeyoung Lee; Wonjin Lee; Soobin Kwak; Darae Jeong; Junseok Kim
Journal: Journal of Function Spaces
 (2021)

---

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


---

## Body

## 1. Introduction

We consider the following 3D version of the time-fractional Black–Scholes (BS) model [1]:
(1)∂αu∂tαx,y,z,t+LBSux,y,z,t=0forx,y,z,t∈Ω×0,T,(2)ux,y,z,T=uTx,y,z,where ux,y,z,t is the option value at time t and uTx,y,z is the payoff function at time t=T,
(3)∂αu∂tαx,y,z,t=1Γ1−αddt∫tTux,y,z,ξ−ux,y,z,Tξ−tαdξ,where 0<α<1 and
(4)LBS=12σx2x2∂2u∂x2+12σy2y2∂2u∂y2+12σz2z2∂2u∂z2+ρxyσxσyxy∂2u∂x∂y+ρyzσyσzyz∂2u∂y∂z+ρzxσzσxzx∂2u∂z∂x+rx∂u∂x+ry∂u∂y+rz∂u∂z−ru.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–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 [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 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 two-dimensional 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 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 [25]. 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 [26], 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 [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 L1 -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–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,
(5)∂αu∂tαx,y,z,t=1Γ1−αddt∫tTux,y,z,ξ−ux,y,z,Tξ−tαdξ=−1Γ1−αddτ∫T−τTux,y,z,ξ−ux,y,z,Tξ−T+ταdξ=−1Γ1−αddτ∫0τux,y,z,T−η−ux,y,z,Tτ−ηαdη,where η=T−ξ is used. Let Ux,y,z,τ=ux,y,z,T−τ; then, Equation (5) becomes
(6)−1Γ1−αddτ∫0τUx,y,z,η−Ux,y,z,0τ−ηαdη=−1Γ1−αddτ∫0τUx,y,z,ητ−ηαdη−ddτ∫0τUx,y,z,0τ−ηαdη=−1Γ1−αddτ∫0τUx,y,z,ητ−ηαdη−Ux,y,z,0τα=−1Γ1−αddτ∫0τ∂Ux,y,z,η∂ητ−η1−α1−αdη=−1Γ1−α∫0τ∂Ux,y,z,η∂ητ−η−αdη,where we have used the integration by parts and the Leibniz integral rule. Therefore, after the change of variables, Equation (1) becomes
(7)1Γ1−α∫0ττ−η−α∂Ux,y,z,η∂ηdη=LBSUx,y,z,τ,with the initial condition Ux,y,z,0=uTx,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 Section2, 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Ω=0,Lx×0,Ly×0,Lz be the computational domain discretized in nonuniform intervals hix=xi+1−xi,hjy=yj+1−yj, and hkz=zk+1−zk for i=1,⋯,Nx−1, j=1,⋯,Ny−1, and k=1,⋯,Nz−1. Here, x1=y1=z1=0,xNx=Lx,yNy=Ly, and zNz=Lz. Δτ=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.Figure 1
Three-dimensional nonuniform mesh.LetUijkn be the numerical approximation of Uxi,yj,zk,nΔτ and τp=pΔτ. The left hand side term in Equation (7) can be approximated by the following numerical quadrature formula:
(8)1Γ1−α∫0τn+1τn+1−η−α∂Uxi,yj,zk,η∂ηdη=1Γ1−α∑p=1n+1∫τp−1τpτn+1−η−α∂Uxi,yj,zk,η∂ηdη≈1Γ2−α∑p=1n+1τn+1−τp−11−α−τn+1−τp1−αUijkp−Uijkp−1Δτ=1ΔταΓ2−α∑p=1n+1n+2−p1−α−n+1−p1−αUijkp−Uijkp−1.Therefore, we propose the following discretization of Equation (7) using Equation (8).
(9)1ΔταΓ2−α∑p=1n+1n+2−p1−α−n+1−p1−αUijkp−Uijkp−1=LBSxUijkn+1/3+LBSyUijkn+2/3+LBSzUijkn+1,where
(10)LBSxUijkn+1/3=σxxi22DxxUijkn+1/3+rxiDxUijkn+1/3+σxσyρxyxiyjDxyUijkn+σyσzρyzyjzkDyzUijkn+σzσxρzxzkxiDzxUijkn−13rUijkn+1/3,LBSyUijkn+2/3=σyyj22DyyUijkn+2/3+ryjDyUijkn+2/3−13rUijkn+2/3,LBSzUijkn+1=σzzk22DzzUijkn+1+rzkDzUijkn+1−13rUijkn+1.The numerical derivatives are defined as(11)DxUijk=−hixUi−1,jkhi−1xhi−1x+hix+hix−hi−1xUijkhi−1xhix+hi−1xUi+1,jkhixhi−1x+hix,DxxUijk=2Ui−1,jkhi−1xhi−1x+hix−2Uijkhi−1xhix+2Ui+1,jkhixhi−1x+hix,DxyUijk=Ui+1,j+1,k−Ui−1,j+1,k−Ui+1,j−1,k+Ui−1,j−1,khixhjy+hi−1xhjy+hixhj−1y+hi−1xhj−1y,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
(12)1ΔταΓ2−α∑p=1n+1n+2−p1−α−n+1−p1−αUijkp−Uijkp−1=Fijkn+Uijkn+1−UijknΔταΓ2−α,where
(13)Fijkn=1ΔταΓ2−α∑p=1nn+2−p1−α−n+1−p1−αUijkp−Uijkp−1.Letδτ=ΔταΓ2−α for simplicity of exposition; then we sequentially solve the following equations [34]:
(14)Uijkn+1/3−Uijknδτ=LBSxUijkn+1/3−Fijkn,(15)Uijkn+2/3−Uijkn+1/3δτ=LBSyUijkn+2/3,(16)Uijkn+1−Uijkn+2/3δτ=LBSzUijkn+1,for 1≤i≤Nx,1≤j≤Ny, and 1≤k≤Nz. 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):
(17)U1jkn+1/3=2U2jkn+1/3−U3jkn+1/3,UNxjkn+1/3=2UNx−1,jkn+1/3−UNx−2,jkn+1/3,forj=2,⋯,Ny−1,k=2,⋯,Nz−1,Ui1kn+1/3=2Ui2kn+1/3−Ui3kn+1/3,UiNykn+1/3=2Ui,Ny−1,kn+1/3−Ui,Ny−2,kn+1/3,fori=1,⋯,Nx,k=2,⋯,Nz−1,Uij1n+1/3=2Uij2n+1/3−Uij3n+1/3,UijNzn+1/3=2Uij,Nz−1n+1/3−Uij,Nz−2n+1/3,fori=1,⋯,Nx,j=1,⋯,Ny.Figure 2
Schematic illustration of the linear boundary condition.


## 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
(18)U0x=c,ifC≥K,0,otherwise,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 αUin be the numerical approximation of the solution, where i=1,⋯,Nx, n=0,⋯,Nτ, and 0<α≤1. A linear boundary condition can be applied. Figure 3 illustrates the numerical solutions of cash-or-nothing European call option for different fractional-orders α=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 α.Figure 3
Numerical results for cash-or-nothing option for different fractional-ordersα=0.4,0.6,0.8, and 1.0 with maturity times (a) T=0.1 and (b) T=2.5.
(a)(b)Figures4(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.0UNτ−αUNτ, 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 4
Differences in the numerical solutions betweenα=1.0 and α=0.4,0.6,0.8, and 1.0 for maturity times (a) T=0.1 and (b) T=2.5.
(a)(b)Figure5 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.Figure 5
Numerical solution atx=100 for maturity 0≤T≤2.5 with α=0.4,0.6,0.8, and 1.0.


### 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:(19)U0x,y,z=c,ifx≥K1,y≥K2,z≥K3,0,otherwise,where the strike prices are K1=K2=K3=100 and the cash is c=100. The parameter values are T=0.1, r=0.03, L=M=N=200, σx=0.3, σy=0.3, σz=0.3, and ρxy=0.5, ρyz=0.5, and ρzx=0.5. 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 x2<h, then we reset x2=0.5x3 so that we can apply the linear boundary condition. Similarly, if xNx−xNx−1<h, then we reset xNx−1=0.5xNx+xNx−2.Figure 6
Mesh with a mesh sizeh for the convergence test.Table1 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 Uxi,yj,zk,T which uses h=1 and time step Δτ=0.1/80.Table 1
Three-asset cash-or-nothing option prices with varyingh and time step Δτ.
h=8h=4h=2h=1Δτ=0.1/1021.530122.604622.907622.9860Δτ=0.1/2021.457222.543322.854822.9358Δτ=0.1/4021.468922.589422.917623.0039Δτ=0.1/8021.507022.667523.015723.1082Δτ=0.1/16021.549022.747923.116723.2161We consider the piecewise-uniform meshΩx=Ω1∪Ω2∪Ω3∪Ω4∪Ω5 with
(20)Ω1=D1−,D1−+h,D1−+2h,⋯,D1+,Ω2=D2−,D2−+2h,⋯,D1−∪D1+,D1++2h,⋯,D2+,Ω3=D3−,D3−+4h,⋯,D2−∪D2+,D2++4h,⋯,D3+,Ω4=D4−,D4−+8h,⋯,D3−∪D3+,D3++8h,⋯,D4+,Ω5=0,0.5D4−,D4−∪D4+,D4++0.5L−D4+,L.Here,Di± are the upper and lower bounds of each uniform mesh and are defined as follows:
(21)D1−=K1−0.5+m1h,D1+=K1+0.5+m1h,D2−=D1−−2hm2,D2+=D1++2hm2,D3−=D2−−4hm3,D3+=D2++4hm3,D4−=D3−−8hm4,D4+=D3++8hm4,where 2×mi is the number of points in mesh Ωi for i=1,2,3,4. In particular, m4=D3−/8h−1 where x is the maximum integer not greater than x. From now on, we use m1=5, m2=5, and m3=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. di±, which is defined in Figure 7, is an index of the point x with the Di± values defined above.Figure 7
Piecewise-uniform meshΩx.Given the same option, Table2 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.Table 2
Option prices with respect toα.
α0.10.20.30.40.5Price20.653521.381822.040222.607623.0786α0.60.70.80.91.0Price23.459723.764524.008624.205324.3633Figure8 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 8
Plot of option prices with respect toα.Figure9 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 piecewise-uniform mesh with a step size 16h. Likewise, we define D5±=D4±±16hm5, and m5=D4−/16h−1. We compute the CPU time and price of using the nonuniform mesh with m1=1,m2=2,m3=3, and m4=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.Figure 9
CPU time and price of three-asset cash-or-nothing option for nonuniform mesh withm1=1,m2=2, and m3=3,m4=1 and uniform mesh with respect to h=2,3,⋯,8.Figure 10
Nonuniform mesh in a comparison test.


#### 3.2.2. Equity-Linked Security

We consider a three-asset equity-linked security (ELS) option that contains knock-in-barrier (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 (K1,K2,⋯,K6) on the early exercise date before maturity, the contract provides the specific returns (c1,c2,⋯,c6) 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.Figure 11
ELS payoff structure.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,K1,K2,⋯,K6,x0,L) and m1=5,m2=20, where x0 is the current underlying price.Figure 12
Nonuniform mesh for the three-asset ELS.When using a nonuniform grid as shown in Figure12, 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.


## 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
(18)U0x=c,ifC≥K,0,otherwise,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 αUin be the numerical approximation of the solution, where i=1,⋯,Nx, n=0,⋯,Nτ, and 0<α≤1. A linear boundary condition can be applied. Figure 3 illustrates the numerical solutions of cash-or-nothing European call option for different fractional-orders α=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 α.Figure 3
Numerical results for cash-or-nothing option for different fractional-ordersα=0.4,0.6,0.8, and 1.0 with maturity times (a) T=0.1 and (b) T=2.5.
(a)(b)Figures4(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.0UNτ−αUNτ, 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 4
Differences in the numerical solutions betweenα=1.0 and α=0.4,0.6,0.8, and 1.0 for maturity times (a) T=0.1 and (b) T=2.5.
(a)(b)Figure5 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.Figure 5
Numerical solution atx=100 for maturity 0≤T≤2.5 with α=0.4,0.6,0.8, and 1.0.


## 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:(19)U0x,y,z=c,ifx≥K1,y≥K2,z≥K3,0,otherwise,where the strike prices are K1=K2=K3=100 and the cash is c=100. The parameter values are T=0.1, r=0.03, L=M=N=200, σx=0.3, σy=0.3, σz=0.3, and ρxy=0.5, ρyz=0.5, and ρzx=0.5. 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 x2<h, then we reset x2=0.5x3 so that we can apply the linear boundary condition. Similarly, if xNx−xNx−1<h, then we reset xNx−1=0.5xNx+xNx−2.Figure 6
Mesh with a mesh sizeh for the convergence test.Table1 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 Uxi,yj,zk,T which uses h=1 and time step Δτ=0.1/80.Table 1
Three-asset cash-or-nothing option prices with varyingh and time step Δτ.
h=8h=4h=2h=1Δτ=0.1/1021.530122.604622.907622.9860Δτ=0.1/2021.457222.543322.854822.9358Δτ=0.1/4021.468922.589422.917623.0039Δτ=0.1/8021.507022.667523.015723.1082Δτ=0.1/16021.549022.747923.116723.2161We consider the piecewise-uniform meshΩx=Ω1∪Ω2∪Ω3∪Ω4∪Ω5 with
(20)Ω1=D1−,D1−+h,D1−+2h,⋯,D1+,Ω2=D2−,D2−+2h,⋯,D1−∪D1+,D1++2h,⋯,D2+,Ω3=D3−,D3−+4h,⋯,D2−∪D2+,D2++4h,⋯,D3+,Ω4=D4−,D4−+8h,⋯,D3−∪D3+,D3++8h,⋯,D4+,Ω5=0,0.5D4−,D4−∪D4+,D4++0.5L−D4+,L.Here,Di± are the upper and lower bounds of each uniform mesh and are defined as follows:
(21)D1−=K1−0.5+m1h,D1+=K1+0.5+m1h,D2−=D1−−2hm2,D2+=D1++2hm2,D3−=D2−−4hm3,D3+=D2++4hm3,D4−=D3−−8hm4,D4+=D3++8hm4,where 2×mi is the number of points in mesh Ωi for i=1,2,3,4. In particular, m4=D3−/8h−1 where x is the maximum integer not greater than x. From now on, we use m1=5, m2=5, and m3=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. di±, which is defined in Figure 7, is an index of the point x with the Di± values defined above.Figure 7
Piecewise-uniform meshΩx.Given the same option, Table2 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.Table 2
Option prices with respect toα.
α0.10.20.30.40.5Price20.653521.381822.040222.607623.0786α0.60.70.80.91.0Price23.459723.764524.008624.205324.3633Figure8 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 8
Plot of option prices with respect toα.Figure9 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 piecewise-uniform mesh with a step size 16h. Likewise, we define D5±=D4±±16hm5, and m5=D4−/16h−1. We compute the CPU time and price of using the nonuniform mesh with m1=1,m2=2,m3=3, and m4=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.Figure 9
CPU time and price of three-asset cash-or-nothing option for nonuniform mesh withm1=1,m2=2, and m3=3,m4=1 and uniform mesh with respect to h=2,3,⋯,8.Figure 10
Nonuniform mesh in a comparison test.


### 3.2.2. Equity-Linked Security

We consider a three-asset equity-linked security (ELS) option that contains knock-in-barrier (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 (K1,K2,⋯,K6) on the early exercise date before maturity, the contract provides the specific returns (c1,c2,⋯,c6) 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.Figure 11
ELS payoff structure.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,K1,K2,⋯,K6,x0,L) and m1=5,m2=20, where x0 is the current underlying price.Figure 12
Nonuniform mesh for the three-asset ELS.When using a nonuniform grid as shown in Figure12, 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.


## 3.2.1. Cash-or-Nothing Option

We investigated three-asset cash-or-nothing European call option with the following payoff function:(19)U0x,y,z=c,ifx≥K1,y≥K2,z≥K3,0,otherwise,where the strike prices are K1=K2=K3=100 and the cash is c=100. The parameter values are T=0.1, r=0.03, L=M=N=200, σx=0.3, σy=0.3, σz=0.3, and ρxy=0.5, ρyz=0.5, and ρzx=0.5. 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 x2<h, then we reset x2=0.5x3 so that we can apply the linear boundary condition. Similarly, if xNx−xNx−1<h, then we reset xNx−1=0.5xNx+xNx−2.Figure 6
Mesh with a mesh sizeh for the convergence test.Table1 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 Uxi,yj,zk,T which uses h=1 and time step Δτ=0.1/80.Table 1
Three-asset cash-or-nothing option prices with varyingh and time step Δτ.
h=8h=4h=2h=1Δτ=0.1/1021.530122.604622.907622.9860Δτ=0.1/2021.457222.543322.854822.9358Δτ=0.1/4021.468922.589422.917623.0039Δτ=0.1/8021.507022.667523.015723.1082Δτ=0.1/16021.549022.747923.116723.2161We consider the piecewise-uniform meshΩx=Ω1∪Ω2∪Ω3∪Ω4∪Ω5 with
(20)Ω1=D1−,D1−+h,D1−+2h,⋯,D1+,Ω2=D2−,D2−+2h,⋯,D1−∪D1+,D1++2h,⋯,D2+,Ω3=D3−,D3−+4h,⋯,D2−∪D2+,D2++4h,⋯,D3+,Ω4=D4−,D4−+8h,⋯,D3−∪D3+,D3++8h,⋯,D4+,Ω5=0,0.5D4−,D4−∪D4+,D4++0.5L−D4+,L.Here,Di± are the upper and lower bounds of each uniform mesh and are defined as follows:
(21)D1−=K1−0.5+m1h,D1+=K1+0.5+m1h,D2−=D1−−2hm2,D2+=D1++2hm2,D3−=D2−−4hm3,D3+=D2++4hm3,D4−=D3−−8hm4,D4+=D3++8hm4,where 2×mi is the number of points in mesh Ωi for i=1,2,3,4. In particular, m4=D3−/8h−1 where x is the maximum integer not greater than x. From now on, we use m1=5, m2=5, and m3=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. di±, which is defined in Figure 7, is an index of the point x with the Di± values defined above.Figure 7
Piecewise-uniform meshΩx.Given the same option, Table2 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.Table 2
Option prices with respect toα.
α0.10.20.30.40.5Price20.653521.381822.040222.607623.0786α0.60.70.80.91.0Price23.459723.764524.008624.205324.3633Figure8 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 8
Plot of option prices with respect toα.Figure9 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 piecewise-uniform mesh with a step size 16h. Likewise, we define D5±=D4±±16hm5, and m5=D4−/16h−1. We compute the CPU time and price of using the nonuniform mesh with m1=1,m2=2,m3=3, and m4=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.Figure 9
CPU time and price of three-asset cash-or-nothing option for nonuniform mesh withm1=1,m2=2, and m3=3,m4=1 and uniform mesh with respect to h=2,3,⋯,8.Figure 10
Nonuniform mesh in a comparison test.


## 3.2.2. Equity-Linked Security

We consider a three-asset equity-linked security (ELS) option that contains knock-in-barrier (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 (K1,K2,⋯,K6) on the early exercise date before maturity, the contract provides the specific returns (c1,c2,⋯,c6) 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.Figure 11
ELS payoff structure.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,K1,K2,⋯,K6,x0,L) and m1=5,m2=20, where x0 is the current underlying price.Figure 12
Nonuniform mesh for the three-asset ELS.When using a nonuniform grid as shown in Figure12, 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.


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



---

*Source: 9984473-2021-12-24.xml*