Numerical Method for Multi-Dimensional Coupled Forward-Backward Stochastic Differential Equations Based on Fractional Fourier Fast Transform

: Forward-backward stochastic differential equations (FBSDEs) have received more and more attention in the past two decades. FBSDEs can be applied to many ﬁelds, such as economics and ﬁnance, engineering control, population dynamics analysis, and so on. In most cases, FBSDEs are nonlinear and high-dimensional and cannot be obtained as analytic solutions. Therefore, it is necessary and important to design their numerical approximation methods. In this paper, a novel numerical method of multi-dimensional coupled FBSDEs is proposed based on a fractional Fourier fast transform (FrFFT) algorithm, which is used to compute the Fourier and inverse Fourier transforms. For the forward component of FBSDEs, time discretization is used as well as the backward equation to yield a recursive system with terminal conditions. For the numerical experiments to be successful, three types of numerical methods were used to solve the problem, which ensured the efﬁciency and speed of computation. Finally, the numerical methods for different examples are veriﬁed.


Introduction
The study of effective numerical calculation methods for FBSDEs has become a hot topic in recent years due to their continuous development in both theory and practical applications. Because the solutions of FBSDEs are rarely found in real problems, numerical methods are a good choice. In 1990, Pardoux and Peng [1] proved the existence and uniqueness of the nonlinear BSDE solution under standard conditions. After that, the study of FBSDEs became more and more popular. In contrast, the research results of coupled FBSDEs are lower than those of decoupled FBSDEs.
Currently, because of the computational complexity in the field of numerical experiments, coupled FBSDEs develop more slowly than decoupled FBSDEs. Because the analytic solutions of FBSDEs are rarely found in real problems, numerical solutions have become a good choice. Many researchers have tried to solve this problem. Generally, it can be divided into the following three main types of strategies.
However, coupled FBSDEs are a concept that has a wider range (e.g., coupled FBSDEs were used to simulate large investors who influence the stock price of financial markets).
where X t = X 1 t , · · · , X l t T ∈ R l , (0 ≤ t ≤ T) is an ι-dimensional forward stochastic differential equation (FSDE), and Y t ∈ R, (0 ≤ t ≤ T) is a one-dimensional backward stochastic differential equation(BSDE). X 0 = x 0 and Y T = g(X T ) are the initial condition of the FSDE and the end condition of the BSDE, respectively. Z T t = (Z t,1 , · · · , Z t,k ) T ∈ R k is a control process. µ(t, X t , Y t , Z t ) = µ 1 t, X 1 t , Y t , Z t , · · · , µ n t, X l t , Y t , Z t are the drift terms of the FSDE, Σ(t, X t , Y t ) = Σ i,j (t, X t , Y t ) , 1 ≤ i ≤ l, 1 ≤ j ≤ k are diffusion terms of the FSDE, and W t = W 1 t , · · · , W k t T is defined on the probability space Ω, F, (F t ) 0≤t≤T , P with F t the natural σ-algebra generated by W t , which is a standard k-dimensional normal Brownian motion process and satisfies the usual conditions. Here, f (t, X t , Y t , Z t ) is called the generator function. The symbol (·) T denotes the transpose of a vector.
We first introduce some definitions and assumptions that underpin the proposed algorithm. Then, we elaborate the discretization of the normal coupled FBSDEs, which is the foundation for the following content.

Definitions and Assumptions
Following the survey papers [18,24,25], we provide the following standard definitions and assumptions. We define M 2 0, T; R l as the set of F T -adapted processes v(·) : Ω → R l that are square integrable, i.e., where |·| is the standard Euclidean norm in the Euclidean space. Denote L 1 T R l as the set of absolutely integrable processes. A triple set of processes (X,Y,Z) is an adapted solution of Equation (1) if it satisfies (1) almost surely and belongs to M 2 0, T; R l . In order to ensure the existence of the solution to Equation (1) and the validity of the numerical algorithm, we provide some assumptions [5] as follows: (1) The functions µ, Σ, f , and g are uniformly Lipschitz in (x, y, z), for t ∈ [0, T]. For example, µ should satisfy where t ∈ [0, T] and L Lipz > 0 is the Lipschitz constant. (2) The functions µ, f , and g satisfy the linear growth conditions as follows: where k 1 , k 2 , k 3 > 0 are constants.
(3) The diffusion Σ has bounded second derivatives with positive lower and upper bounds.
is a function of x and D x is the weak derivative symbol for x.
Assumptions (1)-(4) can ensure that the solution (x t , y t , z t ) to the coupled FBSDE (1) is existential and unique. If the function V(x) or D x V(x) does not satisfy assumption (4), we can truncate in bounds (for details, see [7]).
In the next subsection, we will discretize the FSDE in Equation (1) by using the Euler scheme and the BSDE in Equation (1) by θ-scheme.

Discretization of the General Coupled FBSDEs
It is classical to approximate the real solution of the FSDE by using the scheme of the Euler discretization. Consider a partition ∆ of time grid 0 = t 0 < t 1 < t 2 < · · · < t M = T with time steps ∆t m+1 : t m+1 − t m , for m = 1, 2, · · · , M − 1. For illustrative convenience, we denote X m = X t m , Y m = Y t m , Z m = Z t m and ∆W m = W t m+1 − W t m . We rewrite the local evolution of the FSDE in the time interval [t m , t m+1 ) to the following form: The discretization of the FSDE is: The Equation (3) can also be explained as a left point approximation of Equation (2).
where s T j,(·) is the j th row vector of s. For convenience and clarity, we consider the following substitutions: We rewrite Equations (7)-(9): We consider the Fourier transform of a function V(x) as: and the Fourier inverse transform in the same way: We then need to numerically calculate the above two equations, truncating the integral with a sufficiently large region ∏ ∈ R ι and ∏ ∈ R ι with: We get the approximation of integrals by the midpoint rule: where , j = 1, · · · , ι, n j ∈ 0, · · · , N − 1; , j = 1, · · · , ι, n j ∈ 0, · · · , N − 1.
We computeV(ξ) by Equation (14) on ξ-grid and V(x) by Equation (15) on x-grid. For efficiency, Equations (14) and (15) are computed by the FrFFT, which can be reviewed in [7,27]. We first introduce a one-dimensional situation and then generalize to a multidimensional one with the help of the one-dimensional situation. Let . The appendix in [7] has been demonstrated by the ι-multidimensional case of Equation (14). We will generalize Equation (15) to a ι-multidimensional situation. Consider the ι-dimensional integral, where To review this ι-dimensional sum, we can carry out the caculation as follows. where Each problem is able to be calculated by the one-dimensional FrFFT [28]. The Fourier transform is performed on both sides of Equations (11)- (13). In addition, we can get that (ŷ(t m , ξ),ẑ(t m , ξ)) is represented by To approximate the two equations, we introduce Lemma 1 and Theorem 1.
However, it still does not provide an implementable solution for the decoupling case. In the next section, we will provide the detailed numerical algorithms to implement the above theory.

Numerical Algorithms of Coupled FBSDE Combined with FrFFT
In this section, to verify the effectiveness of Theorem 1, we refer to three algorithms in [18]. In order to iterate, three algorithms need to be transformed.

Explicit Algorithm
We apply the technique from [30]. Suppose the time step ∆t m+1 = t m+1 − t m to be small enough. For the processes Y, Z, there is an approximation of the positive process: x, y(t m+1 , x), z(t m+1 , x))∆t m+1 + Σ(t m+1 , x, y(t m+1 , x))∆W m+1 . (23) Equation (23) differs from Equation (3), where the right-point approximation is applied. Because of the approximation of the different forward process X t , the characteristic functions of X ∆ m+1 − x are also different. However, the differences are minimal because we just need to replace t m in Equation (18) with t m+1 : Furthermore, we have to replace t m with t m+1 in Equations (10)- (13) and Equations (19)- (22). where Based on the above analysis, we give the following detailed Algorithm 1 to implement the numerical solution to Equation (1):

Initialization:
1: Calculateŷ(t m , ξ),ẑ(t m , ξ) andf (t m , ξ) applying Equation (14) by FrFFT on the ξ-grid. Main loop: 2: for m = M − 1 to 1 do 3: Calculateẑ(t m , ξ) andγ(t m+1 , ξ) via Equations (31) and (32). 4: Calculate z(t m , x) and γ(t m+1 , x) applying Equation (15) by FrFFT on the x-grid. 5: Do Picard iteration P times for Equation (26) to get y(t m , x) on the x-grid. Calculateŷ(t m , ξ) andγ(t m , ξ) applying Equation (14) by FrFFT on the ξ-grid. 9: end if 10: end for 11: if x 0 is not on the x-grid then 12: The cubic spline is interpolated on the x-grid to get the value at x 0 . 13: end if The following calculations are done when calculating the approximate values of the time t m : The overall complexity of the explicit algorithm is Different: Each iteration is made on the ξ-grid. Finally, at t 0 , the cubic spline is used to get the value at point x 0 .

Local Algorithm
We consider the iteration of each time interval [t m , t m+1 ) for m = 0, · · · , M − 1. Then one obtains an iteration that is similar to the method of [17].
Unlike in the literature [18], since the algorithm in this paper is performed on the grid for each iteration, each iteration needs to use the cubic spline to get the value of the corresponding point x at t m . This allows the next iteration to proceed smoothly.
We denote by P local the number of local iterations per time interval and k = 0, · · · , P local the current local iteration. Furthermore, we denote by , and γ(t m , x) in iteration k. Finally, we denote the values by y * (t m , x) and z * (t m , x) at X m on the x-grid using the cubic spline. The cubic spline is interpolated on the x-grid to get the value of y k * (t m , x) and z k * (t m , x) in local iteration. Then, Equation (18) will be rewritten by: However, we do not have functions with y k (t m , x) and z k (t m , x) at the iteration k and time point t m , so we apply the previous iteration for approximating Equation (33), i.e., (34) Furthermore, their derivations do not change. whereγ Based on the above analysis, we give the following detailed Algorithm 2 to implement the numerical solution to Equation (1): (14) by FrFFT on the ξ-grid. 2: Calculate y P local * (t M , x) and z P local * (t M , x) by terminal conditions. Main loop: (41) and (42). 8: Calculate z(t m , x) and γ(t m+1 , x) applying Equation (15) by FrFFT on the x-grid. 9: Do Picard iteration P times for Equation (36) to get y k (t m , x) on the x-grid. 10: Get f k (t m , x) on the x-grid.

11:
if x m is not on the x-grid then 12: The cubic spline is interpolated on the x-grid to get the value of y k * (t m , x) and z k Break.
Different: Each iteration uses the cubic spline to get the desired value (i.e.,y k * (t m , x)).

Global Algorithm
Finally, we can use one of the techniques proposed in [16] that uses an iterative process over the entire time domain. We denote by P global the number of global iterations, with k = 0, · · · , P global as the current global iteration. Furthermore, we denote by x) at X m on the x-grid using the cubic spline in the k th iteration. Then the Equation (18) will be rewritten by: However, we do not have functions with y k (t m , x) and z k (t m , x) at iteration k and time point t m , so we also apply the previous iteration: (44) Using same logic as in Section 4.2, we have: where whereγ Then the process of Algorithm 3 is as follows: Calculateẑ(t m , ξ) andγ(t m+1 , ξ) via Equations (51) and (52). 9: Calculate z(t m , x) and γ(t m+1 , x) applying Equation (15) by FrFFT on the x-grid. 10: Do Picard iteration P times for Equation (36) to get y k (t m , x) on the x-grid.

11:
Get f k (t m , x) on the x-grid.

12:
if x m is not on the x-grid then 13: The cubic spline is interpolated on the x-grid to get the value of y k * (t m , x) and z k * (t m , x) at x m . 14: end if 15: if m = 0 then 16: Calculateŷ k (t m , ξ) andf k (t m , ξ) applying Equation (14) by FrFFT on the ξ-grid. 17: end if 18: end for Break. 21: end if 21: end for The overall complexity of the global algorithm is O(MP global (DN + DN ι + DN ι log(N)+ PN ι + 1)).
Different: Like the local algorithm, each iteration uses the cubic spline to get the desired value (i.e., y k * (t m , x)).

Numerical Experiments
In this section, we analyze several examples using the numerical methods developed above (Algorithm 1, Algorithm 2 and Algorithm 3) and analyze their results. In principle, the method can be applied to any dimension. However, due to the curse of dimensionality, it can compute intuitively in low dimensions. So, in order to display the experimental results, we only consider one-dimensional components; that is, X t has one component. MATLAB 9.12.0.1884302 was used for the computations, with an Intel(R) Pentium(R) CPU 4415U @ 2.30GHz and 4.00 GB of RAM.
For the local algorithm and the global algorithm, we use stopping criteria to determine when iterations should be stopped. ε 1 and ε 2 are set to 10 −3 . Here we apply the same type of stopping criterion for Picard iterations. The errors in all experiments in this section are absolute errors, which are defined as Due to the Fourier method used in this paper, the integration interval [a, b] is referenced by [10,31] and the integration interval a, b is referenced by [7].
For each method, we tested the following θ schemes:

Example 1: A Decoupled Example
We consider the following decoupled FBSDE: where

The solution of Equation (53) is given by:
(y(t, x), z(t, x)) = sin 1 2 In this experiment, we use T = 1 and X 0 = 0 and set a = −10, b = 10, a = −45, and b = 45. It is known that the exact solution is (y(0, x 0 ), z(0, x 0 )) = 0, 1 2 . Since µ and Σ are constants, the result of this example can be obtained directly using ordinary numerical methods without three algorithms mentioned in Section 4. First, we analyze the error in the approximations of y(0, x 0 ) and z(0, x 0 ) of different values of N = 2 d . The error is listed in Table 1, where the approximations are obtained by the scheme A. Furthermore, the error in the approximation of y(0, x 0 ) and z(0, x 0 ) can be found in Figure 1.
In Table 1 and Figure 1, we find that the error is accepted. What is more, with the increase of M and N, the error gradually becomes smaller. The method is stable and the calculation time is fast.

Example 2: Forward Process Depending on X t
We consider the following coupled FBSDE: where  In this experiment, we use = 1 and 0 = 0 and set = −10, = 10,̃= −45, and ̃= 45. It is known that the exact solution is ( (0, 0 ), (0, 0 )) = (0, 1 2 ). Since and are constants, the result of this example can be obtained directly using ordinary numerical methods without three algorithms mentioned in Section 4. First, we analyze the error in the approximations of (0, 0 ) and (0, 0 ) of different values of = 2 . The error is listed in Table 1, where the approximations are obtained by the scheme A. Furthermore, the error in the approximation of (0, 0 ) and (0, 0 ) can be found in Figure 1.   The solution of Equation (54) is given by: For the experiment, we use T = 0.1 and X 0 = 0 and set a = −8, b = 8, a = −20, and b = 20. It is known that the exact solution is (y(0, x 0 ), z(0, x 0 )) = 1 2 , 1 8 . Because in the case of this example, the effect of the local algorithm and the global algorithm is the same, we only show the result of one of them as follows: Table 2 shows the error for the approximations of y(0, x 0 ) and z(0, x 0 ) for different values of M and N. In addition, we observe that the error in the approximation is mainly determined by the number of time points M; as the value of N increases, the error stops decreasing at a certain point. Therefore, we pay more attention to the error behavior when the number of time points changes. For all the next numerical tests, we will use N = 2 9 to ensure sufficient accuracy. In Figures 2 and 3, the errors of Y 0 achieve convergence as M becomes larger in all schemes. The errors of Z 0 also achieve convergence in schemes B and C. However, the error of the explicit algorithm is significantly larger than that of the local algorithm. In short, although the processes of schemes A and D of the error of Z 0 are oscillating, they eventually converge at a very small number.

Example 3: Diffusion Depending on Y t
We consider the following coupled FBSDE [5]:
The results of the explicit algorithm, the local algorithm, and the global algorithm can be found separately in Figures 4-6. The error of Y 0 does not differ significantly in the three algorithms and actually converges as the value of M becomes larger. For the error of Z 0 in all schemes, we observe that the explicit method converges. For the local algorithm and global algorithm, scheme B and scheme C are smoother than other schemes. However, they are still convergent. )).
The results of the explicit algorithm, the local algorithm, and the global algorithm can be found separately in Figures 4-6. The error of 0 does not differ significantly in the three algorithms and actually converges as the value of becomes larger. For the error of 0 in all schemes, we observe that the explicit method converges. For the local algorithm and global algorithm, scheme B and scheme C are smoother than other schemes. However, they are still convergent.       Table 3 shows the CPU times for the three algorithms. Since the time of each scheme (B, C, and D) in the same algorithm is approaching, these will not be shown here. Here, we note that the explicit algorithm is remarkably faster than others. The CPU speeds of the local algorithm and the global algorithm are close.
The solution of Equation (56) is given by:  Table 3 shows the CPU times for the three algorithms. Since the time of each scheme (B, C, and D) in the same algorithm is approaching, these will not be shown here. Here, we note that the explicit algorithm is remarkably faster than others. The CPU speeds of the local algorithm and the global algorithm are close.
The results for the three algorithms can be found in Figures 7-9. For y(0, x 0 ), we observe that the three algorithms obviously converge. For z(0, x 0 ), scheme B and scheme C in the three algorithms are smoother than other schemes. Scheme A and scheme D have local oscillations, but overall, their errors are still small.       Table 4 shows the CPU times for the explicit, local, and global algorithms. Here, only the case of scheme A is shown. We noticed that the explicit algorithm is remarkably faster than the others. The CPU speeds of the local algorithm and the global algorithm are similar.    We consider the following coupled FBSDE [17]: where The solution of Equation (57) is given by: (y(t, x), z(t, x)) = sin(t + x), In the experiment, we use T = 0.1 and X 0 = 0 and set a = −4.2307, b = 4.7135, a = −15, and b = 15. It is known that the exact solution is (y(0, x 0 ), z(0, x 0 )) = (0, 0).
The results for the explicit algorithm, the local algorithm, and the global algorithm can be found separately in Figures 10-12. For y(0, x 0 ), we observe that the effects of the three algorithms do not differ significantly. Although the error from the figures increases when M ≥ 200, they converge at a certain number. For z(0, x 0 ), scheme B and scheme C are obviously smoother than other schemes. Interestingly, all schemes are smooth, and the error converges to a number for the explicit algorithm. For the local and global algorithms, the error has local oscillations, and their results are not as accurate as those of scheme B and scheme C. Table 5 shows the CPU times for the explicit, local, and global algorithms. Since the time of each scheme (B, C, and D) in the same algorithm is approaching, only the case of scheme A is shown. Here, we also observe that the explicit algorithm is obviously faster than others. The CPU speeds of the local algorithm and the global algorithm are close.
can be found separately in Figures 10-12. For (0, 0 ), we observe that the effects of the three algorithms do not differ significantly. Although the error from the figures increases when ≥ 200, they converge at a certain number. For (0, 0 ), scheme B and scheme C are obviously smoother than other schemes. Interestingly, all schemes are smooth, and the error converges to a number for the explicit algorithm. For the local and global algorithms, the error has local oscillations, and their results are not as accurate as those of scheme B and scheme C.     Table 5 shows the CPU times for the explicit, local, and global algorithms. Since the time of each scheme (B, C, and D) in the same algorithm is approaching, only the case of scheme A is shown. Here, we also observe that the explicit algorithm is obviously faster    Table 5 shows the CPU times for the explicit, local, and global algorithms. Since the time of each scheme (B, C, and D) in the same algorithm is approaching, only the case of scheme A is shown. Here, we also observe that the explicit algorithm is obviously faster than others. The CPU speeds of the local algorithm and the global algorithm are close.  We consider the following coupled FBSDE: where This FBSDE satisfies assumptions (1)~(4). The solution of Equation (58) is given by: (y(t, x), z(t, x)) = sin(t + x), cos 2 (t + x) .
As the diffusion coefficient now depends on the process Z t , the nonlinear Feynman-Kac theorem is not applied. Therefore, we take scheme B and C in the first time step for the local algorithm and the global algorithm. The results for the three algorithms can be found separately in Figures 13-15. For y(0, x 0 ) in the explicit algorithm, we observe that all schemes converge and that the effect of the local algorithm is better than those of the explicit and global methods. For z(0, x 0 ), scheme B and scheme C are smoother than other schemes in the explicit algorithm. Scheme A and scheme D are undulant. However, the approximations of z(0, x 0 ) using the local and global methods have a divergent tendency. Their precision needs to be improved. For the experiment, we use = 1 and 0 = 0 and set = −2 , = 2 ,̃= −2 , and ̃= 2 . It is known that the exact solution is ( (0, 0 ), (0, 0 )) = (0,1).
As the diffusion coefficient now depends on the process , the nonlinear Feynman-Kac theorem is not applied. Therefore, we take scheme B and C in the first time step for the local algorithm and the global algorithm. The results for the three algorithms can be found separately in Figures 13-15. For (0, 0 ) in the explicit algorithm, we observe that all schemes converge and that the effect of the local algorithm is better than those of the explicit and global methods. For (0, 0 ), scheme B and scheme C are smoother than other schemes in the explicit algorithm. Scheme A and scheme D are undulant. However, the approximations of (0, 0 ) using the local and global methods have a divergent tendency. Their precision needs to be improved. Figure 13. Plot of the error of the explicit algorithm for example 6. Figure 13. Plot of the error of the explicit algorithm for example 6. Table 6 shows the CPU times for the explicit, local, and global algorithms. The speed of all algorithms is still relatively great.   Table 6 shows the CPU times for the explicit, local, and global algorithms. The speed of all algorithms is still relatively great.     Table 6 shows the CPU times for the explicit, local, and global algorithms. The speed of all algorithms is still relatively great.

Conclusions
In this paper, based on three different numerical algorithms, we extend the fractional Fourier fast transform to propose methods for the numerical solution of multi-dimensional coupled FBSDEs. In this method, we prove a theorem in which the complicated conditional expectation containing Brown motion can be switched to expectation with time points and simple information. In order to verify the validity, we give three algorithms: the explicit algorithm, the local algorithm, and the global algorithm. Some examples show that the explicit algorithm performs best over the local and global algorithms because it is the fastest and most accurate. The calculated cost of this paper is cheaper than that of the case without FrFFT. The cost of the algorithms for the local algorithm and the global algorithm is more expensive than that of the explicit algorithm. However, the results of the local algorithm and the global algorithm are better than those of the explicit algorithm, especially in the case of example 2 (i.e., forward process depending on X t ). In a word, in the case of diffusion depending on Z t , the accuracy of the numerical simulation of value of Z 0 must be improved.
In future research, the case of diffusion depending on Z t should be considered to improve accuracy, e.g., using Richardson extrapolation to accelerate the speed of convergence. In addition, FBSDEs with jumps can also be considered by applying the Fourier transform.