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

Abstract

Forward-backward stochastic differential equations (FBSDEs) have received more and more attention in the past two decades. FBSDEs can be applied to many fields, such as economics and finance, 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 efficiency and speed of computation. Finally, the numerical methods for different examples are verified.

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

The first strategy is to use fixed point theory (one can refer to the literature [2,3]). The second strategy utilizes the parabolic PDE to solve the FBSDE based on the general stochastic parabolic PDE discretization (for details, see [4,5,6]). The third strategy is based on Fourier transform (such as [7,8,9,10]). Among these numerical methods, some are higher-order numerical methods that have high convergence rates, such as [11,12,13,14,15].

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). There are also some methods to solve this problem. Bender and Zhang [16], who proved the convergence of time discretization and Markov iteration, laid the foundation for the numerical algorithm simulation of multidimensional coupled FBSDEs. In addition, based on FBSDE theory, Zhao, Fu, and Zhou [17] investigated a new type of high-order multi-step scheme for coupled FBSDEs. Recently, Huijskens, Ruijter, and Oosterlee [18] expanded the Fourier cos method to coupled FBSDEs with three numerical methods (i.e., explicit method, local method and global method), which had expected convergence and efficiency. Nam and Xu [19] proved the solution of coupled FBSDEs is existential by using purely probabilistic methods to obtain semilinear parabolic PDEs. Inspired by these papers, we expand the idea that multidimensional coupled FBSDEs can be solved using the fractional Fourier transform to obtain the numerical solution of Equation (1). In addition, in recent years, deep learning methods have also been introduced in the numerical field [20,21,22]. These algorithms are highly accurate, but at the same time, they are expensive to run. Therefore, this paper differs from other papers in the following aspects:

(1)

Proving the theorem in different ways to ensure that the algorithm of fractional fast Fourier transform can be retained during numerical experiments.

(2)

The numerical experiments were verified with coupled FBSDEs in different situations.

The organization of this paper is as follows: Section 2 will provide some assumptions about coupled FBSDEs. We use the classical Euler scheme for the forward process of FBSDE, which was applied to decoupled FBSDEs by Zhao et al. [23], to obtain the flexibility of the proposed calculation methods. For the backward process, the theta scheme will be used. In Section 3, by incorporating the Fourier variations, we expand a discretization of general coupled FBSDEs. In Section 4, three different approaches for dealing with the part of the coupling are described. In Section 5, six numerical experiments are discussed to verify the numerical methods. Finally, we summarize some conclusions in Section 6.

2. Assumptions and Discretization of the Coupled FBSDEs

Now, from this section, we studied the numerical calculation method of coupled FBSDEs (1) investigated by the following:

$$\{\begin{array}{c}{X}_t=X_0+{{\displaystyle \int }}_0^t\mu \left(s,X_s,Y_s,Z_s\right)ds+{{\displaystyle \int }}_0^t\Sigma \left(s,X_s,Y_s\right)d{W}_s\hfill \\ Y_t=Y_T+{{\displaystyle \int }}_t^{T}f\left(s,X_s,Y_s,Z_s\right)ds-{{\displaystyle \int }}_t^T{Z}_s^{T}d{W}_s\hfill \end{array}$$

(1)

where $X_t={\left(X_t^1,\cdots ,X_t^l\right)}^T\in {ℝ}^l,\left(0\le t\le T\right)$ is an $\iota$-dimensional forward stochastic differential equation (FSDE), and $Y_t\in ℝ,\left(0\le t\le T\right)$ is a one-dimensional backward stochastic differential equation(BSDE). $X_0=x_0$ and $Y_T=g\left(X_T\right)$ are the initial condition of the FSDE and the end condition of the BSDE, respectively. $Z_t^T={\left(Z_{t,1},\cdots ,Z_{t,k}\right)}^T\in {ℝ}^k$ is a control process. $\mu \left(t,X_t,Y_t,Z_t\right)=\left({\mu }_1\left(t,X_t^1,Y_t,Z_t\right),\cdots ,{\mu }_n\left(t,X_t^l,Y_t,Z_t\right)\right)$ are the drift terms of the FSDE, $\Sigma \left(t,X_t,Y_t\right)=\left({\Sigma }_{i,j}\left(t,X_t,Y_t\right)\right),1\le i\le l,1\le j\le k$ are diffusion terms of the FSDE, and $W_t={\left(W_t^1,\cdots ,W_t^k\right)}^T$ is defined on the probability space $\left(\mathsf{\Omega },\mathbb{F},{\left({ℱ}_t\right)}_{0\le t\le T},\mathbb{P}\right)$ with ${ℱ}_t$ the natural $\sigma$-algebra generated by $W_t$, which is a standard $k$-dimensional normal Brownian motion process and satisfies the usual conditions. Here, $f\left(t,X_t,Y_t,Z_t\right)$ 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.

2.1. Definitions and Assumptions

Following the survey papers [18,24,25], we provide the following standard definitions and assumptions. We define ${ℳ}^2\left(0,T;{ℝ}^l\right)$ as the set of ${ℱ}_T$-adapted processes $v(·):\Omega \to {ℝ}^l$ that are square integrable, i.e.,

$$E\left[{{\displaystyle \int }}_0^T{\left|v\left(s\right)\right|}^{2}ds\right]<\infty$$

where $|·|$ is the standard Euclidean norm in the Euclidean space. Denote ${\mathbb{L}}_T^1\left({ℝ}^l\right)$ 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 ${ℳ}^2\left(0,T;{ℝ}^l\right)$. 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 $\mu ,\Sigma ,f,$ and $g$ are uniformly Lipschitz in $\left(x,y,z\right)$, for $t\in \left[0,T\right]$. For example, $\mu$ should satisfy

$$\left|\mu \left(t,x_1,y_1,z_1\right)-\mu \left(t,x_2,y_2,z_2\right)\right|\le L^{Lipz}\left(\left|x_1-x_2\right|+\left|y_1-y_2\right|+\left|z_1-z_2\right|\right),$$

where $t\in \left[0,T\right]$ and $L^{Lipz}>0$ is the Lipschitz constant.

(2)

The functions $\mu$, $f,$ and $g$ satisfy the linear growth conditions as follows:

$$\left|f\left(t,x,y,z\right)\right|\le k_1\left(1+\left|y\right|+\left|z\right|\right),$$

$$\left|\mu \left(t,x,y,z\right)\right|\le k_2,$$

$$\left|g\left(x\right)\right|\le k_3,$$

where $k_1,k_2,k_3>0$ are constants.

(3)

The diffusion $\Sigma$ has bounded second derivatives with positive lower and upper bounds.

(4)

$V\left(x\right)\in {\mathbb{L}}_T^1\left({ℝ}^l\right)$ and $D_{x}V\left(x\right)\in {\mathbb{L}}_T^1\left({ℝ}^l\right)$, where $V\left(x\right)$ is a function of $x$ and $D_x$ is the weak derivative symbol for $x$.

Assumptions (1)–(4) can ensure that the solution $\left(x_t,y_t,z_t\right)$ to the coupled FBSDE (1) is existential and unique. If the function $V\left(x\right)$ or $D_{x}V\left(x\right)$ 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 $\theta$-scheme.

2.2. 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 $\mathsf{\Delta }$ of time grid $0=t_0<t_1<t_2<\cdots <t_M=T$ with time steps $\Delta t_{m+1}:t_{m+1}-t_m$, for $m=1,2,\cdots ,M-1$. For illustrative convenience, we denote $X_m=X_{t_m},Y_m=Y_{t_m},Z_m=Z_{t_m}$ and $\Delta 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:

$$X_{m+1}=X_m+{{\displaystyle \int }}_{t_m}^{t_{m+1}}\mu \left(s,X_s,Y_s,Z_s\right)ds+{{\displaystyle \int }}_{t_m}^{t_{m+1}}\Sigma \left(s,X_s,Y_s\right)d{W}_s.$$

(2)

The discretization of the FSDE is:

$$X_{m+1}^{\Delta }=X_m^{\Delta }+\mu \left(t_m,X_m^{\Delta },Y_m^{\Delta },Z_m^{\Delta }\right)\Delta t_{m+1}+\Sigma \left(t_m,X_m^{\Delta },Y_m^{\Delta }\right)\Delta W_{m+1}$$

(3)

The Equation (3) can also be explained as a left point approximation of Equation (2).

For the BSDE, we obtain:

$$Y_{m+1}=Y_m+{{\displaystyle \int }}_{t_m}^{t_{m+1}}f\left(s,X_s,Y_s,Z_s\right)ds-{{\displaystyle \int }}_{t_m}^{t_{m+1}}{Z}_s^{T}d{W}_s.$$

(4)

$Y_t$ is an ${ℱ}_t$-adapted process, so for both sides of Equation (4), we should adopt conditional expectations and use the $\theta$-scheme method to approximate the integral:

$$\begin{array}{c}{Y}_m=E_m^x\left[Y_{m+1}\right]+{{\displaystyle \int }}_{t_m}^{t_{m+1}}{E}_m^x\left[f\left(s,X_s,Y_s,Z_s\right)\right]ds\\ \approx E_m^x\left[Y_{m+1}\right]+{\theta }_{1}f\left(t_m,X_m,Y_m,Z_m\right)\Delta t_{m+1}+\left(1-{\theta }_1\right)\\ E_m^x\left[f\left(t_{m+1},X_{m+1},Y_{m+1},Z_{m+1}\right)\right]\Delta t_{m+1},\end{array}$$

(5)

where ${\theta }_1\in \left[0,1\right],$ and $E_m^x[·]$ denotes $E_m^x\left[·|X_m^{\Delta }=x\right]$. We first multiply both sides of Equation (4) by $\Delta W_{m+1}$. Then, also taking conditional expectations, similar to Equation (5), we get:

$$\begin{array}{c}0=E_m^x\left[Y_{m+1}\Delta W_{m+1}\right]+{{\displaystyle \int }}_{t_m}^{t_{m+1}}{E}_m^x\left[f\left(s,X_s,Y_s,Z_s\right)\Delta W_{m+1}\right]ds-{{\displaystyle \int }}_{t_m}^{t_{m+1}}{E}_m^x\left[Z_s^T\right]ds\\ \approx E_m^x\left[Y_{m+1}\Delta W_{m+1}\right]+\left(1-{\theta }_2\right)\left[f\left(s,X_s,Y_s,Z_s\right)\Delta W_{m+1}\right]\Delta t_{m+1}\\ -{\theta }_2{Z}_m^T\Delta t_{m+1}-\left(1-{\theta }_2\right)E_m^x\left[Z_{m+1}^T\right]\Delta t_{m+1}.\end{array}$$

(6)

The terminal values $Y_M$ and $Z_M$ can be obtained by using $X_M^{\Delta }$ because of a result of the Nonlinear Feynman-Kac Theorem [25,26]. Then we have:

$$Y_M=g\left(X_M\right)\hspace{1em}\hspace{1em}{Z}_M={\Sigma }^T\left(t_M,X_M,Y_M\right)\nabla g\left(X_M\right)={\Sigma }^T\left(t_M,X_M,Y_M\right)D_{x}g\left(X_M\right)$$

with $\nabla$ as a gradient operator. The above equations are the basis of the following numerical methods. Since $X^{\Delta }$ in Equation (3) is a Markov process, we use a discrete-time approximation $\left(y\left(t_m,X_m^{\Delta }\right),z\left(t_m,X_m^{\Delta }\right)\right)$ for $\left(Y_m,Z_m\right)$ in the BSDE. Then we have:

$$y\left(t_M,x\right)=g\left(x\right),\hspace{1em}y\left(t_M,x\right)={\Sigma }^T\left(t_M,x,y\left(t_M,x\right)\right)D_{x}g\left(x\right),$$

(7)

$$\begin{array}{cc}\hfill y\left(t_m,x\right)& =E_m^x\left[y\left(t_{m+1},X_{m+1}^{\Delta }\right)\right]+{\theta }_{1}f\left(t_m,x,y\left(t_m,x\right),z\left(t_m,x\right)\right)\Delta t_{m+1}\hfill \\ & +\left(1+{\theta }_1\right)E_m^x\left[f\left(t_{m+1},X_{m+1}^{\Delta },y\left(t_{m+1},X_{m+1}^{\Delta }\right),z\left(t_{m+1},X_{m+1}^{\Delta }\right)\right)\right]\Delta t_{m+1},\hfill \end{array}$$

(8)

$$\begin{array}{cc}\hfill z\left(t_m,x\right)& =-\frac{1-{\theta }_2}{{\theta }_2}{E}_m^x\left[z\left(t_{m+1},X_{m+1}^{\Delta }\right)\right]+\frac{1}{{\theta }_2}{E}_m^x\left[y\left(t_{m+1},X_{m+1}^{\Delta }\right)\Delta W_{m+1}\right]\frac{1}{\Delta t_{m+1}}\hfill \\ & +\frac{1-{\theta }_2}{{\theta }_2}{E}_m^x\left[f\left(t_{m+1},X_{m+1}^{\Delta },y\left(t_{m+1},X_{m+1}^{\Delta }\right),z\left(t_{m+1},X_{m+1}^{\Delta }\right)\right)\Delta W_{m+1}\right]{,}_{}\hfill \end{array}$$

(9)

where $m=M-1,\cdots ,0,{\theta }_1\in \left[0,1\right]$ and ${\theta }_2\in \left(0,1\right]$. In Equation (8), when setting ${\theta }_1=0$, we will get an explicit scheme for $y\left(t_m,x\right)$. Based on observations of Equations (8) and (9), the numerical iteration of the BSDE lies in conditional expectations. In the next section, we will introduce the fractional fast Fourier transform (FrFFT) method to transfer the iteration to the Fourier space so that Equations (8) and (9) iterate smoothly.

3. Fractional Fast Fourier Transform Method

In this section, we derive a Fourier transform for Equations (7)–(9) to transfer the iteration to the Fourier space. Because the three existing algorithms differ from this scheme in the next section, the final FrFFT method used in the three algorithms is different in this section.

The conditional expectations are formed like the two equations below:

$$E\left[V\left(t_{m+1},X_{m+1}^{\Delta }\right)\right],$$

$$E\left[V\left(t_{m+1},X_{m+1}^{\Delta }\right)\Delta W_{m+1}\right].$$

Since the evolution of $X_{m+1}^{\Delta }$ from Equation (3) is known, let $m=:\mu \left(t_m,x,y\left(t_m,x\right),z\left(t_m,x\right)\right)$ and $s=:\Sigma (t_m,x,y\left(t_m,x\right)$). Then, we have:

$$E_m^x\left[V\left(t_{m+1},X_{m+1}^{\Delta }\right)\Delta W_{m+1,j}\right]=\Delta t_{m+1}{s}_{j,(·)}^T{E}_m^x\left[D_{x}V\left(t_{m+1},X_{m+1}^{\Delta }\right)\right]{,}_{}$$

where $s_{j,(·)}^T$ is the $j^{th}$ row vector of $s$.

For convenience and clarity, we consider the following substitutions:

$${𝒫}_{\Delta }:=E_m^x[·]{,}_{}$$

$$\mu :=\mu \left(t_m,x,y\left(t_m,x\right),z\left(t_m,x\right)\right),$$

$$\Sigma :=\mu \left(t_m,x,y\left(t_m,x\right)\right),$$

$$f\left(t_m,x\right):=f\left(t_m,x,y\left(t_m,x\right),z\left(t_m,x\right)\right),$$

$$\gamma \left(t_{m+1},X_{m+1}^{\Delta }\right):={𝒫}_{\Delta }y\left(t_{m+1},X_{m+1}^{\Delta }\right)+\left(1+{\theta }_1\right){𝒫}_{\Delta }f\left(t_{m+1},X_{m+1}^{\Delta }\right)\Delta t_{m+1}.$$

(10)

We rewrite Equations (7)–(9):

$$y\left(t_M,x\right)=g\left(x\right),\hspace{1em}\hspace{1em}z\left(t_M,x\right)={\Sigma }^T{D}_{x}g\left(x\right),$$

(11)

$$y\left(t_m,x\right)=\gamma \left(t_{m+1},X_{m+1}^{\Delta }\right)+{\theta }_{1}f\left(t_m,x\right)\Delta t_{m+1},$$

(12)

$$\begin{array}{c}z\left(t_m,x\right)=-\frac{1-{\theta }_2}{{\theta }_2}{𝒫}_{\Delta }z\left(t_{m+1},X_{m+1}^{\Delta }\right)+\frac{1}{{\theta }_2}{\Sigma }^T{𝒫}_{\Delta }y\left(t_{m+1},X_{m+1}^{\Delta }\right)\\ +\frac{1-{\theta }_2}{{\theta }_2}{\Sigma }^T{𝒫}_{\Delta }f\left(t_{m+1},X_{m+1}^{\Delta }\right)\Delta t_{m+1}.\end{array}$$

(13)

We consider the Fourier transform of a function $V\left(x\right)$ as:

$${\mathbb{F}}_{\iota }V\left(\xi \right)={{\displaystyle \int }}_{{ℝ}^{\iota }}^{}V\left(x\right)e^{i{\xi }^{T}x}dx,\hspace{1em}\hspace{1em}\xi \in {ℝ}^{\iota }$$

and the Fourier inverse transform in the same way:

$$V\left(x\right)={\mathbb{F}}_{\iota }^{-1}\widehat{V}\left(x\right)=\frac{1}{{\left(2\pi \right)}^{\iota }}{{\displaystyle \int }}_{{ℝ}^{\iota }}^{}\widehat{V}\left(\xi \right)e^{-i{\xi }^{T}x}d\xi ,$$

where $\widehat{V}\left(\xi \right):={\mathbb{F}}_{\iota }V\left(\xi \right)$ and $V\left(x\right)\in \left\{y\left(t_m,x\right),z\left(t_m,x\right),f\left(t_m,x\right)\right\}$. We then need to numerically calculate the above two equations, truncating the integral with a sufficiently large region $\prod \in {ℝ}^{\iota }$ and ${\displaystyle \tilde{\prod }}\in {ℝ}^{\iota }$ with:

$$\prod =\left[a,b\right]:=\left[a_1,b_1\left]\times \cdots \right[a_{\iota },b_{\iota }\right],$$

$$\tilde{\prod }=\left[\tilde{a},\tilde{b}\right]:=\left[{\tilde{a}}_1,{\tilde{b}}_1\right]\times \cdots \left[{\tilde{a}}_{\iota },{\tilde{b}}_{\iota }\right].$$

We get the approximation of integrals by the midpoint rule:

$$\widehat{V}\left(\xi \right)\approx h_x{\displaystyle \sum }_{n_1=0}^{N-1}\cdots {\displaystyle \sum }_{n_{\iota }=0}^{N-1}{e}^{i{\xi }^T{x}_n}V\left(x_n\right),$$

(14)

$$V\left(x\right)\approx \frac{1}{{\left(2\pi \right)}^{\iota }}{h}_{\xi }{\displaystyle \sum }_{{\tilde{n}}_1=0}^{N-1}\cdots {\displaystyle \sum }_{{\tilde{n}}_{\iota }=0}^{N-1}{e}^{-i{\xi }_{\tilde{n}}x}\widehat{V}\left({\xi }_{\tilde{n}}\right),$$

(15)

where $h_{x,j}=\frac{b_j-a_j}{N},h_x={\prod }_{j=1}^{\iota }{h}_{x,j},h_{\xi ,j}=\frac{{\tilde{b}}_j-{\tilde{a}}_j}{N},h_{\xi }={\prod }_{j=1}^{\iota }{h}_{\xi ,j}$, and

$$x_n={\left(a_j+\left(n_j+\frac{1}{2}\right)h_{x,j}\right)}_{1\times \iota },\hspace{1em}\hspace{1em}j=1,\cdots ,\iota ,\hspace{1em}\hspace{1em}{n}_j\in 0,\cdots ,N-1;$$

$${\xi }_{\tilde{n}}={\left({\tilde{a}}_j+\left({\tilde{n}}_j+\frac{1}{2}\right)h_{\xi ,j}\right)}_{1\times \iota },\hspace{1em}\hspace{1em}j=1,\cdots ,\iota ,\hspace{1em}\hspace{1em}{\tilde{n}}_j\in 0,\cdots ,N-1.$$

We compute $\widehat{V}\left(\xi \right)$ by Equation (14) on $\mathsf{\xi }$-grid and $V\left(x\right)$ 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

$$x_n=a+n{h}_x,\hspace{1em}\hspace{1em}{h}_x=\frac{b-a}{N},\hspace{1em}\hspace{1em}n=0,\cdots ,N-1,$$

$${\xi }_{\tilde{n}}=\tilde{a}+\tilde{n}{h}_{\xi },\hspace{1em}\hspace{1em}{h}_{\xi }=\frac{\tilde{b}-\tilde{a}}{N},\hspace{1em}\hspace{1em}\tilde{n}=0,\cdots ,N-1,$$

then

$$\begin{array}{cc}\hfill \widehat{V}\left({\xi }_{\tilde{n}}\right)& \approx h_x{\displaystyle {\displaystyle \sum }_{n=0}^{N-1}}{e}^{i{\xi }_{\tilde{n}}{x}_n}V\left(x_n\right)\hfill \\ & =h_x{\displaystyle {\displaystyle \sum }_{n=0}^{N-1}}{e}^{i\left(\tilde{a}+\tilde{n}{h}_{\xi }\right)\left(a+n{h}_x\right)}V\left(x_n\right)\hfill \\ & =h_x{e}^{i\left(a\tilde{a}+a\tilde{n}{h}_{\xi }\right)}{\displaystyle {\displaystyle \sum }_{n=0}^{N-1}}{e}^{i\tilde{a}n{h}_x}{e}^{in\tilde{n}{h}_{\xi }{h}_x}V\left(x_n\right)\hfill \\ & =h_x{e}^{i\left(a\tilde{a}+a\tilde{n}{h}_{\xi }\right)}{\displaystyle {\displaystyle \sum }_{n=0}^{N-1}}{e}^{-2\pi in\tilde{n}\tilde{\vartheta }}{V}_n^{*},\hfill \end{array}$$

and

$$\begin{array}{cc}\hfill V\left(x_n\right)& \approx \frac{1}{2\pi }{h}_{\xi }{\displaystyle {\displaystyle \sum }_{\tilde{n}=0}^{N-1}}{e}^{-i{\xi }_{\tilde{n}}{x}_n}\widehat{V}\left({\xi }_{\tilde{n}}\right)\hfill \\ & =\frac{1}{2\pi }{h}_{\xi }{\displaystyle {\displaystyle \sum }_{\tilde{n}=0}^{N-1}}{e}^{-i\left(\tilde{a}+\tilde{n}{h}_{\xi }\right)\left(a+n{h}_x\right)}\widehat{V}\left({\xi }_{\tilde{n}}\right)\hfill \\ & =\frac{1}{2\pi }{h}_{\xi }{e}^{-i\left(\tilde{a}a+\tilde{a}n{h}_x\right)}{\displaystyle {\displaystyle \sum }_{\tilde{n}=0}^{N-1}}{e}^{-ia\tilde{n}{h}_{\xi }}{e}^{-i\tilde{n}{h}_{\xi }{h}_x}\widehat{V}\left({\xi }_{\tilde{n}}\right)\hfill \\ & =\frac{1}{2\pi }{h}_{\xi }{e}^{-i\left(\tilde{a}a+\tilde{a}n{h}_x\right)}{\displaystyle {\displaystyle \sum }_{\tilde{n}=0}^{N-1}}{e}^{-2\pi i\tilde{n}n\vartheta }{\widehat{V}}_{\tilde{n}}^{*},\hfill \end{array}$$

where $\tilde{\vartheta }=-\frac{h_{\xi }{h}_x}{N},\vartheta =\frac{h_{\xi }{h}_x}{N},V_n^{*}=e^{i\tilde{a}n{h}_x}V\left(x_n\right)$, and ${\widehat{V}}_{\tilde{n}}^{*}=e^{-ia\tilde{n}{h}_{\xi }}\widehat{V}\left({\xi }_{\tilde{n}}\right)$. The appendix in [7] has been demonstrated by the ι-multidimensional case of Equation (14). We will generalize Equation (15) to a $\iota$-multidimensional situation. Consider the $\iota$-dimensional integral,

$$\begin{array}{cc}\hfill V\left(x_n\right)& =\frac{1}{{\left(2\pi \right)}^{\iota }}{{\displaystyle \int }}_{\left[\tilde{a},\tilde{b}\right]}^{}{e}^{-i{\xi }^T{x}_n}\widehat{V}\left(\xi \right)d\xi \hfill \\ & \approx \frac{1}{{\left(2\pi \right)}^{\iota }}{h}_{{\xi }_1}\cdots h_{{\xi }_2}{\displaystyle {\displaystyle \sum }_{{\tilde{n}}_1=0}^{N-1}}\cdots {\displaystyle {\displaystyle \sum }_{{\tilde{n}}_{\iota }=0}^{N-1}}{e}^{-i\left({\xi }_{{\tilde{n}}_1}^1{x}_{n_1}^1+\cdots +{\xi }_{{\tilde{n}}_{\iota }}^{\iota }{x}_{n_{\iota }}^{\iota }\right)}\widehat{V}\left({\xi }_{{\tilde{n}}_1}^1,\cdots ,{\xi }_{{\tilde{n}}_{\iota }}^{\iota }\right),\hfill \end{array}$$

where

$$x_{n_j}^j=a_j+n_j{h}_{x_j},\hspace{1em}\hspace{1em}{h}_{x_j}=\frac{b_j-a_j}{N},\hspace{1em}{n}_j=0,\cdots ,N-1,$$

$${\xi }_{{\tilde{n}}_j}^j={\tilde{a}}_j+{\tilde{n}}_j{h}_{{\xi }_j},\hspace{1em}\hspace{1em}{h}_{{\xi }_j}=\frac{{\tilde{b}}_j-{\tilde{a}}_j}{N},\hspace{1em}\hspace{1em}{\tilde{n}}_j=0,\cdots ,N-1.$$

To review this $\iota$-dimensional sum, we can carry out the caculation as follows.

$$V\left(x_{n_1}^1,\cdots ,x_{n_{\iota }}^{\iota }\right)\approx \frac{1}{{\left(2\pi \right)}^{\iota }}{h}_{{\xi }_1}{\displaystyle \sum }_{{\tilde{n}}_1=0}^{N-1}{e}^{-i{\xi }_{{\tilde{n}}_1}^1{x}_{n_1}^1}G\left({\xi }_{{\tilde{n}}_1}^1,x_{n_2}^2,\cdots ,x_{n_{\iota }}^{\iota }\right),$$

(16)

where

$$G\left({\xi }_{{\tilde{n}}_1}^1,x_{n_2}^2,\cdots ,x_{n_{\iota }}^{\iota }\right)=h_{{\xi }_2}\cdots h_{{\xi }_{\iota }}{\displaystyle \sum }_{{\tilde{n}}_2=0}^{N-1}\cdots {\displaystyle \sum }_{{\tilde{n}}_{\iota }=0}^{N-1}{e}^{-i\left({\xi }_{{\tilde{n}}_2}^2{x}_{n_2}^2+\cdots +{\xi }_{{\tilde{n}}_{\iota }}^{\iota }{x}_{n_{\iota }}^{\iota }\right)}\widehat{V}\left({\xi }_{{\tilde{n}}_1}^1,\cdots ,{\xi }_{{\tilde{n}}_{\iota }}^{\iota }\right).$$

(17)

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 $\left(\widehat{y}\left(t_m,\xi \right),\widehat{z}\left(t_m,\xi \right)\right)$ is represented by

$${\mathbb{F}}_{\iota }\left({\Sigma }^{\tau }{𝒫}_{\Delta }{D}_{x}V\left(x\right)\right)\left(\xi \right),\hspace{1em}\hspace{1em}{\mathbb{F}}_{\iota }\left({𝒫}_{\Delta }V\left(x\right)\right)\left(\xi \right)$$

for some function $V\left(x\right)$.

To approximate the two equations, we introduce Lemma 1 and Theorem 1.

Lemma 1 ([29]).

If $V\in L^1\left({ℝ}^n\right)$, then

$${𝒫}_{\Delta }V\left(x\right)=\frac{1}{{\left(2\pi \right)}^{\iota }}{{\displaystyle \int }}_{{ℝ}^{\iota }}^{}{e}^{-i{\xi }^{T}x}{\varphi }_{\Delta }\left(-\xi ;x\right)\widehat{V}\left(\xi \right)d\xi ,$$

where${\varphi }_{\Delta }(·)$is the characteristic function of$X_{m+1}^{\Delta }-x$and${\varphi }_{\Delta }\left(-\xi ;x\right)$can be denoted by

$${\varphi }_{\Delta }\left(-\xi ;x\right)=exp\left(-\Delta t_m\left(\frac{1}{2}{\xi }^T\Sigma {\Sigma }^T\xi +i{\xi }^T\mu \right)\right).$$

(18)

Theorem 1.

If $V\in L^1\left({ℝ}^n\right)$, then

$${\mathbb{F}}_{\iota }\left({𝒫}_{\Delta }V\left(x\right)\right)\left(\xi \right)={\varphi }_{\Delta }\left(-\xi ;x\right)\widehat{V}\left(\xi \right),$$

$${\mathbb{F}}_{\iota }\left({\Sigma }^T{𝒫}_{\Delta }{D}_{x}V\left(x\right)\right)\left(\xi \right)=-i\xi {\Sigma }^T{\varphi }_{\Delta }\left(-\xi ;x\right)\widehat{V}\left(\xi \right).$$

Proof. 

For the first equation, since

$$\begin{array}{cc}\hfill {𝒫}_{\Delta }V\left(x\right)& =\frac{1}{{\left(2\pi \right)}^{\iota }}{{\displaystyle \int }}_{{ℝ}^{\iota }}^{}{e}^{-i{\xi }^{T}x}{\varphi }_{\Delta }\left(-\xi ;x\right)\widehat{V}\left(\xi \right)d\xi \hfill \\ & ={\mathbb{F}}_{\iota }^{-1}\left[{\varphi }_{\Delta }\left(-\xi ;x\right)\widehat{V}\left(\xi \right)\right]\left(\mathrm{x}\right),\hfill \end{array}$$

then we have

$${\mathbb{F}}_{\iota }\left({𝒫}_{\Delta }V\left(x\right)\right)\left(\xi \right)={\mathbb{F}}_{\iota }{\mathbb{F}}_{\iota }^{-1}\left[{\varphi }_{\Delta }\left(-\xi ;x\right)\widehat{V}\left(\xi \right)\right]={\varphi }_{\Delta }\left(-\xi ;x\right)\widehat{V}\left(\xi \right).$$

For the second equation:

$$\begin{array}{cc}\hfill D_{x}V\left(x\right)& =D_x\left[\frac{1}{{\left(2\pi \right)}^{\iota }}{{\displaystyle \int }}_{{ℝ}^{\iota }}^{}\widehat{V}\left(\xi \right)e^{-i{\xi }^{T}x}d\xi \right]\hfill \\ & =\frac{1}{{\left(2\pi \right)}^{\iota }}{{\displaystyle \int }}_{{ℝ}^{\iota }}^{}{D}_x\left[\widehat{V}\left(\xi \right)e^{-i{\xi }^{T}x}\right]d\xi \hfill \\ & =\frac{1}{{\left(2\pi \right)}^{\iota }}{{\displaystyle \int }}_{{ℝ}^{\iota }}^{}-i\xi \widehat{V}\left(\xi \right)e^{-i{\xi }^{T}x}d\xi .\hfill \end{array}$$

Let ${\widehat{V}}_1\left(\xi \right)=-i\xi \widehat{V}\left(\xi \right)$. Then we get $D_{x}V\left(x\right)=\frac{1}{{\left(2\pi \right)}^{\iota }}{{\displaystyle \int }}_{{ℝ}^{\iota }}^{}{\widehat{V}}_1\left(\xi \right)e^{-i{\xi }^{T}x}d\xi =V_1\left(x\right)$. Since

$$\begin{array}{cc}\hfill {\Sigma }^T{𝒫}_{\Delta }{D}_{x}V\left(x\right)& ={\Sigma }^T{𝒫}_{\Delta }{V}_1\left(x\right)\hfill \\ & =\frac{1}{{\left(2\pi \right)}^{\iota }}{{\displaystyle \int }}_{{ℝ}^{\iota }}^{}{e}^{-i{\xi }^{T}x}{\Sigma }^T{\varphi }_{\Delta }\left(-\xi ;x\right){\widehat{V}}_1\left(\xi \right)d\xi \hfill \\ & ={\mathbb{F}}_{\iota }^{-1}\left({\Sigma }^T{\varphi }_{\Delta }\left(-\xi ;x\right)\widehat{V}\left(\xi \right)\right)\left(\mathrm{x}\right),\hfill \end{array}$$

we have

$${\mathbb{F}}_{\iota }\left({\Sigma }^T{\varphi }_{\Delta }\left(-\xi ;x\right)V\left(x\right)\right)\left(\xi \right)={\mathbb{F}}_{\iota }{\mathbb{F}}_{\iota }^{-1}\left({\Sigma }^T{\varphi }_{\Delta }\left(-\xi ;x\right)\widehat{V}\left(\xi \right)\right)=-i\xi {\Sigma }^T{\varphi }_{\Delta }\left(-\xi ;x\right)\widehat{V}\left(\xi \right).$$

□

We now have the approximation of $\left(\widehat{y}\left(t_m,\xi \right),\widehat{z}\left(t_m,\xi \right)\right)$.

$$\widehat{y}\left(t_M,\xi \right)={\mathbb{F}}_{\iota }\left(g\left(x\right)\right)\left(\xi \right),\hspace{1em}\hspace{1em}\widehat{z}\left(t_M,\xi \right)={\Sigma }^T{\mathbb{F}}_{\iota }\left(D_{x}g\left(x\right)\right)\left(\xi \right),$$

(19)

$$\widehat{y}\left(t_m,\xi \right)=\widehat{\gamma }\left(t_{m+1},\xi \right)+{\theta }_1\widehat{f}\left(t_m,\xi \right)\Delta t_{m+1},$$

(20)

$$\begin{array}{cc}\hfill \widehat{z}\left(t_m,\xi \right)=& -\frac{1-{\theta }_2}{{\theta }_2}{\varphi }_{\Delta }\left(-\xi ;x\right)\widehat{z}\left(t_{m+1},\xi \right)-i\frac{1}{{\theta }_2}{\Sigma }^T{\varphi }_{\Delta }\left(-\xi ;x\right)\widehat{y}\left(t_{m+1},\xi \right)\hfill \\ & -i\frac{1-{\theta }_2}{{\theta }_2}{\Sigma }^T{\varphi }_{\Delta }\left(-\xi ;x\right)\widehat{f}\left(t_{m+1},\xi \right)\Delta t_{m+1},\hfill \end{array}$$

(21)

where

$$\widehat{\gamma }\left(t_{m+1},\xi \right):={\varphi }_{\Delta }\left(-\xi ;x\right)\widehat{y}\left(t_{m+1},\xi \right)+\left(1-{\theta }_1\right){\varphi }_{\Delta }\left(-\xi ;x\right)\widehat{f}\left(t_{m+1},\xi \right)\Delta t_{m+1}.$$

(22)

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.

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

4.1. Explicit Algorithm

We apply the technique from [30]. Suppose the time step $\Delta 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_{m+1}=x+\mu \left(t_{m+1},x,y\left(t_{m+1},x\right),z\left(t_{m+1},x\right)\right)\Delta t_{m+1}+\Sigma \left(t_{m+1},x,y\left(t_{m+1},x\right)\right)\Delta 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}^{\Delta }-x$ are also different. However, the differences are minimal because we just need to replace $t_m$ in Equation (18) with $t_{m+1}$:

$$\begin{array}{c}{\varphi }_{\Delta }\left(-\xi ;x\right)=exp(-\Delta t_m(\frac{1}{2}{\xi }^T\Sigma \left(t_{m+1},x,y\left(t_{m+1},x\right)\right){\Sigma }^T\left(t_{m+1},x,y\left(t_{m+1},x\right)\right)\xi \\ +i{\xi }^T\mu \left(t_{m+1},x,y\left(t_{m+1},x\right),z\left(t_{m+1},x\right)\right)\left)\right).\end{array}$$

(24)

Furthermore, we have to replace $t_m$ with $t_{m+1}$ in Equations (10)–(13) and Equations (19)–(22).

$$\mathrm{y}\left(t_M,x\right)=g\left(x\right),\hspace{1em}\hspace{1em}z\left(t_M,x\right)={\Sigma }^T\left(t_M,x,y\left(t_M,x\right)\right)D_{x}g\left(x\right),$$

(25)

$$y\left(t_m,x\right)=\gamma \left(t_{m+1},X_{m+1}\right)+\frac{1}{{\theta }_1}f\left(t_m,x\right)\Delta t_{m+1},$$

(26)

$$\begin{array}{cc}\hfill z\left(t_m,x\right)=& -\frac{1-{\theta }_2}{{\theta }_2}{𝒫}_{\Delta }z\left(t_{m+1},X_{m+1}\right)+\frac{1}{{\theta }_2}{\Sigma }^T\left(t_{m+1},x,y\left(t_{m+1},x\right)\right){𝒫}_{\Delta }y\left(t_{m+1},X_{m+1}\right)\hfill \\ & +\frac{1-{\theta }_2}{{\theta }_2}{\Sigma }^T\left(t_{m+1},x,y\left(t_{m+1},x\right)\right){𝒫}_{\Delta }f\left(t_{m+1},X_{m+1}\right)\Delta t_{m+1},\hfill \end{array}$$

(27)

where

$$\gamma \left(t_{m+1},X_{m+1}\right)={𝒫}_{\Delta }y\left(t_{m+1},X_{m+1}\right)+\left(1-{\theta }_1\right){𝒫}_{\Delta }f\left(t_{m+1},X_{m+1}\right)\Delta t_{m+1},$$

(28)

and

$$\widehat{y}\left(t_M,\xi \right)={\mathbb{F}}_{\iota }\left(g\left(x\right)\right)\left(\xi \right),\hspace{1em}\hspace{1em}\widehat{z}\left(t_M,\xi \right)={\Sigma }^T\left(t_M,x,y\left(t_M,x\right)\right){\mathbb{F}}_{\iota }\left(D_{x}g\left(x\right)\right)\left(\xi \right),$$

(29)

$$\widehat{y}\left(t_m,\xi \right)=\widehat{\gamma }\left(t_{m+1},\xi \right)+{\theta }_1\widehat{f}\left(t_m,\xi \right)\Delta t_{m+1}.$$

(30)

$$\begin{array}{c}\widehat{z}\left(t_m,\xi \right)=-\frac{1-{\theta }_2}{{\theta }_2}{\varphi }_{\Delta }\left(-\xi ;x\right)\widehat{z}\left(t_{m+1},\xi \right)-\\ i\frac{1}{{\theta }_2}{\Sigma }^T\left(t_{m+1},x,y\left(t_{m+1},x\right)\right){\varphi }_{\Delta }\left(-\xi ;x\right)\widehat{y}\left(t_{m+1},\xi \right)\\ -i\frac{1-{\theta }_2}{{\theta }_2}{\Sigma }^T\left(t_{m+1},x,y\left(t_{m+1},x\right)\right){\varphi }_{\Delta }\left(-\xi ;x\right)\widehat{f}\left(t_{m+1},\xi \right)\Delta t_{m+1},\end{array}$$

(31)

where

$$\widehat{\gamma }\left(t_{m+1},\xi \right)={\varphi }_{\Delta }\left(-\xi ;x\right)\widehat{y}\left(t_{m+1},\xi \right)+\left(1-{\theta }_1\right){\varphi }_{\Delta }\left(-\xi ;x\right)\widehat{f}\left(t_{m+1},\xi \right)\Delta t_{m+1}.$$

(32)

Based on the above analysis, we give the following detailed Algorithm 1 to implement the numerical solution to Equation (1):

Algorithm 1. Explicit algorithm

Initialization:

1: Calculate $\widehat{y}\left(t_m,\xi \right),\widehat{z}\left(t_m,\xi \right)$ and $\widehat{f}\left(t_m,\xi \right)$ applying Equation (14) by FrFFT on the $\xi$-grid.

Main loop:

2: for $m=M-1$ to 1 do

3:  Calculate $\widehat{z}\left(t_m,\xi \right)$ and $\widehat{\gamma }\left(t_{m+1},\xi \right)$ via Equations (31) and (32).

4:  Calculate $z\left(t_m,x\right)$ and $\gamma \left(t_{m+1},x\right)$ applying Equation (15) by FrFFT on the $x$-grid.

5:  Do Picard iteration $P$ times for Equation (26) to get y(tm, x) on the $x$-grid.

6:  Get $f\left(t_m,x\right)$ on the $x$-grid.

7:  if $m\ne 0$ then

8:   Calculate $\widehat{y}\left(t_m,\xi \right)$ and $\widehat{\gamma }\left(t_m,\xi \right)$ applying Equation (14) by FrFFT on the $\xi$-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$:

• Calculation of the function ${\varphi }_{\Delta }\left(-\xi ;x\right)$ on the $x$-grid in $O\left(DN\right)$ operations.
• Computation of $\widehat{z}\left(t_m,\xi \right),\widehat{\gamma }\left(t_m,\xi \right)$, and $\widehat{y}\left(t_m,\xi \right)$ on the $\xi$-grid in $O\left(D{N}^{\iota }\right)$ operations.
• Computation of $z\left(t_m,x\right)$ and $\gamma \left(t_m,x\right)$ on the $x$-grid in $O\left(D{N}^{\iota }log\left(N\right)\right)$ operations.
• Computation of $y\left(t_m,x\right)$ on the $x$-grid in $O\left(P{N}^{\iota }\right)$ operations.
• Computation of the cubic spline on the $x$-grid in $O\left(1\right)$ operations.

The overall complexity of the explicit algorithm is $O\left(M\left(DN+D{N}^{\iota }+D{N}^{\iota }log\left(N\right)+P{N}^{\iota }\right)\right)$, where $D$ is the number of components for $\widehat{z}\left(t_m,\xi \right)$.

Different: Each iteration is made on the $\xi$-grid. Finally, at $t_0$, the cubic spline is used to get the value at point $x_0$.

4.2. Local Algorithm

We consider the iteration of each time interval $\left[t_m,t_{m+1}\right)$ for $m=0,\cdots ,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,\cdots ,P_{local}$ the current local iteration. Furthermore, we denote by $X_m^{\Delta ,k},{\widehat{y}}^k\left(t_m,\xi \right),{\widehat{z}}^k\left(t_m,\xi \right),$ ${\widehat{f}}^k\left(t_m,\xi \right),{\widehat{\gamma }}^k\left(t_m,\xi \right),y^k\left(t_m,x\right),z^k\left(t_m,x\right),f^k\left(t_m,x\right),{\gamma }^k\left(t_m,x\right)$, and ${\gamma }^k\left(t_m,x\right),$ respectively, the value of $X_m^{\Delta },\widehat{y}\left(t_m,\xi \right),\widehat{z}\left(t_m,\xi \right),\widehat{f}\left(t_m,\xi \right),\widehat{\gamma }\left(t_m,\xi \right),y\left(t_m,x\right),z\left(t_m,x\right),f\left(t_m,x\right),$ and $\gamma \left(t_m,x\right)$ in iteration $k$. Finally, we denote the values by $y_{*}\left(t_m,x\right)$ and $z_{*}\left(t_m,x\right)$ 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\left(t_m,x\right)$ and $z_{*}^k\left(t_m,x\right)$ in local iteration. Then, Equation (18) will be rewritten by:

$$\begin{array}{c}{\varphi }_{\Delta }^{\mathrm{k}}\left(-\xi ;x\right)=exp(-\Delta t_m(\frac{1}{2}{\xi }^T\Sigma \left(t_{m+1},x,y_{*}^k\left(t_{m+1},x\right)\right){\Sigma }^T\left(t_{m+1},x,y_{*}^k\left(t_{m+1},x\right)\right)\xi +\\ i{\xi }^T\mu \left(t_{m+1},x,y_{*}^k\left(t_{m+1},x\right),z_{*}^k\left(t_{m+1},x\right)\right)\left)\right).\end{array}$$

(33)

However, we do not have functions with $y^k\left(t_m,x\right)$ and $z^k\left(t_m,x\right)$ at the iteration $k$ and time point $t_m$, so we apply the previous iteration for approximating Equation (33), i.e.,

$$\begin{array}{c}{\varphi }_{\Delta }^{\mathrm{k}}\left(-\xi ;x\right)\approx exp(-\Delta t_m(\frac{1}{2}{\xi }^T\Sigma \left(t_{m+1},x,y_{*}^{k-1}\left(t_{m+1},x\right)\right){\Sigma }^T\left(t_{m+1},x,y_{*}^{k-1}\left(t_{m+1},x\right)\right)\xi +\\ i{\xi }^T\mu \left(t_{m+1},x,y_{*}^{k-1}\left(t_{m+1},x\right),z_{*}^{k-1}\left(t_{m+1},x\right)\right)\left)\right).\end{array}$$

(34)

Furthermore, their derivations do not change.

$$y^{P_{local}}\left(t_M,x\right)=g\left(x\right),\hspace{1em}\hspace{1em}{z}^{P_{local}}\left(t_M,x\right)={\Sigma }^T\left(t_M,x,y_{*}^k\left(t_M,x\right)\right)D_{x}g\left(x\right),$$

(35)

$$y^k\left(t_m,x\right)={\gamma }^k\left(t_{m+1},X_{m+1}\right)+\frac{1}{{\theta }_1}{f}^k\left(t_m,x\right)\Delta t_{m+1},$$

(36)

$$\begin{array}{c}{z}^k\left(t_m,x\right)=-\frac{1-{\theta }_2}{{\theta }_2}{𝒫}_{\Delta }{z}^k\left(t_{m+1},X_{m+1}\right)+\frac{1}{{\theta }_2}{\Sigma }^T\left(t_{m+1},x,y_{*}^{k-1}\left(t_m,x\right)\right){𝒫}_{\Delta }{y}^k\left(t_{m+1},X_{m+1}\right)\\ +\frac{1-{\theta }_2}{{\theta }_2}{\Sigma }^T\left(t_m,x,y_{*}^{k-1}\left(t_m,x\right)\right){𝒫}_{\Delta }{f}^k\left(t_{m+1},X_{m+1}\right)\Delta t_{m+1},\end{array}$$

(37)

where

$${\gamma }^k\left(t_{m+1},X_{m+1}\right)={𝒫}_{\Delta }{y}^k\left(t_{m+1},X_{m+1}\right)+\left(1-{\theta }_1\right){𝒫}_{\Delta }{f}^k\left(t_{m+1},X_{m+1}\right)\Delta t_{m+1},$$

(38)

and

$$\begin{array}{cc}\hfill {\widehat{y}}^{P_{local}}\left(t_M,\xi \right)& ={\mathbb{F}}_{\iota }\left(g\left(x\right)\right)\left(\xi \right),\hspace{1em}\hspace{1em}{\widehat{z}}^{P_{local}}\left(t_M,\xi \right)\hfill \\ & ={\Sigma }^T\left(t_M,x,y_{*}^k\left(t_M,x\right)\right){\mathbb{F}}_{\iota }\left(D_{x}g\left(x\right)\right)\left(\xi \right),\hfill \end{array}$$

(39)

$${\widehat{y}}^k\left(t_m,\xi \right)={\widehat{\gamma }}^k\left(t_{m+1},\xi \right)+{\theta }_1{\widehat{f}}^k\left(t_m,\xi \right)\Delta t_{m+1},$$

(40)

$$\begin{array}{c}{\widehat{z}}^k\left(t_m,\xi \right)=-\frac{1-{\theta }_2}{{\theta }_2}{\varphi }_{\Delta }\left(-\xi ;x\right){\widehat{z}}^k\left(t_{m+1},\xi \right)-\\ i\frac{1}{{\theta }_2}{\Sigma }^T\left(t_m,x,y_{*}^{k-1}\left(t_m,x\right)\right){\varphi }_{\Delta }\left(-\xi ;x\right){\widehat{y}}^k\left(t_{m+1},\xi \right)\\ -i\frac{1-{\theta }_2}{{\theta }_2}{\Sigma }^T\left(t_m,x,y_{*}^{k-1}\left(t_m,x\right)\right){\varphi }_{\Delta }\left(-\xi ;x\right){\widehat{f}}^k\left(t_{m+1},\xi \right)\Delta t_{m+1},\end{array}$$

(41)

where

$${\widehat{\gamma }}^k\left(t_{m+1},\xi \right)={\varphi }_{\Delta }\left(-\xi ;x\right){\widehat{y}}^k\left(t_{m+1},\xi \right)+\left(1-{\theta }_1\right){\varphi }_{\Delta }\left(-\xi ;x\right){\widehat{f}}^k\left(t_{m+1},\xi \right)\Delta t_{m+1}.$$

(42)

Based on the above analysis, we give the following detailed Algorithm 2 to implement the numerical solution to Equation (1):

Algorithm 2. Local algorithm

Initialization:

1: Calculate ${\widehat{y}}^{P_{local}}\left(t_M,\xi \right),{\widehat{z}}^{P_{local}}\left(t_M,\xi \right)$ and ${\widehat{f}}^{P_{local}}\left(t_M,\xi \right)$ applying Equation (14) by FrFFT on the $\xi$-grid.

2: Calculate $y_{*}^{P_{local}}\left(t_M,x\right)$ and $z_{*}^{P_{local}}\left(t_M,x\right)$ by terminal conditions.

Main loop:

3: for $m=M-1$ to 1 do

4:  $y_{*}^0\left(t_m,x\right)\leftarrow y_{*}^{P_{local}}\left(t_{m+1},x\right)$

5:  $z_{*}^0\left(t_m,x\right)\leftarrow z_{*}^{P_{local}}\left(t_{m+1},x\right)$

6:  for $k=1$ to $P_{local}$ do

7:   Calculate $\widehat{z}\left(t_m,\xi \right)$ and $\widehat{\gamma }\left(t_{m+1},\xi \right)$ via Equations (41) and (42).

8:   Calculate $z\left(t_m,x\right)$ and $\gamma \left(t_{m+1},x\right)$ applying Equation (15) by FrFFT on the $x$-grid.

9:   Do Picard iteration $P$ times for Equation (36) to get $y^k\left(t_m,x\right)$ on the $x$-grid.

10:   Get $f^k\left(t_m,x\right)$ 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\left(t_m,x\right)$ and $z_{*}^k\left(t_m,x\right)$ at $x_m$.

13:   end if

14:   if $\max \left(\left|y_{*}^k\left(t_m,x\right)-y_{*}^{k-1}\left(t_m,x\right)\right|,\left|z_{*}^k\left(t_m,x\right)-z_{*}^{k-1}\left(t_m,x\right)\right|\right)\le {\mathsf{\epsilon }}_1$ then

15:    Break.

16:   end if

17:   if $m\ne 0$ then

18:    Calculate ${\widehat{y}}^k\left(t_m,\xi \right)$ and ${\widehat{f}}^k\left(t_m,\xi \right)$ applying Equation (14) by FrFFT on the $\xi$-grid.

19:   end if

20:  end for

21: end for

The total complexity of the local algorithm is $O\left(M{P}_{local}\left(DN+D{N}^{\iota }+D{N}^{\iota }log\left(N\right)+P{N}^{\iota }+1\right)\right)$.

Different: Each iteration uses the cubic spline to get the desired value (i.e.,$y_{*}^k\left(t_m,x\right)$).

4.3. 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,\cdots ,P_{global}$ as the current global iteration. Furthermore, we denote by $X_m^{\Delta ,k},{\widehat{y}}^k\left(t_m,\xi \right),{\widehat{z}}^k\left(t_m,\xi \right),{\widehat{f}}^k\left(t_m,\xi \right),{\widehat{\gamma }}^k\left(t_m,\xi \right),y^k\left(t_m,x\right),z^k\left(t_m,x\right),f^k\left(t_m,x\right),{\gamma }^k\left(t_m,x\right)$ and ${\gamma }^k\left(t_m,x\right)$ respectively the value of $X_m^{\Delta },\widehat{y}\left(t_m,\xi \right),\widehat{z}\left(t_m,\xi \right),\widehat{f}\left(t_m,\xi \right),\widehat{\gamma }\left(t_m,\xi \right),y\left(t_m,x\right),z\left(t_m,x\right),f\left(t_m,x\right)$ and $\gamma \left(t_m,x\right)$ in iteration $k^{th}$. Finally, denote by the values of $y_{*}^k\left(t_m,x\right)$ and $z_{*}^k\left(t_m,x\right)$ at $X_m$ on the $x$-grid using the cubic spline in the $k^{th}$ iteration. Then the Equation (18) will be rewritten by:

$$\begin{array}{c}{\varphi }_{\Delta }^{\mathrm{k}}\left(-\xi ;x\right)=exp(-\Delta t_m(\frac{1}{2}{\xi }^T\Sigma \left(t_{m+1},x,y_{*}^k\left(t_{m+1},x\right)\right){\Sigma }^T\left(t_{m+1},x,y_{*}^k\left(t_{m+1},x\right)\right)\xi +\\ i{\xi }^T\mu \left(t_{m+1},x,y_{*}^k\left(t_{m+1},x\right),z_{*}^k\left(t_{m+1},x\right)\right)\left)\right).\end{array}$$

(43)

However, we do not have functions with $y^k\left(t_m,x\right)$ and $z^k\left(t_m,x\right)$ at iteration $k$ and time point $t_m$, so we also apply the previous iteration:

$$\begin{array}{c}{\varphi }_{\Delta }^{\mathrm{k}}\left(-\xi ;x\right)\approx exp(-\Delta t_m(\frac{1}{2}{\xi }^T\Sigma \left(t_{m+1},x,y_{*}^{k-1}\left(t_{m+1},x\right)\right){\Sigma }^T\left(t_{m+1},x,y_{*}^{k-1}\left(t_{m+1},x\right)\right)\xi +\\ i{\xi }^T\mu \left(t_{m+1},x,y_{*}^{k-1}\left(t_{m+1},x\right),z_{*}^{k-1}\left(t_{m+1},x\right)\right)\left)\right).\end{array}$$

(44)

Using same logic as in Section 4.2, we have:

$$\begin{array}{c}{y}^{P_{global}}\left(t_M,x\right)=g\left(x\right),\hspace{1em}\hspace{1em}{z}^{P_{global}}\left(t_M,x\right)=\hfill \\ {\Sigma }^T\left(t_M,x,y_{*}^k\left(t_M,x\right)\right)D_{x}g\left(x\right),\hfill \end{array}$$

(45)

$$y^k\left(t_m,x\right)={\gamma }^k\left(t_{m+1},X_{m+1}\right)+\frac{1}{{\theta }_1}{f}^k\left(t_m,x\right)\Delta t_{m+1},$$

(46)

$$\begin{array}{c}{z}^k\left(t_m,x\right)=-\frac{1-{\theta }_2}{{\theta }_2}{𝒫}_{\Delta }{z}^k\left(t_{m+1},X_{m+1}\right)+\frac{1}{{\theta }_2}{\Sigma }^T\left(t_{m+1},x,y_{*}^{k-1}\left(t_m,x\right)\right){𝒫}_{\Delta }{y}^k\left(t_{m+1},X_{m+1}\right)\\ +\frac{1-{\theta }_2}{{\theta }_2}{\Sigma }^T\left(t_m,x,y_{*}^{k-1}\left(t_m,x\right)\right){𝒫}_{\Delta }{f}^k\left(t_{m+1},X_{m+1}\right)\Delta t_{m+1},\end{array}$$

(47)

where

$${\gamma }^k\left(t_{m+1},X_{m+1}\right)={𝒫}_{\Delta }{y}^k\left(t_{m+1},X_{m+1}\right)+\left(1-{\theta }_1\right){𝒫}_{\Delta }{f}^k\left(t_{m+1},X_{m+1}\right)\Delta t_{m+1},$$

(48)

and

$$\begin{array}{cc}\hfill {\widehat{y}}^{P_{global}}\left(t_M,\xi \right)& ={\mathbb{F}}_{\iota }\left(g\left(x\right)\right)\left(\xi \right),\hspace{1em}\hspace{1em}{\widehat{z}}^{P_{global}}\left(t_M,\xi \right)\hfill \\ & ={\Sigma }^T\left(t_M,x,y_{*}^k\left(t_M,x\right)\right){\mathbb{F}}_{\iota }\left(D_{x}g\left(x\right)\right)\left(\xi \right),\hfill \end{array}$$

(49)

$${\widehat{y}}^k\left(t_m,\xi \right)={\widehat{\gamma }}^k\left(t_{m+1},\xi \right)+{\theta }_1{\widehat{f}}^k\left(t_m,\xi \right)\Delta t_{m+1},$$

(50)

$$\begin{array}{c}{\widehat{z}}^k\left(t_m,\xi \right)=-\frac{1-{\theta }_2}{{\theta }_2}{\varphi }_{\Delta }\left(-\xi ;x\right){\widehat{z}}^k\left(t_{m+1},\xi \right)-\\ i\frac{1}{{\theta }_2}{\Sigma }^T\left(t_m,x,y_{*}^{k-1}\left(t_m,x\right)\right){\varphi }_{\Delta }\left(-\xi ;x\right){\widehat{y}}^k\left(t_{m+1},\xi \right)\\ -i\frac{1-{\theta }_2}{{\theta }_2}{\Sigma }^T\left(t_m,x,y_{*}^{k-1}\left(t_m,x\right)\right){\varphi }_{\Delta }\left(-\xi ;x\right){\widehat{f}}^k\left(t_{m+1},\xi \right)\Delta t_{m+1},\end{array}$$

(51)

where

$${\widehat{\gamma }}^k\left(t_{m+1},\xi \right)={\varphi }_{\Delta }\left(-\xi ;x\right){\widehat{y}}^k\left(t_{m+1},\xi \right)+\left(1-{\theta }_1\right){\varphi }_{\Delta }\left(-\xi ;x\right){\widehat{f}}^k\left(t_{m+1},\xi \right)\Delta t_{m+1}.$$

(52)

Then the process of Algorithm 3 is as follows:

Algorithm 3. Global algorithm

Initialization:

1: Calculate ${\widehat{y}}^{P_{global}}\left(t_M,\xi \right),{\widehat{z}}^{P_{global}}\left(t_M,\xi \right)$ and ${\widehat{f}}^{P_{global}}\left(t_M,\xi \right)$ applying Equation (14) by FrFFT on the $\xi$-grid.

Main loop:

2: for $m=M-1$ to 1 do

3: $y_{*}^0\left(t_m,x\right)\leftarrow 0$

4: $z_{*}^0\left(t_m,x\right)\leftarrow 0$

5: end for

6: for $k=1$ to $P_{global}$ do

7:  for $m=M-1$ to 1 do

8:   Calculate $\widehat{z}\left(t_m,\xi \right)$ and $\widehat{\gamma }\left(t_{m+1},\xi \right)$ via Equations (51) and (52).

9:   Calculate $z\left(t_m,x\right)$ and $\gamma \left(t_{m+1},x\right)$ applying Equation (15) by FrFFT on the $x$-grid.

10:   Do Picard iteration $P$ times for Equation (36) to get $y^k\left(t_m,x\right)$ on the $x$-grid.

11:   Get $f^k\left(t_m,x\right)$ 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\left(t_m,x\right)$ and $z_{*}^k\left(t_m,x\right)$ at $x_m$.

14:   end if

15:   if $m\ne 0$ then

16:    Calculate ${\widehat{y}}^k\left(t_m,\xi \right)$ and ${\widehat{f}}^k\left(t_m,\xi \right)$ applying Equation (14) by FrFFT on the $\xi$-grid.

17:   end if

18:  end for

19:  if $\underset{1\le \mathrm{m}\le \mathrm{M}}{\max }\left(\left|y_{*}^k\left(t_m,x\right)-y_{*}^{k-1}\left(t_m,x\right)\right|,\left|z_{*}^k\left(t_m,x\right)-z_{*}^{k-1}\left(t_m,x\right)\right|\right)\le {\mathsf{\epsilon }}_2$ then

20:   Break.

21:  end if

21: end for

The overall complexity of the global algorithm is $O\left(M{P}_{global}\left(DN+D{N}^{\iota }+D{N}^{\iota }log\left(N\right)+P{N}^{\iota }+1\right)\right)$.

Different: Like the local algorithm, each iteration uses the cubic spline to get the desired value (i.e., $y_{*}^k\left(t_m,x\right)$).

5. 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. ${\epsilon }_1$ and ${\epsilon }_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 $\left|Y\left(0\right)-Y^{\Delta }\left(0\right)\right|$ for $Y$ and $\left|Z\left(0\right)-Z^{\Delta }\left(0\right)\right|$ for $Z$.

Due to the Fourier method used in this paper, the integration interval $\left[a,b\right]$ is referenced by [10,31] and the integration interval $\left[\tilde{a},\tilde{b}\right]$ is referenced by [7].

For each method, we tested the following $\theta$ schemes:

$$\mathrm{A}:{\theta }_1=\frac{1}{2},{\theta }_2=\frac{1}{2}$$

$$\mathrm{B}:{\theta }_1=1,{\theta }_2=1$$

$$\mathrm{C}:{\theta }_1=0,{\theta }_2=1$$

$$\mathrm{D}:{\theta }_1=1,{\theta }_2=\frac{1}{2}$$

5.1. Example 1: A Decoupled Example

We consider the following decoupled FBSDE:

$$\{\begin{array}{c}d{X}_t=d{W}_t,\\ d{Y}_t=-f\left(t,X_t,Y_t,Z_t\right)dt+Z_{t}d{W}_t,\end{array}$$

(53)

where

$$\begin{array}{c}f\left(t,X_t,Y_t,Z_t\right)=Y_t{Z}_t-Z_t+\frac{9}{8}{Y}_t-\frac{1}{2}sin\left(\frac{1}{2}{X}_t+t\right)cos\left(\frac{1}{2}{X}_t+t\right)-sin\left(\frac{1}{2}{X}_t+t\right)-\\ \frac{1}{2}cos{\left(\frac{1}{2}{X}_t+t\right)}_{.}\end{array}$$

The solution of Equation (53) is given by:

$$\left(y\left(t,x\right),z\left(t,x\right)\right)=\left(sin\left(\frac{1}{2}x+t\right),\frac{1}{2}cos\left(\frac{1}{2}x+t\right)\right).$$

In this experiment, we use $T=1$ and $X_0=0$ and set $a=-10,b=10,\tilde{a}=-45,$ and $\tilde{b}=45$. It is known that the exact solution is $\left(y\left(0,x_0\right),z\left(0,x_0\right)\right)=\left(0,\frac{1}{2}\right)$. Since $\mu$ and $\Sigma$ 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\left(0,x_0\right)$ and $z\left(0,x_0\right)$ of different values of $N=2^d$. The error is listed in

Table 1. Error of results ($Y_0$, $Z_0$) for example 1 by scheme A.

$\mathit{M}$	Error	$\mathit{d}=9$	Time (s)	$\mathit{d}=10$	Time (s)	$\mathit{d}=11$	Time (s)
4	$Y_0$	$2.7\times {10}^{-3}$	0.010592	$1.5\times {10}^{-3}$	0.02287	$2.6\times {10}^{-3}$	0.024719
	$Z_0$	3$.4\times {10}^{-3}$		$8.8\times {10}^{-3}$		$1.19\times {10}^{-2}$
8	$Y_0$	5$.7\times {10}^{-3}$	0.019609	$9\times {10}^{-4}$	0.034164	$4\times {10}^{-4}$	0.041977
	$Z_0$	$1.82\times {10}^{-2}$		$2.1\times {10}^{-3}$		$2\times {10}^{-3}$
16	$Y_0$	$6.5\times {10}^{-3}$	0.03041	$1.5\times {10}^{-3}$	0.062138	$2\times {10}^{-4}$	0.065344
	$Z_0$	$2.21\times {10}^{-2}$		$5\times {10}^{-3}$		$6.32\times {10}^{-4}$
32	$Y_0$	$6.7\times {10}^{-3}$	0.060005	$1.7\times {10}^{-3}$	0.100117	$4\times {10}^{-4}$	0.128925
	$Z_0$	$2.31\times {10}^{-2}$		$5.8\times {10}^{-3}$		$1.3\times {10}^{-3}$
64	$Y_0$	$6.7\times {10}^{-3}$	0.099635	$1.7\times {10}^{-3}$	0.218806	$4\times {10}^{-4}$	0.266788
	$Z_0$	$2.34\times {10}^{-2}$		$5.9\times {10}^{-3}$		$1.5\times {10}^{-3}$
128	$Y_0$	$6.8\times {10}^{-3}$	0.243116	$1.8\times {10}^{-3}$	0.447437	$4\times {10}^{-4}$	0.584943
	$Z_0$	$2.34\times {10}^{-2}$		$6\times {10}^{-3}$		$1.5\times {10}^{-3}$
256	$Y_0$	$6.8\times {10}^{-3}$	0.530708	$1.8\times {10}^{-3}$	1.006999	$4\times {10}^{-4}$	1.499762
	$Z_0$	$2.35\times {10}^{-2}$		$6\times {10}^{-3}$		$1.5\times {10}^{-3}$
512	$Y_0$	$6.8\times {10}^{-3}$	1.455568	$1.7\times {10}^{-3}$	2.657796	$4\times {10}^{-4}$	3.964449
	$Z_0$	$2.36\times {10}^{-2}$		$6\times {10}^{-3}$		$1.5\times {10}^{-3}$

, where the approximations are obtained by the scheme A. Furthermore, the error in the approximation of $y\left(0,x_0\right)$ and $z\left(0,x_0\right)$ 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.

5.2. Example 2: Forward Process Depending on X_t

We consider the following coupled FBSDE:

$$\{\begin{array}{c}d{X}_t=\frac{1}{1+2exp\left(t+X_t\right)}dt+\frac{exp\left(t+X_t\right)}{1+exp\left(t+X_t\right)}d{W}_t,\\ d{Y}_t=-f\left(t,X_t,Y_t,Z_t\right)dt+Z_{t}d{W}_t,\end{array}$$

(54)

where

$$f\left(t,X_t,Y_t,Z_t\right)=-\frac{2Y_t}{1+exp\left(t+X_t\right)}-\frac{1}{2}{\left(\frac{Y_t{Z}_t}{1+exp\left(t+X_t\right)}-Y_t^2{Z}_t\right)}_{.}$$

The solution of Equation (54) is given by:

$$\left(y\left(t,x\right),z\left(t,x\right)\right)={\left(\frac{exp\left(t+x\right)}{1+exp\left(t+x\right)},\frac{{\left(exp\left(t+x\right)\right)}^2}{{\left(1+exp\left(t+x\right)\right)}^3}\right)}_{.}$$

For the experiment, we use $T=0.1$ and $X_0=0$ and set $a=-8,b=8,\tilde{a}=-20,$ and $\tilde{b}=20$. It is known that the exact solution is $\left(y\left(0,x_0\right),z\left(0,x_0\right)\right)=\left(\frac{1}{2},\frac{1}{8}\right)$.

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. Error of results ($Y_0$, $Z_0$) for example 2 by scheme C.

$\mathit{M}$	Error	$\mathit{d}=9$	Time (s)	$\mathit{d}=10$	Time (s)	$\mathit{d}=11$	Time (s)
4	$Y_0$	$8.8\times {10}^{-3}$	0.005422	$8.8\times {10}^{-3}$	0.008908	$8.8\times {10}^{-3}$	0.019683
	$Z_0$	$8.55\times {10}^{-4}$		$1.2\times {10}^{-3}$		$1.3\times {10}^{-3}$
8	$Y_0$	$4.8\times {10}^{-3}$	0.009573	$4.7\times {10}^{-3}$	0.016703	$4.7\times {10}^{-3}$	0.032749
	$Z_0$	$3.54\times {10}^{-4}$		$6.01\times {10}^{-4}$		$6.52\times {10}^{-4}$
16	$Y_0$	$2.5\times {10}^{-3}$	0.025286	$2.4\times {10}^{-3}$	0.028508	$2.4\times {10}^{-3}$	0.052708
	$Z_0$	$8.73\times {10}^{-4}$		$1.3\times {10}^{-3}$		$1.4\times {10}^{-3}$
32	$Y_0$	$1.3\times {10}^{-3}$	0.04016	$1.3\times {10}^{-3}$	0.05021	$1.3\times {10}^{-3}$	0.097333
	$Z_0$	$3.1\times {10}^{-3}$		$4.4\times {10}^{-3}$		$4.8\times {10}^{-3}$
64	$Y_0$	$8.67\times {10}^{-4}$	0.077589	$9.42\times {10}^{-4}$	0.099154	$9.69\times {10}^{-4}$	0.207257
	$Z_0$	$5.7\times {10}^{-3}$		$7.4\times {10}^{-3}$		$7.8\times {10}^{-3}$
128	$Y_0$	$7.26\times {10}^{-4}$	0.159869	$9.34\times {10}^{-4}$	0.20815	$9.94\times {10}^{-4}$	0.432267
	$Z_0$	$7.2\times {10}^{-3}$		$8.4\times {10}^{-3}$		$8.6\times {10}^{-3}$
256	$Y_0$	$9.33\times {10}^{-4}$	0.325549	$1.2\times {10}^{-3}$	0.463633	$1.3\times {10}^{-3}$	0.982044
	$Z_0$	$6.9\times {10}^{-3}$		$6.6\times {10}^{-3}$		$6.3\times {10}^{-3}$
512	$Y_0$	$8.26\times {10}^{-4}$	0.759207	$1\times {10}^{-3}$	1.156449	$1.1\times {10}^{-3}$	2.491836
	$Z_0$	$4.5\times {10}^{-3}$		$2.3\times {10}^{-3}$		$1.6\times {10}^{-3}$

shows the error for the approximations of $y\left(0,x_0\right)$ and $z\left(0,x_0\right)$ 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 Figure 2 and Figure 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.

5.3. Example 3: Diffusion Depending on Y_t

We consider the following coupled FBSDE [5]:

$$\{\begin{array}{c}d{X}_t=\frac{X_t\left(1+X_t^2\right)}{{\left(2+X_t^2\right)}^3}dt+\frac{1+X_t^2}{2+X_t}\sqrt{\frac{1+2Y_t^2}{1+Y_t^2+exp\left(-\frac{2X_t^2}{t+1}\right)}}d{W}_t,\\ d{Y}_t=-f\left(t,X_t,Y_t,Z_t\right)dt+Z_{t}d{W}_t,\end{array}$$

(55)

where

$$f\left(t,X_t,Y_t,Z_t\right)=a\left(t,X_t\right)-b\left(t,X_t,Y_t\right)Z_t.$$

$$a\left(t,x\right)=\frac{1}{t+1}exp\left(-\frac{x^2}{t+1}\right){\left[\frac{4x^2\left(1+x^2\right)}{{\left(2+x^2\right)}^3}+{\left(\frac{1+x^2}{2+x^2}\right)}^2\left(1-\frac{2x^2}{t+1}\right)-\frac{x^2}{t+1}\right]}_{}$$

and

$$b\left(t,x,y\right)=\frac{x}{{\left(2+x^2\right)}^2}\sqrt{\frac{1+y^2+exp\left(-\frac{2x^2}{t+1}\right)}{1+2y^2}}.$$

The solution of Equation (55) is given by:

$$\left(y\left(t,x\right),z\left(t,x\right)\right)={\left(exp\left(-\frac{x^2}{t+1}\right),-\frac{2x\left(1+x^2\right)}{\left(1+t\right)\left(2+x^2\right)}exp\left(-\frac{x^2}{t+1}\right)\right)}_{.}$$

For the experiment, we use $T=1$ and $X_0=0$ and set $a=-5,b=5,\tilde{a}=-30,$ and $\tilde{b}=30$. The correct solution is $\left(y\left(0,x_0\right),z\left(0,x_0\right)\right)=\left(1,0\right)$.

The results of the explicit algorithm, the local algorithm, and the global algorithm can be found separately in Figure 4, Figure 5 and Figure 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.

Table 3. CPU times of scheme A for example 3(s).

$\left(\mathit{d}=9\right)\mathit{M}$	4	8	16	32	64	128	256	512
Explicit	0.012533	0.025716	0.035425	0.066304	0.155958	0.283224	0.613946	1.558809
Local	0.043708	0.070042	0.104299	0.227083	0.442085	0.909928	2.054779	4.157252
Global	0.052723	0.090697	0.182892	0.326598	0.623499	1.280217	2.141616	4.563165

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.

5.4. Example 4: Drift Depending on Z_t

We consider the following coupled FBSDE [17]:

$$\{\begin{array}{c}d{X}_t=cos\left(t+X_t\right)\left(Y_t+Z_t\right)dt+\sqrt{2}\left(Y_{t}sin\left(t+X_t\right)+1\right)d{W}_t,\\ d{Y}_t=-f\left(t,X_t,Y_t,Z_t\right)dt+Z_{t}d{W}_t,\end{array}$$

(56)

where

$$f\left(t,X_t,Y_t,Z_t\right)=-cos\left(t+X_t\right)-co{s}^2\left(t+X_t\right)Z_t+\left(3si{n}^2\left(t+X_t\right)+si{n}^4\left(t+X_t\right)\right)Y_t.$$

The solution of Equation (56) is given by:

$$\left(y\left(t,x\right),z\left(t,x\right)\right)={\left(sin\left(t+x\right),\sqrt{2}cos\left(t+x\right)\left(si{n}^2\left(t+x\right)+1\right)\right)}_{.}$$

For the experiment, we use $T=0.1$ and $X_0=0$ and set $a=-8.5614,b=9.3271,\tilde{a}=-20,$ and $\tilde{b}=20$. It is known that the exact solution is $\left(y\left(0,x_0\right),z\left(0,x_0\right)\right)=\left(0,\sqrt{2}\right)$.

The results for the three algorithms can be found in Figure 7, Figure 8 and Figure 9. For $y\left(0,x_0\right)$, we observe that the three algorithms obviously converge. For $z\left(0,x_0\right)$, 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. CPU times of scheme A for example 4 (s).

$\left(\mathit{d}=9\right)\mathit{M}$	4	8	16	32	64	128	256	512
Explicit	0.009843	0.021973	0.035312	0.061192	0.125936	0.271638	0.629419	1.451987
Local	0.065195	0.105617	0.202525	0.326405	0.538356	1.148581	2.342291	6.382512
Global	0.080436	0.170601	0.310428	0.518835	0.794979	1.540419	3.223398	8.397285

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.

5.5. Example 5: Drift Depending on Z_t

We consider the following coupled FBSDE [17]:

$$\{\begin{array}{c}d{X}_t=cos\left(t+X_t\right)\left(Y_t+Z_t\right)dt+\sqrt{2}{Y}_{t}sin\left(t+X_t\right)d{W}_t,\\ d{Y}_t=-f\left(t,X_t,Y_t,Z_t\right)dt+Z_{t}d{W}_t,\end{array}$$

(57)

where

$$f\left(t,X_t,Y_t,Z_t\right)=-cos\left(t+X_t\right)-Y_t-Z_t+si{n}^2\left(t+X_t\right){\left(Y_t+Z_t+Y_t^3\right)}_{.}$$

The solution of Equation (57) is given by:

$$\left(y\left(t,x\right),z\left(t,x\right)\right)={\left(sin\left(t+x\right),\sqrt{2}cos\left(t+x\right)si{n}^2\left(t+x\right)\right)}_{.}$$

In the experiment, we use $T=0.1$ and $X_0=0$ and set $a=-4.2307,b=4.7135,\tilde{a}=-15$, and $\tilde{b}=15$. It is known that the exact solution is $\left(y\left(0,x_0\right),z\left(0,x_0\right)\right)=\left(0,0\right)$.

The results for the explicit algorithm, the local algorithm, and the global algorithm can be found separately in Figure 10, Figure 11 and Figure 12. For $y\left(0,x_0\right)$, we observe that the effects of the three algorithms do not differ significantly. Although the error from the figures increases when $M\ge 200$, they converge at a certain number. For $z\left(0,x_0\right)$, 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. CPU times of scheme A for the example 5 (s).

$\left(\mathit{d}=9\right)\mathit{M}$	4	8	16	32	64	128	256	512
Explicit	0.009185	0.016861	0.028855	0.051736	0.097597	0.209269	0.503181	1.276852
Local	0.037164	0.054169	0.090534	0.175027	0.339829	0.720882	1.554783	4.157505
Global	0.049146	0.079958	0.160758	0.280679	0.526248	1.199261	2.558428	5.285499

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.

5.6. Example 6: Diffusion Depending on Z_t

We consider the following coupled FBSDE:

$$\{\begin{array}{c}d{X}_t=\frac{1}{2}sin\left(t+X_t\right)cos\left(t+X_t\right)\left(Y_t^2+Z_t\right)dt+\frac{1}{2}cos\left(t+X_t\right)\left(Y_{t}sin\left(t+X_t\right)+Z_t+1\right)d{W}_t,\\ d{Y}_t=-f\left(t,X_t,Y_t,Z_t\right)dt+Z_{t}d{W}_t,\end{array}$$

(58)

where

$$f\left(t,X_t,Y_t,Z_t\right)=Y_t{Z}_t-cos{\left(t+X_t\right)}_{.}$$

This FBSDE satisfies assumptions (1)~(4). The solution of Equation (58) is given by:

$$\left(y\left(t,x\right),z\left(t,x\right)\right)={\left(sin\left(t+x\right),co{s}^2\left(t+x\right)\right)}_{.}$$

For the experiment, we use $T=1$ and $X_0=0$ and set $a=-2\pi ,b=2\pi ,\tilde{a}=-2\pi ,$ and $\tilde{b}=2\pi$. It is known that the exact solution is $\left(y\left(0,x_0\right),z\left(0,x_0\right)\right)=\left(0,1\right)$.

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 Figure 13, Figure 14 and Figure 15. For $y\left(0,x_0\right)$ 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\left(0,x_0\right)$, 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\left(0,x_0\right)$ using the local and global methods have a divergent tendency. Their precision needs to be improved.

Table 6. CPU times of scheme B for example 6 (s).

$\left(\mathit{d}=9\right)\mathit{M}$	4	8	16	32	64	128	256	512
Explicit	0.007739	0.012715	0.024707	0.039612	0.077765	0.187745	0.405907	1.122424
Local	0.066662	0.106036	0.210299	0.403093	0.815283	1.855567	4.125572	11.050587
Global	0.089035	0.201116	0.334744	0.673482	1.176203	2.425878	4.909738	10.333279

shows the CPU times for the explicit, local, and global algorithms. The speed of all algorithms is still relatively great.

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