A Fourier Interpolation Method for Numerical Solution of FBSDEs: Global Convergence, Stability, and Higher Order Discretizations ^†

Abstract

The convolution method for the numerical solution of forward-backward stochastic differential equations (FBSDEs) was originally formulated using Euler time discretizations and a uniform space grid. In this paper, we utilize a tree-like spatial discretization that approximates the BSDE on the tree, so that no spatial interpolation procedure is necessary. In addition to suppressing extrapolation error, leading to a globally convergent numerical solution for the FBSDE, we provide explicit convergence rates. On this alternative grid the conditional expectations involved in the time discretization of the BSDE are computed using Fourier analysis and the fast Fourier transform (FFT) algorithm. The method is then extended to higher-order time discretizations of FBSDEs. Numerical results demonstrating convergence are presented using a commodity price model, incorporating seasonality, and forward prices.

1. Introduction

A variety of numerical methods for backward stochastic differential equations (BSDEs) and forward-backward stochastic differential equations (FBSDEs) have been developed recently. Different applications call for innovative techniques for the efficient resolution of these systems. In finance and economics BSDEs are used for option pricing and hedging El Karoui et al. (1997), reflected BSDEs are used for modeling American options El Karoui et al. (1997), and quadratic BSDEs play an important role in continuous-time recursive utility Duffie and Epstein (1992) and other utility maximization problems.

Following the establishment of the well-posedness of BSDEs by Pardoux and Peng (1992) the first numerical procedures to emerge were partial differential equation (PDE) based methods such as the finite difference approach of Douglas et al. (1996). The PDE method is mainly devoted to coupled problems as is the spectral method of Ma et al. (2008). More recent numerical methods based on machine learning (Beck et al. 2019; Han and Long 2020; E et al. 2017, 2019, 2022) have been applied to the numerical solution of BSDEs which, given the connections between the theory of PDEs and their representations as BSDEs, in turn provides solutions to high dimensional PDEs. Spatial discretization methods were initiated by Chevance (1997) with a quantization approach to conditional expectations. However, it is only since Zhang (2004) and Bouchard and Touzi (2004) that a sound time discretization of (decoupled) FBSDEs is available. The quantization approach was then used in the multidimensional framework of Bally and Pages (2003); Bally et al. (2005) and the coupled FBSDE case of Delarue and Menozzi (2006). The theoretical basis for a multinomial approach for BSDEs was introduced by Briand et al. (2001) and Ma et al. (2002) and followed in practice by Peng and Xu (2011). Monte Carlo methods are the most prolific approach for numerical solutions of (F)BSDEs. They include the backward scheme of Zhang (2004), the Malliavin approach of Bouchard and Touzi (2004) and Crisan et al. (2010), the least-square regression approach of Gobet et al. (2005) and the iterative schemes Bender and Denk (2007) and Bender and Zhang (2008).

Additional approaches to numerical solution of BSDEs include the cubature method Crisan and Manolarakis (2012, 2014), Fourier-cosine expansions Huijskens et al. (2016); Ruijter and Oosterlee (2015, 2016), and the convolution method Hyndman and Oyono Ngou (2017). In introducing the convolution method, Hyndman and Oyono Ngou (2017) developed a local discretization error which includes an extrapolation error. The extrapolation error component is exclusively produced by the fast Fourier transform (FFT) algorithm and the underlying trigonometric interpolation used to compute conditional expectations. To improve the performance of the convolution method it is desirable to eliminate the extrapolation error, and improve the error bound, with an alternative implementation of the FFT algorithm.

In this paper, we propose an alternative space grid for the convolution method, instead of the rectangular grid in Hyndman and Oyono Ngou (2017), which eliminates the extrapolation error, leads to a globally convergent numerical solution for the (F)BSDE, and provides explicit convergence rates. We also apply the numerical method to the Runge–Kutta schemes for FBSDEs proposed by Chassagneux and Crisan (2014). The tree-like nature of the alternative grid avoids extrapolations and leads to a global error bound for the BSDE approximate solutions. Further, the implementation of the convolution method originally presented in Hyndman and Oyono Ngou (2017) is simplified by using an alternative parametric transformation to enforce the necessary periodic boundary conditions.

The paper is organized as follows. Section 2 reviews the explicit Euler time discretization of BSDEs, recalls the convolution method of Hyndman and Oyono Ngou (2017) for computing the necessary conditional expectations, gives a description of an alternative spatial discretization, and provides a generic implementation of the convolution method on this grid using the discrete Fourier transform. The section ends with a global error analysis. Section 3 extends the Fourier interpolation method to higher order time discretizations of FBSDEs and includes the related global error analysis. Finally, Section 4 presents a numerical implementation, in the context of a commodity price model, that illustrates the theoretical results. Section 5 concludes.

2. The Fourier Interpolation Method

In this section, we introduce an alternative grid that removes extrapolation error of Hyndman and Oyono Ngou (2017), after a quick presentation of the Euler scheme for FBSDEs. Section 2.4 presents a global error analysis of the Fourier interpolation method on this alternative grid and under the Euler scheme.

2.1. Time Discretization

Let $\left(\mathsf{\Omega },\mathbf{P},\mathcal{F},{\left\{{\mathcal{F}}_t\right\}}_{t\in [0,T]}\right)$ be a complete filtered probability space generated by a $d-$dimensional Wiener process W. We seek a numerical solution to the FBSDE

$$\left\{\begin{array}{c}d{X}_t=a(t,X_t)dt+\sigma (t,X_t)d{W}_t\hfill \\ -d{Y}_t=f(t,X_t,Y_t,Z_t)dt-Z_t^{*}d{W}_t\hfill \\ X_0=x_0,Y_T=\xi \hfill \end{array}\right.$$

(1)

The forward drift $a:[0,T]\times {\mathbb{R}}^d\to {\mathbb{R}}^d$, the forward volatility $\sigma :[0,T]\times {\mathbb{R}}^d\to {\mathbb{R}}^{d\times d}$, the driver $f:[0,T]\times {\mathbb{R}}^d\times \mathbb{R}\times {\mathbb{R}}^d\to \mathbb{R}$ are deterministic functions. The initial condition is $x_0\in {\mathbb{R}}^d$ and the terminal condition takes the Markovian form $\xi =g\left(X_T\right)$ where $g:{\mathbb{R}}^d\to \mathbb{R}$. The FBSDE coefficients satisfy Assumption 1 so that existence and uniqueness of the FBSDE solution $(X,Y,Z)$ is assured.

Assumption 1.

There exist positive constants $K_1$, $K_2$ $K_3$, and $K_4$ such that the coefficients of the FBSDE (1) satisfy

$$\begin{array}{cc}\hfill \left|a(t,x_1)-a(t,x_2)\right|+{\parallel \sigma (t,x_1)-\sigma (t,x_2)\parallel }_2& \le K_1\left|x_1-x_2\right|\hfill \end{array}$$

(2)

$$\begin{array}{cc}\hfill \left|a(t,x)\right|+{\parallel \sigma (t,x)\parallel }_2& \le K_2\hfill \end{array}$$

(3)

$$\begin{array}{cc}\hfill \left|f(t,x_1,y,z)-f(t,x_2,y,z)\right|& \le K_1\left|x_1-x_2\right|\hfill \end{array}$$

(4)

$$\begin{array}{cc}\hfill \left|f(t,x,y_1,z_1)-f(t,x,y_2,z_2)\right|& \le K_1\left(\left|y_1-y_2\right|+\left|z_1-z_2\right|\right)\hfill \end{array}$$

(5)

$$\begin{array}{cc}\hfill \left|f(t,x,y,z)\right|& \le K_3(1+\left|x\right|+\left|y\right|+\left|z\right|)\hfill \end{array}$$

(6)

for any $t\in [0,T]$, $x,x_1,x_2\in {\mathbb{R}}^d$, $y,y_1,y_2\in \mathbb{R}$, $z,z_1,z_2\in {\mathbb{R}}^d$.

Moreover ${\sigma }^2:=\sigma {\sigma }^{*}$ is (uniformly) invertible, continuous and bounded

$$\begin{array}{c}\hfill \parallel {\left({\sigma }^2(t,x)\right)}^{-1}{\parallel }_2\le K_4\end{array}$$

(7)

for any $t\in [0,T]$, $x\in {\mathbb{R}}^d$.

In addition, the terminal value is square integrable

$${\parallel \xi \parallel }_{L^2}^2:={\mathbf{E}}_{}^{}\left[{\left|g\left(X_T\right)\right|}^2\right]<\infty .$$

(8)

Remark 1.

Assumption 1 makes no explicit assumption on g except square integrability. More restrictive assumptions on g, namely that g is twice continuously differentiable, shall be specified in the main results of Section 2 and Section 3 which provide explicit rates of convergence for the Fourier interpolation algorithms. However, similar to Crisan and Manolarakis (2012), see also Turkedjiev (2015), it should be possible to consider g non-differentiable using a mollification argument on the terminal condition.1 However, we shall not follow this approach in this paper so as to keep the focus on the main contributions which are the overall approach, implementation of the convolution method on the tree-like grid, developing the Runge–Kutta discretization schemes, and explicit convergence rates.

The time discretization of the FBSDE (1) on the time partition $\pi =\{t_0=0<t_1<\dots <t_n=T\}$ consists of the explicit Euler scheme given by

$$\left\{\begin{array}{c}{X}_0^{\pi }=x_0\hfill \\ X_{t_{i+1}}^{\pi }=X_{t_i}^{\pi }+a(t_i,X_{t_i}^{\pi }){\Delta }_i+\sigma (t_i,X_{t_i}^{\pi })\Delta W_i\hfill \\ Z_{t_n}^{\pi }=0,Y_{t_n}^{\pi }={\xi }^{\pi }\hfill \\ Z_{t_i}^{\pi }=\frac{1}{{\Delta }_i}{\mathbf{E}}_{}^{}\left[Y_{t_{i+1}}^{\pi }\Delta W_i|{\mathcal{F}}_{t_i}\right]\hfill \\ Y_{t_i}^{\pi }={\mathbf{E}}_{}^{}\left[Y_{t_{i+1}}^{\pi }|{\mathcal{F}}_{t_i}\right]+f(t_i,X_{t_i}^{\pi },{\mathbf{E}}_{}^{}\left[Y_{t_{i+1}}^{\pi }|{\mathcal{F}}_{t_i}\right],Z_{t_i}^{\pi }){\Delta }_i\hfill \end{array}\right.$$

(9)

with ${\Delta }_i=t_{i+1}-t_i$ and $\Delta W_i=W_{t_{i+1}}-W_{t_i}$. Under an additional Lipschitz condition on the function g we have, from Zhang (2004) and Bouchard and Touzi (2004), that the quadratic discretization error

$$E_{\pi }^2:=\underset{0\le i<n}{max}{\mathbf{E}}_{}^{}\left[\underset{t\in [t_i,t_{i+1}]}{sup}{\left|Y_t-Y_{t_i}^{\pi }\right|}^2\right]+\sum _{i=0}^{n-1}{\mathbf{E}}_{}^{}\left[{\int }_{t_i}^{t_{i+1}}{\left|Z_s-Z_{t_i}^{\pi }\right|}^{2}ds\right]$$

(10)

is of first-order in time, i.e.,

$$E_{\pi }^2=\mathcal{O}\left(\left|\pi \right|\right)$$

(11)

where

$$\left|\pi \right|=\underset{i}{max}{\Delta }_i.$$

(12)

Following Hyndman and Oyono Ngou (2017) and Oyono Ngou (2014), the approximate solution $u_i$ and the approximate gradient ${\dot{u}}_i$ at time node $t_i$, $i=0,1,\dots ,n-1$, are given by

$$\begin{array}{cc}\hfill u_i\left(x\right)& ={\tilde{u}}_i\left(x\right)+{\Delta }_{i}f(t_i,x,{\tilde{u}}_i\left(x\right),\sigma (t_i,x){\dot{u}}_i\left(x\right))\hfill \end{array}$$

(13)

$$\begin{array}{cc}\hfill \sigma (t_i,x){\dot{u}}_i\left(x\right)& =\frac{1}{{\Delta }_i}{\mathbf{E}}_{}^{}\left[Y_{t_{i+1}}^{\pi }\sigma (t_i,X_{t_i}^{\pi })\Delta W_i|X_{t_i}^{\pi }=x\right]\hfill \end{array}$$

$$\begin{array}{cc}\hfill & =\frac{1}{{\Delta }_i}{\int }_{{\mathbb{R}}^d}^{}(y-{\Delta }_{i}a(t_i,x))u_{i+1}(x+y)h_i\left(y\right|x)dy\hfill \end{array}$$

(14)

where the intermediate solution ${\tilde{u}}_i$ is given by

$$\begin{array}{cc}\hfill {\tilde{u}}_i\left(x\right)& ={\mathbf{E}}_{}^{}\left[Y_{t_{i+1}}^{\pi }|X_{t_i}^{\pi }=x\right]\hfill \\ \hfill & ={\int }_{{\mathbb{R}}^d}^{}{u}_{i+1}(y+x)h_i\left(y\right|x)dy,\hfill \end{array}$$

(15)

$u_n=g$ and $h_i$ is a Gaussian density

$$h_i\left(y\right|x)={\left(2\pi \right)}^{-\frac{d}{2}}{\parallel {\Delta }_i{\sigma }^2(t_i,x)\parallel }_2^{-\frac{1}{2}}exp\left(-\frac{1}{2{\Delta }_i}{\mathbf{y}}^{*}{\left({\sigma }^2(t_i,x)\right)}^{-1}\mathbf{y}\right),$$

(16)

and where $\mathbf{y}=y-{\Delta }_{i}a(t_i,x)$ with characteristic function

$${\varphi }_i(\nu ,x)=exp\left(\mathbf{i}{\Delta }_i{\nu }^{*}a(t_i,x)-\frac{1}{2}{\Delta }_i{\nu }^{*}{\sigma }^2(t_i,x)\nu \right).$$

(17)

Let $\mathfrak{F}$ and ${\mathfrak{F}}^{-1}$ denote the Fourier transform operator and the inverse Fourier transform operator, respectively,

$$\begin{array}{cc}\hfill \mathfrak{F}\left[\theta \right]\left(\nu \right)& ={\int }_{{\mathbb{R}}^d}^{}{e}^{-\mathbf{i}{x}^{*}\nu }\theta \left(x\right)dx\hfill \end{array}$$

(18)

$$\begin{array}{cc}\hfill {\mathfrak{F}}^{-1}\left[\theta \right]\left(x\right)& ={\frac{1}{{\left(2\pi \right)}^d}}^{}{\int }_{{\mathbb{R}}^d}{e}^{\mathbf{i}{\nu }^{*}x}\theta \left(\nu \right)d\nu .\hfill \end{array}$$

(19)

The density function $h_i$ and the characteristic function ${\varphi }_i$ satisfy the relations

$$\begin{array}{cc}\hfill h_i\left(y\right|x)& =\frac{1}{{\left(2\pi \right)}^d}{\int }_{{\mathbb{R}}^d}{e}^{-\mathbf{i}{\nu }^{*}y}{\varphi }_i(\nu ,x)d\nu \hfill \\ \hfill y{h}_i\left(y\right|x)& =-\frac{1}{{\left(2\pi \right)}^d}{\int }_{{\mathbb{R}}^d}^{}{e}^{-\mathbf{i}{\nu }^{*}y}\mathbf{i}{\nabla }_{\nu }{\varphi }_i(\nu ,x)d\nu \hfill \end{array}$$

(20)

$$\begin{array}{cc}\hfill & ={\Delta }_{i}a(t_i,x)h_i\left(y\right|x)+\frac{{\Delta }_i{\sigma }^{*}(t_i,x)}{{\left(2\pi \right)}^d}{\int }_{{\mathbb{R}}^d}^{}{e}^{-\mathbf{i}{\nu }^{*}y}\mathbf{i}\nu {\varphi }_i(\nu ,x)d\nu \hfill \end{array}$$

(21)

Using the relationships between the characteristic function and the density function (20) and (21) leads to the representation

$$\begin{array}{cc}\hfill {\tilde{u}}_i\left(x\right)& ={\mathfrak{F}}^{-1}\left[\mathfrak{F}\left[u_{i+1}\right]\left(\nu \right){\varphi }_i(\nu ,x)\right]\left(x\right)\hfill \end{array}$$

(22)

$$\begin{array}{cc}\hfill {\dot{u}}_i\left(x\right)& ={\sigma }^{*}(t_i,x){\mathfrak{F}}^{-1}\left[\mathfrak{F}\left[u_{i+1}\right]\left(\nu \right)\mathbf{i}\nu {\varphi }_i(\nu ,x)\right]\left(x\right)\hfill \end{array}$$

(23)

for Equations (14) and (15) under integrability condition on the approximate solution $u_{i+1}$. In the sequel, we restrict the analysis to the one-dimensional case with $d=1$.

2.2. Space Discretization

The space discretization is performed on a tree-like grid using three parameters: the increment length $l>0$, the even number $N\in {\mathbb{N}}^{*}$ of space steps on the increment length, and the initial number $N_0\in \mathbb{N}$ of increment intervals. Hence, the space step is constant and uniform on the grid

$$\Delta x=\frac{l}{N}.$$

(24)

At the time node $t_i$, $i=0,1,\dots ,n$, the space domain is restricted on an interval of length $N_{i}l$ centred at $x_0$ and discretized uniformly with $N_{i}N$ space steps where

$$N_i=N_0+i$$

(25)

giving the nodes

$$x_{ik}=x_0-\frac{N_{i}l}{2}+k\Delta x,k=0,1,\dots ,N_{i}N.$$

(26)

In particular, the relation

$$x_{ik}=x_{i+1,k+\frac{N}{2}},k=0,1,\dots ,N_{i}N.$$

(27)

holds since the restricted interval at each time node is obtained by evenly increasing the previous one with an interval of length l. If $N_0=0$ then the space grid at mesh time $t_0$ is composed of the single point

$$x_{00}=x_0.$$

(28)

Figure 1 gives examples of alternative grids.

The convolution relations of Equations (22) and (23) call for a discretization of the Fourier space as well. At each mesh time $t_i$, $i=1,2,\dots ,n$, the Fourier space is restricted on an interval of length L centred at zero $\left(0\right)$ and discretized with $N_{i}N$ space steps. The equidistant nodes are thus of the form

$${\nu }_{ik}=-\frac{L}{2}+k\Delta {\nu }_i,k=0,1,\dots ,N_{i}N$$

(29)

where $\Delta {\nu }_i=\frac{L}{N_{i}N}$. The Nyquist relation2 holds whenever L is such that

$$Ll=2\pi N.$$

(30)

2.3. Implementation

In order to compute numerical approximations of Equations (22) and (23) at time node $t_i$, $i=0,1,\dots ,n-1$, we introduce the generic functions ${\theta }_i:\mathbb{R}\to \mathbb{R}$, $\psi :{\mathbb{R}}^2\to \mathbb{C}$ and ${\theta }_{i+1}:\mathbb{R}\to \mathbb{R}$ such that

$${\theta }_i\left(x\right)={\mathfrak{F}}^{-1}\left[\mathfrak{F}\left[{\theta }_{i+1}\right]\left(\nu \right)\psi (\nu ,x)\right]\left(x\right).$$

(31)

We assume that the function ${\theta }_{i+1}$ satisfies the periodicity boundary value equalities of Assumption 2.

Assumption 2.

The generic function ${\theta }_{i+1}$ satisfies

$$\begin{array}{cc}\hfill {\theta }_{i+1}\left(x_{i+1,0}\right)& ={\theta }_{i+1}\left(x_{i+1,N{N}_{i+1}}\right)\hfill \end{array}$$

(32)

$$\begin{array}{cc}\hfill \frac{\partial {\theta }_{i+1}}{\partial x}\left(x_{i+1,0}\right)& =\frac{\partial {\theta }_{i+1}}{\partial x}\left(x_{i+1,N{N}_{i+1}}\right).\hfill \end{array}$$

(33)

Hence, the Fourier integral

$$\mathfrak{F}\left[{\theta }_{i+1}\right]\left(\nu \right)={\int }_{-\infty }^{\infty }{e}^{-\mathbf{i}\nu x}{\theta }_{i+1}\left(x\right)dx$$

(34)

is restricted on the interval $[x_0-\frac{N_{i+1}l}{2},x_0+\frac{N_{i+1}l}{2}]=[x_{i+1,0},x_{i+1,N{N}_{i+1}}]$ and discretized using the grid points ${\left\{x_{i+1,k}\right\}}_{k=0}^{N_{i+1}N}$ with a quadrature rule with weights ${\left\{w_k\right\}}_{k=0}^{N_{i+1}N}$. As to the inverse Fourier integral of Equation (31) we restrict it on the interval $[-\frac{L}{2},\frac{L}{2}]$ and discretize it with lower Riemann sums at the Fourier space grid point ${\left\{{\nu }_{i+1,k}\right\}}_{k=0}^{N_{i+1}N}$.

Let $\mathcal{D}$ and ${\mathcal{D}}^{-1}$ denote the discrete Fourier transform and the inverse discrete Fourier transform, respectively

$$\begin{array}{cc}\hfill \mathfrak{D}{\left[{\left\{x\right\}}_{i=0}^{N-1}\right]}_k& =\frac{1}{N}\sum _{j=0}^{N-1}{e}^{-\mathbf{i}jk\frac{2\pi }{N}}{x}_j\hfill \end{array}$$

(35)

$$\begin{array}{cc}\hfill {\mathfrak{D}}^{-1}{\left[{\left\{x\right\}}_{i=0}^{N-1}\right]}_k& =\sum _{j=0}^{N-1}{e}^{\mathbf{i}jk\frac{2\pi }{N}}{x}_j.\hfill \end{array}$$

(36)

Then the discretization procedure leads to the approximation

$${\theta }_i\left(x_{i+1,k}\right)\approx {(-1)}^k{\mathfrak{D}}^{-1}{\left[{\left\{\psi ({\nu }_{i+1,j},x_{i+1,k})\mathbb{D}{\left[{\theta }_{i+1}\right]}_j\right\}}_{j=0}^{N_{i+1}N-1}\right]}_k$$

(37)

where

$$\mathbb{D}{\left[{\theta }_{i+1}\right]}_j=\mathfrak{D}{\left[{\left\{{(-1)}^s{\tilde{w}}_s{\theta }_{i+1}\left(x_{i+1,s}\right)\right\}}_{s=0}^{N_{i+1}N-1}\right]}_j$$

(38)

and the weights ${\left\{\tilde{w}\right\}}_{k=0}^{N_{i+1}N-1}$ are given by

$${\tilde{w}}_k=w_k+w_{N_{i+1}N}{\delta }_{k,0}.$$

(39)

where $\delta$ stands for the Kronecker delta. The relation of Equation (27) allows us to write

$${\theta }_i\left(x_{ik}\right)\approx {(-1)}^{k+\frac{N}{2}}{\mathfrak{D}}^{-1}{\left[{\left\{\psi ({\nu }_{i+1,j},x_{ik})\mathbb{D}{\left[{\theta }_{i+1}\right]}_j\right\}}_{j=0}^{N_{i+1}N-1}\right]}_{k+\frac{N}{2}}$$

(40)

for $k=0,1,\dots ,N_{i}N$.

In Equation (40), the generic function $\psi$ depends on the space node $x_{ik}$. If the relation generalizes for all space nodes $x_{ik}$, $k=0,1,\dots ,N_{i}N$, the function values ${\theta }_i\left(x_{ik}\right)$, $k=0,1,\dots ,N_{i}N$, can not be computed with a single direct FFT procedure. Instead, a separate FFT procedure using the values of the generic function $\psi$ at $x_{ik}$ is needed to compute the function value ${\theta }_i\left(x_{ik}\right)$. Nonetheless, the vector-matrix representation of the FFT procedure in Equation (40) allows the computation of all function values ${\theta }_i\left(x_{ik}\right)$ with a matrix multiplication. In the vector-matrix representation, Equation (40) is

$${\theta }_i\left(x_{ik}\right)={(-1)}^{k+\frac{N}{2}}{\widehat{F}}_{k+\frac{N}{2}}\mathsf{\Psi }\left(x_{ik}\right)\mathbb{D}\left[{\theta }_{i+1}\right]$$

where ${\widehat{F}}_{k+\frac{N}{2}}$ is the $(k+\frac{N}{2})$th row of the $N_{i+1}N$ dimension inverse FFT matrix $\widehat{F}$ and $\mathsf{\Psi }\left(x_{ik}\right)$ is the $N_{i+1}N$ dimension diagonal matrix built with the values ${\left\{\psi ({\nu }_{i+1,j},x_{ik})\right\}}_{j=0}^{N_{i+1}N-1}$. Let ${\mathsf{\Theta }}^{\left(i\right)}$ be the $N_{i}N$ dimension vector of the function values ${\theta }_i\left(x_{ik}\right)$ such that

$${\mathsf{\Theta }}_{1+k}^{\left(i\right)}={\theta }_i\left(x_{ik}\right)$$

(41)

for $k=0,1,\dots ,N_{i}N$. The matrix representation gives

$${\mathsf{\Theta }}^{\left(i\right)}={\widehat{\mathsf{\Psi }}}^{\left(i\right)}\mathbb{D}\left[{\theta }_{i+1}\right]$$

(42)

where ${\widehat{\mathsf{\Psi }}}^{\left(i\right)}$ is the $(N_{i}N+1)\times N_{i+1}N$ matrix such that

$${\widehat{\mathsf{\Psi }}}_{1+k,1+j}^{\left(i\right)}={(-1)}^{k+\frac{N}{2}}{\overline{\omega }}_i^{j(k+\frac{N}{2})}\psi ({\nu }_{i+1,j},x_{ik})$$

(43)

with ${\overline{\omega }}_i=e^{\mathbf{i}2\pi {\left(N_{i+1}N\right)}^{-1}}$, $k=0,1,\dots ,N_{i}N$ and $j=0,1,\dots ,N_{i+1}N-1$.

The requirements of Assumption 2 can easily be satisfied. Given a function $\eta :[a,b]\to \mathbb{R}$ and $\eta \in {\mathcal{C}}^1$, if we consider the transformation

$$\begin{array}{c}\hfill {\eta }^{\alpha ,\beta }\left(x\right)=\eta \left(x\right)+\alpha x^2+\beta x\end{array}$$

(44)

then the parameters $\alpha$ and $\beta$ can be chosen such that the transform function and its derivative have equal values at the boundaries of any interval. The following lemma gives a method to select the coefficients $\alpha$ and $\beta$ for the transform of Equation (44) such that Assumption 2 holds in general.

Lemma 1.

Suppose the real function $\eta \in {\mathcal{C}}^1[a,b]$ is differentiable and let ${\eta }^{\alpha ,\beta }$ be its transformed function as defined in Equation (44). Then

$$\begin{array}{cc}\hfill \alpha & =\frac{\frac{\partial \eta }{\partial x}\left(a\right)-\frac{\partial \eta }{\partial x}\left(b\right)}{2(b-a)},\hfill \end{array}$$

(45)

$$\begin{array}{cc}\hfill \beta & =\frac{\eta \left(a\right)-\eta \left(b\right)}{(b-a)}-\alpha (b+a)\hfill \end{array}$$

(46)

solve the system of linear equations defined by

$$\left\{\begin{array}{c}{\eta }^{\alpha ,\beta }\left(a\right)={\eta }^{\alpha ,\beta }\left(b\right)\hfill \\ \frac{\partial {\eta }^{\alpha ,\beta }}{\partial x}\left(a\right)=\frac{\partial {\eta }^{\alpha ,\beta }}{\partial x}\left(b\right).\hfill \end{array}\right.$$

(47)

Proof. 

The second equation of the system (47) gives Equation (45) in a straightforward manner. Equation (46) is given by the first equation of the system.    □

Hence, the numerical discretization may be applied on the transformation $u_{i+1}^{\alpha ,\beta }$ at time node $t_i$ but a correction must be performed so to recover the values of the intermediate solution ${\tilde{u}}_i$ and the approximate gradient ${\dot{u}}_i$. The next theorem gives the representation under the transform of Equation (44).

Theorem 1.

Let $u_{i+1}^{\alpha ,\beta }$ be the alternative transform defined in Equation (44) of the approximate solution $u_{i+1}$. Then the intermediate solution ${\tilde{u}}_i$ and the approximate gradient ${\dot{u}}_i$ in Equations (22) and (23) satisfy

$$\begin{array}{ll}{\tilde{u}}_i\left(x\right)\hfill & ={\mathfrak{F}}^{-1}\left[\mathfrak{F}\left[u_{i+1}^{\alpha ,\beta }\right]\left(\nu \right)\varphi (\nu ,x)\right]\left(x\right)-\alpha [{(x+{\Delta }_{i}a(t_i,x))}^2+{\Delta }_i{\sigma }^2(t_i,x)]\hfill \\ & -\beta (x+{\Delta }_{i}a(t_i,x))\end{array}$$

(48)

$$\begin{array}{cc}{\dot{u}}_i\left(x\right)& =\sigma (t_i,x){\mathfrak{F}}^{-1}\left[\mathfrak{F}\left[u_{i+1}^{\alpha ,\beta }\right]\left(\nu \right)\mathbf{i}\nu \varphi (\nu ,x)\right]\left(x\right)-\sigma (t_i,x)[2\alpha (x+{\Delta }_{i}a(t_i,x))+\beta ].\end{array}$$

(49)

Proof. 

The proof follows the steps of Theorem 3.1 of Hyndman and Oyono Ngou (2017) using the transformation introduced in Equation (44) and the relations (20) and (21) between the density function $h_i$ and the characteristic function ${\varphi }_i$.    □

Algorithm 1 details the numerical procedure on the space grid and produces numerical solutions ${\left\{u_{ik}\right\}}_{k=0}^{N_{i}N}$, ${\left\{{\tilde{u}}_{ik}\right\}}_{k=0}^{N_{i}N}$ and ${\left\{{\dot{u}}_{ik}\right\}}_{k=0}^{N_{i}N}$ for the approximate solution $u_i$, the intermediate solution ${\tilde{u}}_i$ and the approximate gradient ${\dot{u}}_i$, respectively, $i=0,1,\dots ,n-1$. We next consider error estimates under the alternative discretization.

Algorithm 1 Fourier Interpolation Method on Alternative Grid

Discretize the restricted real space $[x_0-\frac{N_{n}l}{2},x_0+\frac{N_{n}l}{2}]$ and the restricted Fourier space $[-\frac{L}{2},\frac{L}{2}]$ with $N_{n}N$ space steps so to have the real space nodes ${\left\{x_{nk}\right\}}_{k=0}^{N_{n}N}$ and ${\left\{{\nu }_{nk}\right\}}_{k=0}^{N_{n}N}$ Value $u_n\left(x_{nk}\right)=g\left(x_{nk}\right)$ For any i from $n-1$ to 0 (a) Compute $\alpha$ and $\beta$ defining the transform of Equation (44), such that $${\theta }_{i+1}=u_{i+1}^{\alpha ,\beta }$$ (50) and ${\theta }_{i+1}$ satisfies the boundary conditions of Equations (32) and (33). (b) Compute ${\theta }_i\left(x_{ik}\right)$ through Equation (40) for $k=0,1,\dots ,N_{i}N$ with $$\psi (\nu ,x)={\varphi }_i(\nu ,x)$$ (51) and retrieve the values ${\tilde{u}}_{ik}$ as $$\begin{array}{cc}\hfill {\tilde{u}}_{ik}& ={\theta }_i\left(x_{ik}\right)-\alpha [{(x_{ik}+{\Delta }_{i}a(t_i,x_{ik}))}^2+{\Delta }_i{\sigma }^2(t_i,x_{ik})]\hfill \\ \hfill & -\beta (x_{ik}+{\Delta }_{i}a(t_i,x_{ik})).\hfill \end{array}$$ (52) (c) Compute ${\theta }_i\left(x_{ik}\right)$ through Equation (40) for $k=0,1,\dots ,N_{i}N$ with $$\psi (\nu ,x)=\mathbf{i}\nu \sigma (\nu ,x){\varphi }_i(\nu ,x)$$ (53) and retrieve the values ${\dot{u}}_{ik}$ as $${\dot{u}}_{ik}={\theta }_i\left(x_{ik}\right)-\sigma (t_i,x_{ik})[2\alpha (x_{ik}+{\Delta }_{i}a(t_i,x_{ik}))+\beta ].$$ (54) (d) Compute the values $u_{ik}$ as $$u_{ik}={\tilde{u}}_{ik}+{\Delta }_{i}f(t_i,x_{ik},{\tilde{u}}_{ik},{\dot{u}}_{ik})$$ (55) for $k=0,1,\dots ,N_{i}N$ through Equation (13). (e) Update the real space grid with Equation (27) and the Fourier space grid by discretizing the interval $[-\frac{L}{2},\frac{L}{2}]$ with $N_{i}N$ space steps so to have the real space nodes ${\left\{x_{ik}\right\}}_{k=0}^{N_{i}N}$ and ${\left\{{\nu }_{ik}\right\}}_{k=0}^{N_{i}N}$.

2.4. Spatial Discretization Error Analysis

Let ${\left\{{\mathbf{u}}_{ik}\right\}}_{k=0}^{N_{i}N}$, ${\left\{{\tilde{\mathbf{u}}}_{ik}\right\}}_{k=0}^{N_{i}N}$ and ${\left\{{\dot{\mathbf{u}}}_{ik}\right\}}_{k=0}^{N_{i}N}$ denote the numerical solutions obtained from the convolution method at time node $t_i$ given the solution $u_{i+1}$ at time $t_{i+1}$. For the Fourier interpolation method on the alternative grid, we defined the local discretization error as

$$E_{ik}:=\left|u_i\left(x_k\right)-{\mathbf{u}}_{ik}\right|+\left|{\dot{u}}_i\left(x_k\right)-{\dot{\mathbf{u}}}_{ik}\right|$$

(56)

for $i=0,1,\dots ,n-1$ and $k=0,1,\dots ,N_{i}N$.

Theorem 2.

Suppose that the driver is $f\in {\mathcal{C}}^{1,2}([0,T]\times {\mathbb{R}}^3)$, the terminal condition is $g\in {\mathcal{C}}^2\left(\mathbb{R}\right)$, and Assumption 1 is satisfied. Then the convolution method yields a local discretization error of the form

$$E_{ik}=\mathcal{O}\left(\Delta x\right)+\mathcal{O}\left(e^{-K{\left|{\Delta }_i\right|}^{-1}{l}^2}\right)$$

(57)

for some constant $K>0$ on the alternative grid and under the trapezoidal quadrature rule with weights

$$w_j=1-\frac{1}{2}({\delta }_{j,0}+{\delta }_{j,N_{i+1}N}).$$

Proof. 

We suppose the solution $u_{i+1}$ at time $t_{i+1}$ is known. The solution $u_{i+1}\in {\mathcal{C}}^2$ is twice differentiable since $f\in {\mathcal{C}}^{1,2}$ and $g\in {\mathcal{C}}^2$. Furthermore, $u_{i+1}$ is square integrable with respect to the Gaussian density.

By Theorem 1, we limit ourselves to the case where

$$\begin{array}{cc}\hfill u_{i+1}\left(x_0-\frac{N_{i+1}l}{2}\right)& =u_{i+1}\left(x_0+\frac{N_{i+1}l}{2}\right)\hfill \\ \hfill \frac{\partial u_{i+1}}{\partial x}\left(x_0-\frac{N_{i+1}l}{2}\right)& =\frac{\partial u_{i+1}}{\partial x}\left(x_0+\frac{N_{i+1}l}{2}\right)\hfill \end{array}$$

so that the coefficients of the transform are $\alpha =\beta =0$. Let $T_i$ be the Fourier polynomial interpolating $u_{i+1}$ on $\left[x_0-\frac{N_{i+1}l}{2},x_0+\frac{N_{i+1}l}{2}\right]$. Then

$$\begin{array}{cc}\hfill T_i\left(x\right)& :=\sum _{k=-\frac{N_{i+1}N}{2}}^{\frac{N_{i+1}N}{2}-1}{d}_j{e}^{\mathbf{i}k\frac{2\pi }{N_{i+1}l}x}\hfill \end{array}$$

(58)

$$\begin{array}{cc}\hfill & =u_{i+1}\left(x\right)+\mathcal{O}(\Delta x),\forall x\in \left[x_0-\frac{N_{i+1}l}{2},x_0+\frac{N_{i+1}l}{2}\right]\hfill \end{array}$$

(59)

where

$${(-1)}^{j-\frac{N_{i+1}N}{2}}{d}_{j-\frac{N_{i+1}N}{2}}=\mathbb{D}{\left[u_{i+1}\right]}_j,j=0,1,\dots ,N_{1+i}N-1$$

(60)

when using the trapezoidal quadrature rule. We have that

$$\begin{array}{cc}\hfill {\tilde{u}}_i\left(x_{ik}\right)& ={\int }_{\left|y\right|\le \frac{l}{2}}{u}_{i+1}(x_{ik}+y)h_i\left(y\right|x_{ik})dy+{\int }_{\left|y\right|>\frac{l}{2}}{u}_{i+1}(x_{ik}+y)h_i\left(y\right|x_{ik})dy.\hfill \end{array}$$

Note that

$$\begin{array}{cc}\hfill {\int }_{\left|y\right|>\frac{l}{2}}{h}_i\left(y\right|x_{ik})dy& =\mathbf{P}\left[|\Delta X_i^{\pi }|>\frac{l}{2}|X_{t_i}^{\pi }=x_{ik}\right]\hfill \\ & =\mathbf{P}\left[{\left|\frac{\Delta X_i^{\pi }}{\sigma (t,x_{ik})\sqrt{{\Delta }_i}}\right|}^2>\frac{l^2}{4{\sigma }^2(t_i,x_{ik}){\Delta }_i}|X_{t_i}^{\pi }=x_{ik}\right]\hfill \end{array}$$

and that the random variable ${\left(\frac{\Delta X_i^{\pi }}{\sigma (t,x_{ik})\sqrt{{\Delta }_i}}\right)}^2$ has a non-central chi-square distribution with 1 degree of freedom and non-centrality parameter $\lambda ={\left(\frac{a(t_i,x_{ik})}{\sigma (t_i,x_{ik})}\right)}^2{\Delta }_i$. Then for some constant $K>0$ which is inversely proportional to ${\Delta }_i$ we have

$$\begin{array}{cc}\hfill {\int }_{\left|y\right|>\frac{l}{2}}{u}_{i+1}(x_{ik}+y)h_i\left(y\right|x_{ik})dy& =\mathcal{O}\left(e^{-K{l}^2}\right)\hfill \end{array}$$

by the Cauchy-Schwarz and Chernoff inequalities since the solution $u_{i+1}$ is square integrable. Hence

$$\begin{array}{cc}\hfill & {\tilde{u}}_i\left(x_{ik}\right)={\int }_{\left|y\right|\le \frac{l}{2}}{u}_{i+1}(x_{ik}+y)h_i\left(y\right|x_{ik})dy+\mathcal{O}\left(e^{-K{l}^2}\right)\hfill \\ \hfill & ={\int }_{\left|y\right|\le \frac{l}{2}}{T}_i(x_{ik}+y)h_i\left(y\right|x_{ik})dy+\mathcal{O}\left(\Delta x\right)+\mathcal{O}\left(e^{-K{l}^2}\right)(\mathrm{by}\mathrm{Equation}(59))\hfill \\ \hfill & ={\int }_{\mathbb{R}}{T}_i(x_{ik}+y)h_i\left(y\right|x_{ik})dy+\mathcal{O}\left(\Delta x\right)+\mathcal{O}\left(e^{-K{l}^2}\right)\hfill \\ \hfill & \hfill (\mathrm{by}\mathrm{Chernoff}\mathrm{inequality},\mathrm{since}{T}_i\mathrm{is}\mathrm{bounded})\\ \hfill & ={\int }_{\mathbb{R}}\sum _{j=-\frac{N_{i+1}N}{2}}^{\frac{N_{i+1}N}{2}-1}{d}_j{e}^{\mathbf{i}j\frac{2\pi }{N_{i+1}l}(x_{i,k}+y)}{h}_i\left(y\right|x_{ik})dy+\mathcal{O}(\Delta x)+\mathcal{O}\left(e^{-K{l}^2}\right)\hfill \\ \hfill & =\sum _{j=-\frac{N_{i+1}N}{2}}^{\frac{N_{i+1}N}{2}-1}{d}_j{e}^{\mathbf{i}j\frac{2\pi }{N_{i+1}l}{x}_{i,k}}{\varphi }_i\left(j\frac{2\pi }{N_{i+1}l},x_{ik}\right)+\mathcal{O}(\Delta x)+\mathcal{O}\left(e^{-K{l}^2}\right)\hfill \\ \hfill & ={(-1)}^{k+\frac{N}{2}}\sum _{j=0}^{N_{i+1}N-1}\varphi ({\nu }_{i+1,j},x_{ik})\mathbb{D}{\left[u_{i+1}\right]}_j{e}^{\mathbf{i}\frac{2\pi }{N_{i+1}N}j(k+\frac{N}{2})}+\mathcal{O}(\Delta x)\hfill \\ & +\mathcal{O}\left(e^{-K{l}^2}\right)(\mathrm{by}\mathrm{Equation} (60)\mathrm{when}\mathrm{using}\mathrm{the}\mathrm{trapezoidal}\mathrm{quadrature}\mathrm{rule})\hfill \\ \hfill & ={\tilde{\mathbf{u}}}_{ik}+\mathcal{O}(\Delta x)+\mathcal{O}\left(e^{-K{l}^2}\right).\hfill \end{array}$$

Similar techniques show that

$${\dot{u}}_i\left(x_k\right)={\dot{\mathbf{u}}}_{ik}+\mathcal{O}\left(\Delta x\right)+\mathcal{O}\left(e^{-K{l}^2}\right)$$

(61)

where $K>0$ is inversely proportional to ${\Delta }_i$. The Lipschitz property of the driver f completes the proof.    □

As expected, the alternative discretization improves the local error bound by eliminating extrapolation errors in Hyndman and Oyono Ngou (2017). The result of Theorem 2 establishes the consistency of the convolution method with respect to the approximate functions $u_i$ and gradients ${\dot{u}}_i$. Furthermore, the absence of extrapolation errors in the local discretization allows us to develop a tighter bound for the global discretization error. The following corollary proves helpful when deriving the global discretization error bound.

Corollary 1.

Under the conditions of Theorem 2, there is $C>0$ such that

$$\underset{i,k}{sup}{E}_{i,k}=\mathcal{O}(\Delta x)+\mathcal{O}\left(e^{-C{\left|\pi \right|}^{-1}{l}^2}\right).$$

(62)

We define the global error as

$$E_{l,\Delta x}:=\underset{i,k}{sup}{e}_{ik}+\underset{i,k}{sup}{\dot{e}}_{ik}$$

(63)

where

$$\begin{array}{cc}\hfill e_{ik}& =\left|u_{n-i}\left(x_k\right)-u_{n-i,k}\right|\hfill \end{array}$$

(64)

$$\begin{array}{cc}\hfill {\dot{e}}_{ik}& =\left|{\dot{u}}_{n-i}\left(x_k\right)-{\dot{u}}_{n-i,k}\right|\hfill \end{array}$$

(65)

for $i=1,\dots ,n$ with $e_{0,k}={\dot{e}}_{0,k}=0$. The next theorem describes the stability and convergence properties of the convolution method.

Theorem 3.

Suppose the conditions of Theorem 2 are satisfied. If the space discretization is such that

$$\underset{i}{sup}max\left(\frac{K_4^{\frac{1}{2}}\Delta x}{\sqrt{2\pi {\Delta }_i}},\frac{K_4\Delta x}{\pi {\Delta }_i}\right)\le 1$$

(66)

then the Fourier interpolation method is stable and the global discretization error $E_{l,\Delta x}$ satisfies

$$E_{l,\Delta x}=\mathcal{O}(\Delta x)+\mathcal{O}\left(e^{-C{\left|\pi \right|}^{-1}{l}^2}\right)$$

(67)

where $C>0$ and $K_4$ is the upper bound of Equation (7).

Proof. 

Let us first notice that

$$\begin{array}{cc}\hfill e_{ik}& \le E_{n-i,k}+\left|{\mathbf{u}}_{n-i,k}-u_{n-i,k}\right|\hfill \\ \hfill & \le E_{n-i,k}+(1+{\Delta }_{i}K)\left|{\tilde{\mathbf{u}}}_{n-i,k}-{\tilde{u}}_{n-i,k}\right|+{\Delta }_{i}K\left|{\dot{\mathbf{u}}}_{n-i,k}-{\dot{u}}_{n-i,k}\right|\hfill \end{array}$$

(68)

where $K>0$ is the Lipschitz constant of the driver f. Furthermore, we have that

$${\dot{e}}_{ik}\le E_{n-i,k}+\left|{\dot{\mathbf{u}}}_{n-i,k}-{\dot{u}}_{n-i,k}\right|.$$

(69)

Furthermore, the construction of the Fourier interpolation method gives

$$\begin{array}{cc}\hfill \left|{\tilde{\mathbf{u}}}_{i,k}-{\tilde{u}}_{i,k}\right|& \le \left|{\mathfrak{D}}^{-1}{\left[{\left\{\varphi ({\nu }_{i+1,j},x_{ik})\mathbb{D}{[u_{i+1}-u_{i+1,s}]}_j\right\}}_{j=0}^{N_{i+1}N-1}\right]}_{k+\frac{N}{2}}\right|\hfill \\ \hfill & \le \frac{1}{N_{i+1}N}\left(\sum _{j=0}^{N_{i+1}N-1}\left|\varphi ({\nu }_{i+1,j},x_{ik})\right|\right)\underset{k}{sup}\left|u_{i+1}\left(x_{ik}\right)-u_{i+1,k}\right|\hfill \\ \hfill & \hfill (\mathrm{using}\mathrm{the}\mathrm{matrix}-\mathrm{vector}\mathrm{representation}\mathrm{of}\mathrm{DFTs})\\ \hfill & \le \frac{1}{N_{i+1}N}\left(\sum _{j=0}^{N_{i+1}N-1}\left|\varphi ({\nu }_{i+1,j},x_{ik})\right|\right)\underset{k}{sup}{e}_{n-i-1,k}\hfill \\ \hfill & \le \frac{{(\Delta {\nu }_{i+1})}^{-1}}{N_{i+1}N}\left({\int }_{\mathbb{R}}\left|\varphi (\nu ,x_{ik})\right|d\nu \right)\underset{k}{sup}{e}_{n-i-1,k}\hfill \end{array}$$

$$\begin{array}{cc}\hfill & \le \frac{K_4^{\frac{1}{2}}\Delta x}{{\left(2\pi {\Delta }_i\right)}^{\frac{1}{2}}}{\displaystyle \underset{k}{sup}}{e}_{n-i-1,k}\hfill \end{array}$$

(70)

where the last inequality holds by Assumption 1. Similarly,

$$\begin{array}{cc}\hfill \left|{\dot{\mathbf{u}}}_{i,k}-{\dot{u}}_{i,k}\right|& \le \left|{\mathfrak{D}}^{-1}{\left[{\left\{\psi ({\nu }_{i+1,j},x_{ik})\mathbb{D}{[u_{i+1}-u_{i+1,s}]}_j\right\}}_{j=0}^{N_{i+1}N-1}\right]}_{k+\frac{N}{2}}\right|\hfill \\ \hfill & \le \frac{1}{N_{i+1}N}\left(\sum _{j=0}^{N_{i+1}N-1}\left|\mathbf{i}{\nu }_{i+1,j}\varphi ({\nu }_{i+1,j},x_{ik})\right|\right)\underset{k}{sup}{e}_{n-i-1,k}\hfill \\ \hfill & \hfill (\mathrm{using}\mathrm{the}\mathrm{matrix}-\mathrm{vector}\mathrm{representation}\mathrm{of}\mathrm{DFTs})\\ \hfill & \le \frac{{(\Delta {\nu }_{i+1})}^{-1}}{N_{i+1}N}\left({\int }_{\mathbb{R}}\left|\nu \varphi (\nu ,x_{ik})\right|d\nu \right)\underset{k}{sup}{e}_{n-i-1,k}\hfill \end{array}$$

$$\begin{array}{cc}\hfill & \le \frac{K_4\Delta x}{\pi {\Delta }_i}{\displaystyle \underset{k}{sup}}{e}_{n-i-1,k}.\hfill \end{array}$$

(71)

The inequalities of Equations (68), (70) and (71) lead to

$$\begin{array}{cc}\hfill e_{i,k}& \le C_0{E}_{i,k}+\left(1+2{\Delta }_{i}K\right)max\left(\frac{K_4^{\frac{1}{2}}\Delta x}{\sqrt{2\pi {\Delta }_i}},\frac{K_4\Delta x}{\pi {\Delta }_i}\right)\underset{k}{sup}{e}_{i-1,k}\hfill \\ \hfill & \le C_0\underset{i,k}{sup}{E}_{i,k}+\left(1+2{\Delta }_{i}K\right)max\left(\frac{K_4^{\frac{1}{2}}\Delta x}{\sqrt{2\pi {\Delta }_i}},\frac{K_4\Delta x}{\pi {\Delta }_i}\right)\underset{k}{sup}{e}_{i-1,k}\hfill \end{array}$$

where $C_0>0$ and $K>0$ is the Lipschitz constant of the driver f. Consequently,

$$\begin{array}{cc}\hfill \underset{k}{sup}{e}_{i,k}& \le C_0\underset{i,k}{sup}{E}_{i,k}+\left(1+2{\Delta }_{i}K\right)max\left(\frac{K_4^{\frac{1}{2}}\Delta x}{\sqrt{2\pi {\Delta }_i}},\frac{K_4\Delta x}{\pi {\Delta }_i}\right)\underset{k}{sup}{e}_{i-1,k}\hfill \\ \hfill & \le C_0\underset{i,k}{sup}{E}_{i,k}+\left(1+2{\Delta }_{i}K\right)\underset{k}{sup}{e}_{i-1,k}\hfill \end{array}$$

(72)

since

$$\underset{i}{sup}max\left(\frac{K_4^{\frac{1}{2}}\Delta x}{\sqrt{2\pi {\Delta }_i}},\frac{K_4\Delta x}{\pi {\Delta }_i}\right)\le 1$$

and the Gronwall’s Lemma yields

$$\underset{k}{sup}{e}_{i,k}\le C_0{e}^{2TK}\underset{i,k}{sup}{E}_{i,k}$$

(73)

from the inequality of Equation (72) for $i=0,1,\dots ,n$ knowing that $e_{0,k}=0$. The last equation establishes the stability of the Fourier interpolation method for the approximate solution $u_i$ since its error at any time step is absolutely bounded.

The inequalities of Equations (69), (71) and (73) lead to

$$\begin{array}{cc}\hfill \underset{k}{sup}{\dot{e}}_{i,k}& \le \left(C_1+\frac{\Delta x}{\pi {\Delta }_i}{C}_0{e}^{2TK}\right)\underset{i,k}{sup}{E}_{i,k}\hfill \\ \hfill & \le \left(C_1+C_0{e}^{2TK}\right)\underset{i,k}{sup}{E}_{i,k}\hfill \end{array}$$

(74)

for a positive constant $C_1>0$. Hence, the convolution method is also stable for the approximate gradient ${\dot{u}}_i$. The result of Equation (67) follows by taking the supremum on the left hand sides of Equations (73) and (74) over time steps and applying Corollary 1.    □

As for most explicit methods for PDEs, the convolution method requires a stability condition as described in Equation (66). In general, Theorem 3 shows that the convolution method converges when the space discretization is relatively as fine as the time discretization. Other numerical methods for BSDEs, and particularly Monte Carlo based methods, have a stability and convergence condition. Indeed, error explosion occurs for fine time discretizations in the backward methods of Gobet et al. (2005) and Bouchard and Touzi (2004). In order to maintain stability and convergence, the space discretization has to be refined by increasing the number of simulated paths.

3. Higher Order Time Discretization for FBSDEs

In this section, we discuss further extensions of the Fourier interpolation method on the alternative grid. In particular, we apply the Fourier interpolation method to Runge–Kutta schemes for FBSDEs proposed by Chassagneux and Crisan (2014). Rather than computing the numerical approximation $(Y_{t_i}^{\pi },Z_{t_i}^{\pi })$ recursively backward in time for $i=0,\dots ,n$ the Runge–Kutta method introduced in Chassagneux and Crisan (2014) for decoupled FBSDEs uses p intermediate stages implemented between time partition points.

3.1. Runge–Kutta Schemes

The FBSDE of Equation (1) is discretized on the time partition $\pi$. Following Definition 1.1 of Chassagneux and Crisan (2014) let $q\in {\mathbb{N}}^{*}$, we consider the q-stage Runge–Kutta scheme giving the following numerical solution at mesh time $t_i$

$$\begin{array}{cc}\hfill Z_{t_i}^{\pi }& ={\mathbf{E}}_{t_i}^{}\left[H_{t_i,{\Delta }_i}^{{\phi }_1}{Y}_{t_{i+1}}^{\pi }+{\Delta }_i\sum _{j=1}^q{\beta }_j{H}_{t_i,(1-{\gamma }_j){\Delta }_i}^{{\phi }_1}f(t_{i,j},X_{t_{i,j}}^{\pi },Y_{i,j}^{\pi },Z_{i,j}^{\pi })\right]\hfill \end{array}$$

(75)

$$\begin{array}{cc}\hfill Y_{t_i}^{\pi }& ={\mathbf{E}}_{t_i}^{}\left[Y_{t_{i+1}}^{\pi }+{\Delta }_i\sum _{j=1}^{q+1}{\alpha }_{j}f(t_{i,j},X_{t_{i,j}}^{\pi },Y_{i,j}^{\pi },Z_{i,j}^{\pi })\right]\hfill \end{array}$$

(76)

for a set positive coefficients ${\left\{{\gamma }_j\right\}}_{j=1}^{q+1}$ such that $0={\gamma }_1<\dots <{\gamma }_{q+1}=1$. The intermediate solutions ${\left\{(Y_{i,j}^{\pi },Z_{i,j}^{\pi })\right\}}_{j=2}^q$ take the form

$$\begin{array}{cc}\hfill Z_{i,j}^{\pi }& ={\mathbf{E}}_{t_{i,j}}^{}\left[H_{t_{i,j},{\gamma }_j{\Delta }_i}^{{\phi }_j}{Y}_{t_{i+1}}^{\pi }+{\Delta }_i\sum _{k=1}^{j-1}{\beta }_{jk}{H}_{t_{i,j},({\gamma }_j-{\gamma }_k){\Delta }_i}^{{\phi }_j}f(t_{i,k},X_{t_{i,k}}^{\pi },Y_{i,k}^{\pi },Z_{i,k}^{\pi })\right]\hfill \end{array}$$

(77)

$$\begin{array}{cc}\hfill Y_{i,j}^{\pi }& ={\mathbf{E}}_{t_{i,j}}^{}\left[Y_{t_{i+1}}^{\pi }+{\Delta }_i\sum _{k=1}^j{\alpha }_{jk}f(t_{i,k},X_{t_{i,k}}^{\pi },Y_{i,k}^{\pi },Z_{i,k}^{\pi })\right]\hfill \end{array}$$

(78)

where

$$t_{i,j}=t_i+(1-{\gamma }_j){\Delta }_i,1\le j\le q+1$$

(79)

with $(Y_{i,1}^{\pi },Z_{i,1}^{\pi })=(Y_{t_{i+1}}^{\pi },Z_{t_{i+1}}^{\pi })$, $(Y_{i,q+1}^{\pi },Z_{i,q+1}^{\pi })=(Y_{t_i}^{\pi },Z_{t_i}^{\pi })$ and terminal condition

$$(Y_{t_n},Z_{t_n})=(g\left(X_T\right),{\sigma }^{*}(T,X_T)\nabla g\left(X_T\right)).$$

(80)

The coefficients ${\left\{{\alpha }_j\right\}}_{j=1}^{q+1}$, ${\left\{{\beta }_j\right\}}_{j=1}^q$, $\{{\alpha }_{jk}:1\le j\le q,1\le k\le j\}$ and $\{{\beta }_{jk}:1\le j\le q,1\le k<j\}$ are all positive and satisfy

$$\begin{array}{cc}\hfill \sum _{j=1}^{q+1}{\alpha }_j& =1\hfill \end{array}$$

(81)

$$\begin{array}{cc}\hfill {\beta }_{jj}& =0,1\le j\le q,\hfill \end{array}$$

(82)

$$\begin{array}{cc}\hfill \sum _{k=1}^j{\alpha }_{jk}& =\sum _{k=1}^{j-1}{\beta }_{jk}={\gamma }_j,1<j\le q.\hfill \end{array}$$

(83)

Let ${\mathcal{B}}^m$ denote the set of continuous and bounded functions on $[0,1]$ such that

$${\mathcal{B}}^m:=\left\{\varphi \in {\mathcal{C}}_b:{\int }_0^1{s}^k\varphi \left(s\right)ds={\delta }_{0,k},k\le m\mathrm{and}k,m\in {\mathbb{N}}^{*}\right\}.$$

(84)

The stochastic coefficient $H_{t,\Delta }^{\phi }$ with $t\in [0,T)$ and $\Delta >0$ is defined as

$$H_{t,\Delta }^{\phi }:=\frac{1}{\Delta }{\int }_t^{t+\Delta }\phi \left(\frac{s-t}{\Delta }\right)d{W}_s$$

(85)

with $\phi \in {\mathcal{B}}^m$ for some $m\in {\mathbb{N}}^{*}$.

The above definitions give the $q-$stage Runge–Kutta and the method has a number of advantages. The global error of the $q-$stage Runge–Kutta scheme ${\mathcal{E}}_{\pi }$ is defined as

$$\begin{array}{cc}\hfill {\mathcal{E}}_{\pi }^2& :=\underset{0\le i<n}{max}\parallel Y_{t_i}-Y_{t_i}^{\pi }{\parallel }_{L^2}^2+\sum _{i=0}^{n-1}{\Delta }_i{\parallel Z_{t_i}-Z_{t_i}^{\pi }\parallel }_{L^2}^2\hfill \\ \hfill & =\underset{0\le i<n}{max}{\mathbf{E}}_{}^{}\left[{\left|Y_{t_i}-Y_{t_i}^{\pi }\right|}^2\right]+\sum _{i=0}^{n-1}{\Delta }_i{\mathbf{E}}_{}^{}\left[{\left|Z_{t_i}-Z_{t_i}^{\pi }\right|}^2\right]\hfill \end{array}$$

(86)

and is hence weaker than the error $E_{\pi }$ considered for the Euler scheme. Nonetheless, the global error ${\mathcal{E}}_{\pi }$ is easier to handle since it is strongly related to the local time discretization error which simplifies the theoretical study in Chassagneux and Crisan (2014).

The scheme can be represented by the following tableau

${\gamma }_1$	${\alpha }_{1,1}$	0	…	0	0	${\beta }_{1,1}$	0	…	0
${\gamma }_2$	${\alpha }_{2,1}$	${\alpha }_{2,2}$	…	0	0	${\beta }_{2,1}$	${\beta }_{2,2}$	…	0
⋮	⋮	⋮	⋱	⋮	⋮	⋮	⋮	⋱	⋮
${\gamma }_q$	${\alpha }_{q,1}$	${\alpha }_{q,2}$	…	${\alpha }_{q,q}$	0	${\beta }_{q,1}$	${\beta }_{q,2}$	…	${\beta }_{q,q}$
${\gamma }_{q+1}$	${\alpha }_1$	${\alpha }_2$	…	${\alpha }_q$	${\alpha }_{q+1}$	${\beta }_1$	${\beta }_2$	…	${\beta }_q$

One can observe that if ${\alpha }_{q+1}=0$ and ${\alpha }_{jj}=0$, $1<j\le q$, then the q-stage Runge–Kutta scheme is explicit. Otherwise, the scheme is implicit. For instance, the Runge–Kutta schemes with tableau

0	0	0	0
1	0	1	1

and the scheme with tableau

0	0	0	0
1	$\frac{1}{2}$	$\frac{1}{2}$	1

known as the Crank-Nicolson scheme constitute $1-$stage implicit Runge–Kutta schemes. The only $1-$stage explicit Runge–Kutta scheme admits the tableau

0	0	0	0
1	1	0	1

In Chassagneux and Crisan (2014), the implicit and the explicit $1-$stage Runge–Kutta schemes are shown to be at least one-half $\left(\frac{1}{2}\right)$ order convergent. The Crank-Nicolson scheme, already studied in Crisan and Manolarakis (2014), presents a first-order of convergence. Notice that the Euler schemes used in the previous section are not $1-$stage Runge–Kutta schemes since they do not lead to any consistent tableau. Nonetheless, their structure is equivalent to the explicit $1-$stage Runge–Kutta scheme and both schemes display the same half $\left(\frac{1}{2}\right)$ order of convergence. The following tableau gives a example of explicit 2-stage Runge–Kutta schemes of first-order of convergence for ${\gamma }_2\in (0,1]$ and ${\beta }_1\in [0,1]$.

0	0	0	0	0	0
${\gamma }_2$	${\gamma }_2$	0	0	${\gamma }_2$	0
1	$1-\frac{1}{2{\gamma }_2}$	$\frac{1}{2{\gamma }_2}$	0	${\beta }_1$	$1-{\beta }_1$

3.2. Further Simplification

From the q-stage Runge–Kutta scheme for BSDEs, one notices that we have at least $2q$ conditional expectations to compute at each time step. These conditional expectations can be simplified and made more suitable for numerical implementation if we consider a reasonable time discretization of the forward SDE. Hence, we make the following assumption.

Assumption 3

(Forward process discretization). The following are assumed throughout this section.

1. 

The forward SDE is discretized with the piecewise constant process $X^{\pi }$ such that for $t\in [t_i,t_{i+1})$ we have $X_t^{\pi }=X_{t_i}^{\pi }$ pathwise.

2. 

The forward SDE time discretization with global error ${\mathcal{E}}_{X,\pi }$ is of order $m>0$, i.e.,

$${\mathcal{E}}_{X,\pi }^2:=\underset{0\le i\le n}{max}{\parallel X_{t_i}-X_{t_i}^{\pi }\parallel }_{L^2}^2=\mathcal{O}\left({\left|\pi \right|}^{2m}\right).$$

(87)

3. 

The forward SDE time discretization admits the conditional characteristic functions ${\varphi }_i:{\mathbb{R}}^d\times {\mathbb{R}}^d\to \mathbb{C}$

$${\varphi }_i(\nu ,x)={\mathbf{E}}_{}^{}\left[e^{\mathbf{i}{\nu }^{*}\left(X_{t_{i+1}}^{\pi }-X_{t_i}^{\pi }\right)}|X_{t_i}^{\pi }=x\right]$$

(88)

and ${\mathsf{\Phi }}_{i,j}:{\mathbb{R}}^d\times {\mathbb{R}}^d\to {\mathbb{C}}^d$

$${\mathsf{\Phi }}_{i,j}(\nu ,x)={\mathbf{E}}_{}^{}\left[H_{t_{i,j},{\gamma }_j{\Delta }_i}^{{\phi }_j}{e}^{\mathbf{i}{\nu }^{*}\left(X_{t_{i+1}}^{\pi }-X_{t_i}^{\pi }\right)}|X_{t_i}^{\pi }=x\right]$$

(89)

for $0\le i<n$ and $1<j\le q+1$ with ${\phi }_{q+1}={\phi }_1$.

4. 

There are positive constants $p_0$, $q_0$, $s_0$, $K_0$ and $C_0>0$ such that

$$max\left(\underset{s\in {\mathbb{R}}_d^{+}}{inf}{e}^{-s*t}{\varphi }_i(\mathbf{i}s,x),\underset{s\in {\mathbb{R}}_d^{+}}{inf}{e}^{-s*t}{\varphi }_i(-\mathbf{i}s,x)\right)\le e^{-K_0{\Delta }_i^{-s_0}{\left|t\right|}^{q_0}},$$

(90)

$\forall t\in {\mathbb{R}}_d^{+}$, hence the discrete version of the forward process has conditional exponential moments. In addition,

$${\int }_{{\mathbb{R}}^d}\left|{\varphi }_i(\nu ,x)\right|d\nu +\underset{1<j\le q+1}{max}{\int }_{{\mathbb{R}}^d}\left|{\mathsf{\Phi }}_{i,j}(\nu ,x)\right|d\nu \le C_0{\Delta }_i^{-p_0}.$$

(91)

Itô-Taylor expansion based schemes are an example of SDE discretization satisfying the conditions of Assumption 3. A more complete presentation of these schemes can be found in Kloeden and Platen (1992). The next theorem gives a simplification of the BSDE time discretization expressions.

Theorem 4.

Under Assumption 3 (1), the solution of the q-stage Runge–Kutta scheme satisfies

$${\left\{(Y_{i,j}^{\pi },Z_{i,j}^{\pi })\right\}}_{j=2}^{q+1}\in {\mathcal{F}}_{t_i}$$

(92)

for $0\le i<n$. Consequently, we can write

$$\begin{array}{cc}\hfill Z_{i,j}^{\pi }& ={\mathbf{E}}_{t_i}^{}\left[H_{t_{i,j},{\gamma }_j{\Delta }_i}^{{\phi }_j}\left(Y_{t_{i+1}}^{\pi }+{\Delta }_i{\beta }_{j,1}f(t_{i+1},X_{t_{i+1}}^{\pi },Y_{t_{i+1}}^{\pi },Z_{t_{i+1}}^{\pi })\right)\right]\hfill \\ \hfill Y_{i,j}^{\pi }& ={\mathbf{E}}_{t_i}^{}\left[Y_{t_{i+1}}^{\pi }+{\Delta }_i{\alpha }_{j,1}f(t_{i+1},X_{t_{i+1}}^{\pi },Y_{t_{i+1}}^{\pi },Z_{t_{i+1}}^{\pi })\right]\hfill \end{array}$$

(93)

$$\begin{array}{cc}\hfill & +{\Delta }_i\sum _{k=2}^j{\alpha }_{jk}f(t_{i,k},X_{t_i}^{\pi },Y_{i,k}^{\pi },Z_{i,k}^{\pi })\hfill \end{array}$$

(94)

for $0\le i<n$ and $1<j\le q+1$ where ${\phi }_{q+1}={\phi }_1$, ${\beta }_{q+1,1}={\beta }_1$ and ${\alpha }_{q+1,k}={\alpha }_k$.

Proof. 

Clearly $(Y_{i,q+1}^{\pi },Z_{i,q+1}^{\pi })=(Y_{t_i}^{\pi },Z_{t_i}^{\pi })\in {\mathcal{F}}_{t_i}$ from Equations (75) and (76). For $1<j\le q$ and $0\le i<n$, we have

$$\begin{array}{cc}\hfill Y_{i,j}^{\pi }& ={\mathbf{E}}_{}^{}\left[Y_{t_{i+1}}^{\pi }+{\Delta }_i\sum _{k=1}^j{\alpha }_{jk}f(t_{i,k},X_{t_{i,k}}^{\pi },Y_{i,k}^{\pi },Z_{i,k}^{\pi })X_{t_{i,j}}^{\pi }\right](\mathrm{starting}\mathrm{from}\mathrm{Equation}(78))\hfill \\ \hfill & ={\mathbf{E}}_{}^{}\left[Y_{t_{i+1}}^{\pi }+{\Delta }_i\sum _{k=1}^j{\alpha }_{jk}f(t_{i,k},X_{t_{i,k}}^{\pi },Y_{i,k}^{\pi },Z_{i,k}^{\pi })X_{t_i}^{\pi }\right]\hfill \\ \hfill & \hfill (\mathrm{by}\mathrm{Assumption}3\mathrm{since}{t}_{i,j}\in [t_i,t_{i+1}))\\ \hfill & ={\mathbf{E}}_{t_i}^{}\left[Y_{t_{i+1}}^{\pi }+{\Delta }_i\sum _{k=1}^j{\alpha }_{jk}f(t_{i,k},X_{t_{i,k}}^{\pi },Y_{i,k}^{\pi },Z_{i,k}^{\pi })\right]\hfill \end{array}$$

so that $Y_{i,j}^{\pi }\in {\mathcal{F}}_{t_i}$. Similar arguments also show that $Z_{i,j}^{\pi }\in {\mathcal{F}}_{t_i}$ starting from Equation (77).

Since ${\left\{(Y_{i,j}^{\pi },Z_{i,j}^{\pi })\right\}}_{j=2}^{q+1}\in {\mathcal{F}}_{t_i}$, we naturally get Equation (94) from Equations (78) and (76). In addition, knowing that

$${\mathbf{E}}_{t_i}^{}\left[H_{t_{i,j},({\gamma }_i-{\gamma }_k){\Delta }_i}^{{\phi }_j}\right]=0,1<k<j$$

(95)

leads to Equation (93) from Equations (77) and (75).    □

As a consequence of Assumption 3, if the $q-$stage Runge–Kutta scheme and the forward SDE time discretization are of order $m>0$ then error of the FBSDE numerical solution defined as ${\mathcal{E}}_{X,\pi }+{\mathcal{E}}_{\pi }$ is of order m. We must hence choose the Runge–Kutta scheme and the SDE scheme accordingly.

3.3. Fourier Representation

Following Theorem 4, the intermediate solutions ${\left\{(u_{i,j},{\dot{u}}_{i,j})\right\}}_{j=2}^{q+1}$ at mesh time $t_i$, $0\le i<n$, are given by

$$\begin{array}{cc}\hfill {\dot{u}}_{i,j}\left(x\right)& ={\mathbf{E}}_{}^{}\left[H_{t_{i,j},{\gamma }_j{\Delta }_i}^{{\phi }_j}{\tilde{u}}_{i+1}(X_{t_{i+1}}^{\pi },{\beta }_{j,1})|X_{t_i}^{\pi }=x\right]\hfill \end{array}$$

(96)

$$\begin{array}{cc}\hfill u_{i,j}\left(x\right)& ={\mathbf{E}}_{}^{}\left[{\tilde{u}}_{i+1}(X_{t_{i+1}}^{\pi },{\alpha }_{j,1})|X_{t_i}^{\pi }=x\right]+{\Delta }_i\sum _{k=2}^j{\alpha }_{jk}f(t_{i,k},x,u_{i,k}\left(x\right),{\dot{u}}_{i,k}\left(x\right))\hfill \end{array}$$

(97)

for $1<j\le q+1$ with ${\phi }_{q+1}={\phi }_1$, ${\beta }_{q+1,1}={\beta }_1$ and ${\alpha }_{q+1,k}={\alpha }_k$. The approximate solution $u_i$ and approximate gradient ${\dot{u}}_i$ at mesh time $t_i$, $0\le i<n$, are then

$$\begin{array}{c}\hfill u_i\left(x\right)=u_{i,q+1}\left(x\right)\end{array}$$

(98)

$$\begin{array}{c}\hfill {\dot{u}}_i\left(x\right)={\dot{u}}_{i,q+1}\left(x\right)\end{array}$$

(99)

with

$${\tilde{u}}_{i+1}(x,\alpha )=u_{i+1}\left(x\right)+{\Delta }_i\alpha f(t_{i+1},x,u_{i+1}\left(x\right),{\dot{u}}_{i+1}\left(x\right))$$

(100)

and

$$\begin{array}{cc}\hfill u_n\left(x\right)& =g\left(x\right)\hfill \end{array}$$

(101)

$$\begin{array}{cc}\hfill {\dot{u}}_n\left(x\right)& ={\sigma }^{*}(T,x)\nabla g\left(x\right).\hfill \end{array}$$

(102)

In this setting, we have that

$$u_{i,j}\left(x\right)={\mathbf{E}}_{}^{}\left[{\tilde{u}}_{i+1}(X_{t_{i+1}}^{\pi },{\alpha }_{j,1})|X_{t_i}^{\pi }=x\right]+{\Delta }_i\sum _{k=2}^j{\alpha }_{jk}f(t_{i,k},x,u_{i,k}\left(x\right),{\dot{u}}_{i,k}\left(x\right))$$

(103)

Note that

$$\begin{array}{cc}\hfill & {\mathbf{E}}_{}^{}\left[{\tilde{u}}_{i+1}(X_{t_{i+1}}^{\pi },{\alpha }_{j,1})|X_{t_i}^{\pi }=x\right]={\mathbf{E}}_{t_i}^x\left[\frac{1}{{\left(2\pi \right)}^d}{\int }_{{\mathbb{R}}^d}{e}^{\mathbf{i}{\nu }^{*}{X}_{t_{i+1}}^{\pi }}\mathfrak{F}\left[{\tilde{u}}_{i+1}(.,{\alpha }_{j,1})\right]\left(\nu \right)d\nu \right]\hfill \end{array}$$

$$\begin{array}{cc}\hfill & =\frac{1}{{\left(2\pi \right)}^d}{\int }_{{\mathbb{R}}^d}{\mathbf{E}}_{t_i}^x\left[e^{\mathbf{i}{\nu }^{*}{X}_{t_{i+1}}^{\pi }}\right]\mathfrak{F}\left[{\tilde{u}}_{i+1}(.,{\alpha }_{j,1})\right]\left(\nu \right)d\nu (\mathrm{using}\mathrm{Fubini}’\mathrm{s}\mathrm{theorem})\hfill \end{array}$$

$$\begin{array}{cc}\hfill & =\frac{1}{{\left(2\pi \right)}^d}{\int }_{{\mathbb{R}}^d}{e}^{\mathbf{i}{\nu }^{*}x}{\varphi }_i(\nu ,x)\mathfrak{F}\left[{\tilde{u}}_{i+1}(.,{\alpha }_{j,1})\right]\left(\nu \right)d\nu .\hfill \end{array}$$

(104)

Therefore, by (103) and (104), we have

$$\begin{array}{cc}\hfill u_{i,j}\left(x\right)& ={\mathfrak{F}}^{-1}\left[\mathfrak{F}\left[{\tilde{u}}_{i+1}(.,{\alpha }_{j,1})\right]\left(\nu \right){\varphi }_i(\nu ,x)\right]\left(x\right)\hfill \\ \hfill & +{\Delta }_i\sum _{k=2}^j{\alpha }_{jk}f(t_{i,k},x,u_{i,k}\left(x\right),{\dot{u}}_{i,k}\left(x\right))\hfill \end{array}$$

(105)

whenever ${\tilde{u}}_{i+1}(.,\alpha )$ is Lebesgue integrable.

As to the intermediate solutions ${\dot{u}}_{i,j}$, $0\le i<n$ and $1<j\le q+1$, we have

$$\begin{array}{cc}\hfill {\dot{u}}_{i,j}\left(x\right)& ={\mathbf{E}}_{}^{}\left[H_{t_{i,j},{\gamma }_j{\Delta }_i}^{{\phi }_j}{\tilde{u}}_{i+1}(X_{t_{i+1}}^{\pi },{\beta }_{j,1})|X_{t_i}^{\pi }=x\right]\hfill \\ \hfill & ={\mathbf{E}}_{t_i}^x\left[H_{t_{i,j},{\gamma }_j{\Delta }_i}^{{\phi }_j}\frac{1}{{\left(2\pi \right)}^d}{\int }_{{\mathbb{R}}^d}{e}^{\mathbf{i}{\nu }^{*}{X}_{t_{i+1}}^{\pi }}\mathfrak{F}\left[{\tilde{u}}_{i+1}(.,{\beta }_{j,1})\right]\left(\nu \right)d\nu \right]\hfill \\ \hfill & =\frac{1}{{\left(2\pi \right)}^d}{\int }_{{\mathbb{R}}^d}{\mathbf{E}}_{t_i}^x\left[H_{t_{i,j},{\gamma }_j{\Delta }_i}^{{\phi }_j}{e}^{\mathbf{i}{\nu }^{*}{X}_{t_{i+1}}^{\pi }}\right]\mathfrak{F}\left[{\tilde{u}}_{i+1}(.,{\beta }_{j,1})\right]\left(\nu \right)d\nu (\mathrm{using}\mathrm{Fubini}’\mathrm{s}\mathrm{theorem})\hfill \\ \hfill & =\frac{1}{{\left(2\pi \right)}^d}{\int }_{{\mathbb{R}}^d}{e}^{\mathbf{i}{\nu }^{*}x}{\mathsf{\Phi }}_{i,j}(\nu ,x)\mathfrak{F}\left[{\tilde{u}}_{i+1}(.,{\beta }_{j,1})\right]\left(\nu \right)d\nu \hfill \end{array}$$

$$\begin{array}{c}={\mathfrak{F}}^{-1}\left[\mathfrak{F}\left[{\tilde{u}}_{i+1}(.,{\alpha }_{j,1})\right]\left(\nu \right){\mathsf{\Phi }}_{i,j}(\nu ,x)\right]\left(x\right)\hfill \end{array}$$

(106)

for an integrable function ${\tilde{u}}_{i+1}(.,\alpha )$.

Even if the expressions in Equations (105) and (106) appear too general, they are implementable with the Fourier interpolation method for $d=1$ in various particular cases. One can retrieve the characteristics ${\varphi }_i$ and ${\mathsf{\Phi }}_{i,j}$ and also perform the corrections due to the transform of Equation (44) for many SDE time discretizations. The following lemma helps in retrieving the conditional characteristics.

Lemma 2.

The conditional characteristics ${\mathsf{\Phi }}_{i,j}$ are

$${\mathsf{\Phi }}_{i,j}(\nu ,x)={\mathbf{E}}_{t_i}^x\left[{\mathbf{H}}_{i,j}^{*}\mathbf{i}\nu e^{\mathbf{i}{\nu }^{*}\left(X_{t_{i+1}}^{\pi }-X_{t_i}^{\pi }\right)}\right]$$

(107)

with

$${\mathbf{H}}_{i,j}=\frac{1}{{\gamma }_j{\Delta }_i}{\int }_{t_{i,j}}^{t_{i+1}}{D}_s{X}_{t_{i+1}}^{\pi }{\phi }_j\left(\frac{s-t_{i,j}}{{\gamma }_j{\Delta }_i}\right)ds$$

(108)

where $D_s{X}_{t_{i+1}}^{\pi }$ is the Malliavin derivative of $X_{t_{i+1}}^{\pi }$ given $X_{t_i}^{\pi }=x$.

Proof. 

The lemma is proved by applying the duality formula and the chain rule successively to Equation (89).    □

3.3.1. Half-Order Itô-Taylor Schemes

The Euler scheme constitutes the main example of half-order Itô-Taylor scheme with step

$$X_{t_{i+1}}^{\pi }=X_{t_i}^{\pi }+a(t_i,X_{t_i}^{\pi }){\Delta }_i+\sigma (t_i,X_{t_i}^{\pi })\Delta W_i.$$

In addition, we have that $D_s{\Delta }_i={\mathbf{0}}_{d\times 1}$ and $D_s\Delta W_i={\mathbf{I}}_{d\times d}$ for $s\in (t_i,t_{i+1})$ where $\mathbf{0}$ and $\mathbf{I}$ are the zero matrix and the identity matrix, respectively. Hence,

$$D_s{X}_{t_{i+1}}^{\pi }=\sigma (t_i,X_{t_i}^{\pi })$$

so we get, from Equation (107), that

$$\begin{array}{cc}\hfill {\mathsf{\Phi }}_{i,j}(\nu ,x)& ={\sigma }^{*}(t_i,x)\mathbf{i}\nu {\mathbf{E}}_{t_i}^x\left[e^{\mathbf{i}{\nu }^{*}\left(X_{t_{i+1}}^{\pi }-X_{t_i}^{\pi }\right)}\right](\mathrm{since}{\phi }_j\in {\mathcal{B}}^0),\hfill \\ \hfill & ={\sigma }^{*}(t_i,x)\mathbf{i}\nu {\varphi }_i(\nu ,x).\hfill \end{array}$$

(109)

The conditional characteristic function is explicitly given by

$${\varphi }_i(\nu ,x)=exp\left\{{\Delta }_i\left(\mathbf{i}{\nu }^{*}a(t_i,x)-\frac{1}{2}{\nu }^{*}{\sigma }^2(t_i,x)\nu \right)\right\}$$

(110)

since the increment has a Gaussian distribution.

Equations (105) and (106) along with the characteristics of Equations (110) and (109) define the Fourier method under half-order Itô-Taylor schemes and the method is implementable in one dimension ($d=1$) with the procedure given in Section 2.3. The following theorem generalizes the result of Theorem 1 to Runge–Kutta schemes under half-order Itô-Taylor schemes.

Theorem 5.

Let ${\tilde{u}}_{i+1}^{\alpha ,\beta }(.,y)$ be the alternative transform defined in Equation (44) of the approximate solution ${\tilde{u}}_{i+1}(.,y)$. Then the intermediate solutions $u_{i,j}$ and ${\dot{u}}_{i,j}$ in Equations (105) and (106) satisfy

$$\begin{array}{cc}\hfill u_{i,j}\left(x\right)& ={\mathfrak{F}}^{-1}\left[\mathfrak{F}\left[{\tilde{u}}_{i+1}^{\alpha ,\beta }(.,{\alpha }_{j,1})\right]\left(\nu \right){\varphi }_i(\nu ,x)\right]\left(x\right)\hfill \\ \hfill & -\alpha [{(x+{\Delta }_{i}a(t_i,x))}^2+{\Delta }_i{\sigma }^2(t_i,x)]-\beta (x+{\Delta }_{i}a(t_i,x))\hfill \\ \hfill & +{\Delta }_i\sum _{k=2}^j{\alpha }_{jk}f(t_{i,k},x,u_{i,k}\left(x\right),{\dot{u}}_{i,k}\left(x\right))\hfill \end{array}$$

(111)

$$\begin{array}{cc}\hfill {\dot{u}}_{i,j}\left(x\right)& =\sigma (t_i,x){\mathfrak{F}}^{-1}\left[\mathfrak{F}\left[{\tilde{u}}_{i+1}^{\alpha ,\beta }(.,{\beta }_{j,1})\right]\left(\nu \right)\mathbf{i}\nu {\varphi }_i(\nu ,x)\right]\left(x\right)\hfill \\ \hfill & -\sigma (t_i,x)[2\alpha (x+{\Delta }_{i}a(t_i,x))+\beta ].\hfill \end{array}$$

(112)

under a half-order Itô-Taylor scheme when $d=1$.

3.3.2. First-Order Itô-Taylor Schemes

Consider the first-order scheme

$$X_{t_{i+1}}^{\pi }=X_{t_i}^{\pi }+a(t_i,X_{t_i}^{\pi }){\Delta }_i+\sigma (t_i,X_{t_i}^{\pi })\Delta W_i+{\sigma }^2(t_i,X_{t_i}^{\pi }){\int }_{t_i}^{t_{i+1}}{\int }_{t_i}^{t}d{W}_{u}d{W}_t.$$

Then knowing that

$$D_s{\int }_{t_i}^{t_{i+1}}{\int }_{t_i}^{t}d{W}_{u}d{W}_t=\mathbf{D}(\Delta W_i),$$

(113)

using the fundamental theorem of calculus where $\mathbf{D}\left(x\right)$ is the diagonal matrix composed with the elements of x, for $s\in (t_i,t_{i+1})$, the Malliavin derivative of the discretized forward process is given by

$$D_s{X}_{t_{i+1}}^{\pi }=\sigma (t_i,x)+{\sigma }^2(t_i,x)\mathbf{D}(\Delta W_i).$$

(114)

Equation (107) leads to

$$\begin{array}{ll}{\mathsf{\Phi }}_{i,j}(\nu ,x)& ={\sigma }^{*}(t_i,x)\mathbf{i}\nu {\mathbf{E}}_{t_i}^x\left[e^{\mathbf{i}{\nu }^{*}\left(X_{t_{i+1}}^{\pi }-X_{t_i}^{\pi }\right)}\right]\hfill \\ \hfill & \hfill +{\mathbf{E}}_{t_i}^x\left[\mathbf{D}(\Delta W_i){\sigma }^2(t_i,x)\mathbf{i}\nu e^{\mathbf{i}{\nu }^{*}\left(X_{t_{i+1}}^{\pi }-X_{t_i}^{\pi }\right)}\right](\mathrm{since}{\phi }_j\in {\mathcal{B}}^0)\\ \hfill & ={\sigma }^{*}(t_i,x)\mathbf{i}\nu {\varphi }_i(\nu ,x)+{\mathbf{E}}_{t_i}^x\left[\mathbf{D}\left({\sigma }^2(t_i,x)\mathbf{i}\nu \right)\Delta W_i{e}^{\mathbf{i}{\nu }^{*}\left(X_{t_{i+1}}^{\pi }-X_{t_i}^{\pi }\right)}\right]\hfill \end{array}$$

$$\begin{array}{cc}\hfill & ={(\mathbf{I}-\mathbf{D}\left({\sigma }^2(t_i,x)\mathbf{i}\nu \right){\Delta }_i)}^{-1}{\sigma }^{*}(t_i,x)\mathbf{i}\nu {\varphi }_i(\nu ,x)\hfill \end{array}$$

(115)

since

$$\begin{array}{cc}\hfill {\mathbf{E}}_{t_i}^x\left[\Delta W_i{e}^{\mathbf{i}{\nu }^{*}\left(X_{t_{i+1}}^{\pi }-X_{t_i}^{\pi }\right)}\right]& ={\sigma }^{*}(t_i,x){\Delta }_i{\mathbf{E}}_{t_i}^x\left[\mathbf{i}\nu e^{\mathbf{i}{\nu }^{*}\left(X_{t_{i+1}}^{\pi }-X_{t_i}^{\pi }\right)}\right]\hfill \\ & +\mathbf{D}\left({\sigma }^2(t_i,x)\mathbf{i}\nu \right){\Delta }_i{\mathbf{E}}_{t_i}^x\left[\Delta W_i{e}^{\mathbf{i}{\nu }^{*}\left(X_{t_{i+1}}^{\pi }-X_{t_i}^{\pi }\right)}\right]\hfill \end{array}$$

using the duality formula, so that

$$\begin{array}{cc}\hfill {\mathbf{E}}_{t_i}^x\left[\Delta W_i{e}^{\mathbf{i}{\nu }^{*}\left(X_{t_{i+1}}^{\pi }-X_{t_i}^{\pi }\right)}\right]& ={\zeta }_i(\nu ,x){\Delta }_i{\sigma }^{*}(t_i,x){\mathbf{E}}_{t_i}^x\left[\mathbf{i}\nu e^{\mathbf{i}{\nu }^{*}\left(X_{t_{i+1}}^{\pi }-X_{t_i}^{\pi }\right)}\right]\hfill \\ \hfill & ={\zeta }_i(\nu ,x){\Delta }_i{\sigma }^{*}(t_i,x)\mathbf{i}\nu {\varphi }_i(\nu ,x)\hfill \end{array}$$

with

$${\zeta }_i(\nu ,x)={(\mathbf{I}-\mathbf{D}\left({\sigma }^2(t_i,x)\mathbf{i}\nu \right){\Delta }_i)}^{-1}.$$

(116)

As to the conditional characteristic ${\varphi }_i$, it can be easily derived as

$${\varphi }_i(\nu ,x)=det{\left({\zeta }_i(\nu ,x)\right)}^{\frac{1}{2}}exp\left(\frac{1}{2}\mathbf{i}{\nu }^{*}{\zeta }_i^{-1}(\nu ,x){\mathbf{1}}_{d\times 1}+\mathbf{i}{\nu }^{*}{\kappa }_i\left(x\right)\right)$$

(117)

where

$${\kappa }_i\left(x\right)=a(t_i,x){\Delta }_i-\frac{1}{2}({\sigma }^2(t_i,x){\Delta }_i+1){\mathbf{1}}_{d\times 1}$$

knowing that $X_{t_{i+1}}^{\pi }-X_{t_i}^{\pi }$, given $X_{t_i}^{\pi }=x$, is an affine function of a multivariate non-central ${\chi }^2$ random variable with 1 degree of freedom and non-centrality parameters 1.

Equations (105) and (106) along with the expressions in Equations (115) and (117) characterize the method under first order discretizations on the forward process. The procedure introduced in Section 2.3 allows us to do the necessary computations given the characteristics ${\varphi }_i$ and ${\mathsf{\Phi }}_{i,j}$ and using the following theorem.

Theorem 6.

Let ${\tilde{u}}_{i+1}^{\alpha ,\beta }(.,y)$ be the alternative transform defined in Equation (44) of the approximate solution ${\tilde{u}}_{i+1}(.,y)$. Then the intermediate solutions $u_{i,j}$ and ${\dot{u}}_{i,j}$ in Equations (105) and (106) satisfy

$$\begin{array}{cc}\hfill u_{i,j}\left(x\right)& ={\mathfrak{F}}^{-1}\left[\mathfrak{F}\left[{\tilde{u}}_{i+1}^{\alpha ,\beta }(.,{\alpha }_{j,1})\right]\left(\nu \right){\varphi }_i(\nu ,x)\right]\left(x\right)\hfill \\ \hfill & -\alpha \left[{(x+{\Delta }_{i}a(t_i,x))}^2+{\Delta }_i{\sigma }^2(t_i,x)+\frac{1}{2}{\Delta }_i^2{\sigma }^4(t_i,x)\right]\hfill \\ \hfill & -\beta (x+{\Delta }_{i}a(t_i,x))+{\Delta }_i\sum _{k=2}^j{\alpha }_{jk}f(u_{i,k}\left(x\right),{\dot{u}}_{i,k}\left(x\right))\hfill \end{array}$$

(118)

$$\begin{array}{cc}\hfill {\dot{u}}_{i,j}\left(x\right)& =\sigma (t_i,x){\mathfrak{F}}^{-1}\left[\mathfrak{F}\left[{\tilde{u}}_{i+1}^{\alpha ,\beta }(.,{\beta }_{j,1})\right]\left(\nu \right)\mathbf{i}\nu {\zeta }_i(\nu ,x){\varphi }_i(\nu ,x)\right]\left(x\right)\hfill \\ \hfill & -\sigma (t_i,x)\left[2\alpha \left(x+{\Delta }_{i}a(t_i,x)+{\Delta }_i{\sigma }^2(t_i,x)\right)+\beta \right]\hfill \end{array}$$

(119)

under a first-order Itô-Taylor scheme when $d=1$.

Proof. 

By the definition of the alternative transform, we must have that

$$\begin{array}{cc}\hfill u_{i,j}\left(x\right)& ={\mathfrak{F}}^{-1}\left[\mathfrak{F}\left[{\tilde{u}}_{i+1}^{\alpha ,\beta }(.,{\alpha }_{j,1})\right]\left(\nu \right){\varphi }_i(\nu ,x)\right]\left(x\right)-{\mathbf{E}}_{t_i}^x\left[\alpha {\left(X_{t_{i+1}}^{\pi }\right)}^2+\beta X_{t_{i+1}}^{\pi }\right]\hfill \\ \hfill & +{\Delta }_i{\mathbf{1}}_{\{j>2\}}\sum _{k=2}^{j-1}{\alpha }_{jk}f(u_{i,k}\left(x\right),{\dot{u}}_{i,k}\left(x\right)).\hfill \end{array}$$

(120)

Notice that

$${\mathbf{E}}_{t_i}^x\left[X_{t_{i+1}}^{\pi }\right]=x+{\Delta }_{i}a(t_i,x)$$

(121)

and

$$\begin{array}{cc}\hfill & {\mathbf{E}}_{t_i}^x\left[{\left(X_{t_{i+1}}^{\pi }\right)}^2\right]={\mathbf{E}}_{t_i}^x{\left[X_{t_{i+1}}^{\pi }\right]}^2+{\mathbf{Var}}_{t_i}^x\left[X_{t_{i+1}}^{\pi }\right]\hfill \\ \hfill & ={(x+{\Delta }_{i}a(t_i,x))}^2+{\mathbf{E}}_{t_i}^x\left[{\left(\sigma (t_i,x)\Delta W_i+{\sigma }^2(t_i,x){\int }_{t_i}^{t_{i+1}}{\int }_{t_i}^{t}d{W}_{u}d{W}_t\right)}^2\right]\hfill \\ \hfill & ={(x+{\Delta }_i\sigma (t_i,x))}^2+{\Delta }_i{\sigma }^2(t_i,x)+\frac{1}{2}{\Delta }_i^2{\sigma }^4(t_i,x).\hfill \end{array}$$

(122)

Equations (120)–(122) lead to the expression for $u_{i,j}$ in Equation (118).

The definition of the alternative transform also requires

$$\begin{array}{cc}\hfill {\dot{u}}_{i,j}\left(x\right)& =\sigma (t_i,x){\mathfrak{F}}^{-1}\left[\mathfrak{F}\left[{\tilde{u}}_{i+1}^{\alpha ,\beta }(.,{\beta }_{j,1})\right]\left(\nu \right)\mathbf{i}\nu {\zeta }_i(\nu ,x){\varphi }_i(\nu ,x)\right]\left(x\right)\hfill \\ \hfill & -{\mathbf{E}}_{t_i}^x\left[H_{t_{i,j},{\gamma }_j{\Delta }_i}^{{\phi }_j}(\alpha {\left(X_{t_{i+1}}^{\pi }\right)}^2+\beta X_{t_{i+1}}^{\pi })\right]\hfill \\ \hfill & =\sigma (t_i,x){\mathfrak{F}}^{-1}\left[\mathfrak{F}\left[{\tilde{u}}_{i+1}^{\alpha ,\beta }(.,{\beta }_{j,1})\right]\left(\nu \right)\mathbf{i}\nu {\zeta }_i(\nu ,x){\varphi }_i(\nu ,x)\right]\left(x\right)\hfill \\ \hfill & -\sigma (t_i,x){\mathbf{E}}_{t_i}^x\left[2\alpha \left(X_{t_{i+1}}^{\pi }\right)+\beta \right]\hfill \\ \hfill & -{\sigma }^2(t_i,x){\mathbf{E}}_{t_i}^x\left[\Delta W_i\left(2\alpha \left(X_{t_{i+1}}^{\pi }\right)+\beta \right)\right](\mathrm{using}\mathrm{the}\mathrm{duality}\mathrm{formula})\hfill \\ \hfill & =\sigma (t_i,x){\mathfrak{F}}^{-1}\left[\mathfrak{F}\left[{\tilde{u}}_{i+1}^{\alpha ,\beta }(.,{\beta }_{j,1})\right]\left(\nu \right)\mathbf{i}\nu {\zeta }_i(\nu ,x){\varphi }_i(\nu ,x)\right]\left(x\right)\hfill \end{array}$$

$$-\sigma (t_i,x)\left[2\alpha \left(x+{\Delta }_{i}a(t_i,x)+{\Delta }_i{\sigma }^2(t_i,x)\right)+\beta \right]$$

(123)

using the duality formula once again.    □

The implementation of higher order time discretization for FBSDEs on the alternative grid is described in the following algorithm. Algorithm 2 produces the numerical intermediate solutions ${\left\{u_{i,j,k}\right\}}_{k=0}^{N_{i}N}$, ${\left\{{\tilde{u}}_{i,j,k}\right\}}_{k=0}^{N_{i}N}$ and ${\left\{{\dot{u}}_{i,j,k}\right\}}_{k=0}^{N_{i}N}$ at time step $t_i$, $0\le i<n$ and stage j, $1\le j\le q+1$ for the approximate solution $u_i$, the intermediate solution ${\tilde{u}}_i$ and the approximate gradient ${\dot{u}}_i$, respectively, $i=0,1,\dots ,n-1$.

Algorithm 2 Fourier Interpolation Method on Alternative Grid for q-stage Runge–Kutta schemes

Discretize the restricted real space $[x_0-\frac{N_{n}l}{2},x_0+\frac{N_{n}l}{2}]$ and the restricted Fourier space $[-\frac{L}{2},\frac{L}{2}]$ with $N_{n}N$ space steps so to have the real space nodes ${\left\{x_{nk}\right\}}_{k=0}^{N_{n}N}$ and ${\left\{{\nu }_{nk}\right\}}_{k=0}^{N_{n}N}$ Value $u_n\left(x_{nk}\right)=g\left(x_{nk}\right)$ For any i from $n-1$ to 0 (a) For any j, $1<j\le q+1$ Compute $\alpha$ and $\beta$ defining the transform of Equation (44), such that $${\theta }_{i+1}={\tilde{u}}_{i+1}^{\alpha ,\beta }(.,{\alpha }_{j,1})$$ (124) and ${\theta }_{i+1}$ satisfies the boundary conditions of Equations (32) and (33). Compute ${\theta }_i\left(x_{ik}\right)$ through Equation (40) for $k=0,1,\dots ,N_{i}N$ with $$\psi (\nu ,x)={\varphi }_i(\nu ,x)$$ (125) and retrieve the values ${\tilde{u}}_{i,j,k}$ with the appropriate correction. Compute $\alpha$ and $\beta$ defining the transform of Equation (44), such that $${\theta }_{i+1}={\tilde{u}}_{i+1}^{\alpha ,\beta }(.,{\beta }_{j,1})$$ (126) and ${\theta }_{i+1}$ satisfies the boundary conditions of Equations (32) and (33). Compute ${\theta }_i\left(x_{ik}\right)$ through Equation (40) for $k=0,1,\dots ,N_{i}N$ with $$\psi (\nu ,x)={\mathsf{\Phi }}_{i,j}(\nu ,x)$$ (127) and retrieve the values ${\dot{u}}_{i,j,k}$ with the appropriate correction. Compute the values $u_{i,j,k}$ as $$u_{i,j,k}={\tilde{u}}_{i,j,k}+{\Delta }_i\sum _{s=2}^j{\alpha }_{js}f(t_{i,s},x_{i,k},u_{i,s,k},{\dot{u}}_{i,s,k})$$ (128) for $k=0,1,\dots ,N_{i}N$ through Equation (13). Update the real space grid with Equation (27) and the Fourier space grid by discretizing the interval $[-\frac{L}{2},\frac{L}{2}]$ with $N_{i}N$ space steps so to have the real space nodes ${\left\{x_{ik}\right\}}_{k=0}^{N_{i}N}$ and ${\left\{{\nu }_{ik}\right\}}_{k=0}^{N_{i}N}$. (b) Set $u_{i,k}=u_{i,q+1,k}$ and ${\dot{u}}_{i,k}={\dot{u}}_{i,q+1,k}$

Algorithm 2, as compared to Algorithm 1, contains intermediate solutions between time discretization points, resulting from the the intermediate stages of the Runge–Kutta method of Chassagneux and Crisan (2014). These intermediate solutions are computed using the Fourier interpolation method on the alternative grid presented in Section 2. We close this section by examining the spatial discretization error.

3.4. Spatial Discretization Error Analysis

We denote by ${\left\{u_{i,j,k}\right\}}_{k=0}^{N_{i}N}$ and ${\left\{{\dot{u}}_{i,j,k}\right\}}_{k=0}^{N_{i}N}$ the intermediate numerical solutions obtained at time step $t_i$, $i=0,1,\dots ,n-1$ and stage j, $1<j\le q+1$, from the Fourier interpolation method on the alternative grid when using a $q-$stage Runge–Kutta scheme. In addition, ${\left\{{\mathbf{u}}_{i,j,k}\right\}}_{k=0}^{N_{i}N}$ and ${\left\{{\dot{\mathbf{u}}}_{i,j,k}\right\}}_{k=0}^{N_{i}N}$ are the intermediate numerical solutions obtained at the intermediate stage j, $1<j\le q+1$, of time step $t_i$ given the exact solutions $u_{i+1}$ and ${\dot{u}}_{i+1}$ at $t_{i+1}$. We have from the notation previously used that the numerical solutions at $t_i$ write

$$\begin{array}{cc}\hfill u_{i,k}& =u_{i,q+1,k}\hfill \end{array}$$

(129)

$$\begin{array}{cc}\hfill {\dot{u}}_{i,k}& ={\dot{u}}_{i,q+1,k}\hfill \end{array}$$

(130)

and are computed from the intermediate solutions${\left\{{\tilde{u}}_{i,k}\right\}}_{k=0}^{N_{i}N}$, $0<i\le n$ where ${\tilde{u}}_{n,k}={\tilde{u}}_n\left(x_{n,k}\right)$. When the exact solutions $u_{i+1}$ and ${\dot{u}}_{i+1}$ are known at $t_{i+1}$, we also write

$$\begin{array}{cc}\hfill {\mathbf{u}}_{i,k}& ={\mathbf{u}}_{i,q+1,k}\hfill \end{array}$$

(131)

$$\begin{array}{cc}\hfill {\dot{\mathbf{u}}}_{i,k}& ={\dot{\mathbf{u}}}_{i,q+1,k}.\hfill \end{array}$$

(132)

The local (space) discretization error has the form

$$E_{ik}:=\left|u_i\left(x_k\right)-{\mathbf{u}}_{i,k}\right|+\left|{\dot{u}}_i\left(x_k\right)-{\dot{\mathbf{u}}}_{i,k}\right|$$

(133)

for $i=0,1,\dots ,n-1$ and $k=0,1,\dots ,N_{i}N$. The next theorem gives a description of the local (space) discretization error bound.

Theorem 7.

Suppose that the driver $f\in {\mathcal{C}}^{1,2}([0,T]\times {\mathbb{R}}^2)$ and the terminal condition $g\in {\mathcal{C}}^2\left(\mathbb{R}\right)$ and Assumptions 1 and 3 are satisfied, then the Fourier interpolation method yields a local space discretization error of the form

$$\underset{i,k}{sup}{E}_{ik}=\mathcal{O}\left(\Delta x\right)+\mathcal{O}\left(e^{-K{\Delta }_i^{-s_0}{l}^{q_0}}\right)$$

(134)

for some constant $K>0$ on the alternative grid and under the trapezoidal quadrature rule for any explicit q-stage Runge–Kutta scheme.

Proof. 

We follow the steps in the proof of Theorem 2. The truncation error when computing the numerical solutions ${\dot{\mathbf{u}}}_{i,j,k}$ is

$$\begin{array}{cc}\hfill & {\displaystyle {\mathbf{E}}_{t_i}^{x_{ik}}\left[H_{t_{i,j},{\gamma }_j{\Delta }_i}^{{\phi }_j}{\tilde{u}}_{i+1}(t_{i+1},X_{t_{i+1}}^{\pi };{\beta }_{j,1}){\mathbf{1}}_{\left|\Delta X_i^{\pi }\right|>\frac{l}{2}}\right]}\hfill \\ \hfill & <K{\mathbf{E}}_{t_i}^{x_{ik}}{\left[{\left(H_{t_{i,j},{\gamma }_j{\Delta }_i}^{{\phi }_j}\right)}^4\right]}^{\frac{1}{4}}{\mathbf{E}}_{t_i}^{x_{ik}}{\left[{\mathbf{1}}_{\left|\Delta X_i^{\pi }\right|>\frac{l}{2}}\right]}^{\frac{1}{4}}\hfill \\ \hfill & \hfill (\mathrm{using}\mathrm{Cauchy}-\mathrm{Schwarz}\mathrm{inequality}\mathrm{twice}\mathrm{since}{\tilde{u}}_{i+1}(t_{i+1},X_{t_{i+1}}^{\pi };.)\mathrm{is}\mathrm{sq}.\mathrm{int}.)\\ \hfill & <K{\Delta }_i^{-\frac{1}{2}}{\mathbf{E}}_{t_i}^{x_{ik}}{\left[{\mathbf{1}}_{\left|\Delta X_i^{\pi }\right|>\frac{l}{2}}\right]}^{\frac{1}{4}}(\mathrm{since}{H}_{t_{i,j},{\gamma }_j{\Delta }_i}^{{\phi }_j}\mathrm{is}\mathrm{of}\mathrm{Gaussian}\mathrm{distribution})\hfill \\ \hfill & \le K{\Delta }_i^{-\frac{1}{2}}{\left(\underset{s>0}{inf}{e}^{-s\frac{l}{2}}{\varphi }_i(-\mathbf{i}s)+\underset{s>0}{inf}{e}^{-s\frac{l}{2}}{\varphi }_i\left(\mathbf{i}s\right)\right)}^{\frac{1}{4}}(\mathrm{by}\mathrm{Chernoff}’\mathrm{s}\mathrm{inequality})\hfill \\ \hfill & <K{\Delta }_i^{-\frac{1}{2}}{e}^{-K_0{\Delta }_i^{-s_0}{l}^{q_0}}(\mathrm{by}\mathrm{Assumption}3)\hfill \\ \hfill & <K{e}^{-C{\Delta }_i^{-s_0}{l}^{q_0}}.\hfill \end{array}$$

The Fourier interpolation leads to a first-order space discretization error when computing the numerical solutions ${\dot{\mathbf{u}}}_{i,j,k}$ since the driver f and the terminal condition g are twice differentiable.

The same statements hold for the numerical solutions ${\mathbf{u}}_{i,2,k}$ using similar arguments. By recursion and using the Lipschitz property of the driver f, the statements hold for ${\mathbf{u}}_{i,j,k}$, $1<j\le q+1$. Since the time step $t_i$ and the space node $x_{ik}$ are arbitrary, the space truncation and discretization error bounds hold for any i and k.   □

Locally, the truncation error remains spectral. Nonetheless, it is of a unspecified index $q_0$ in this general setting where the conditional characteristic function ${\varphi }_i$ is itself unspecified. For higher order time discretizations, one can expect $q_0\le 2$ since the forward process increment $X_{t_{i+1}}^{\pi }-X_{t_i}^{\pi }$ has a heavy tail distribution. Indeed, the Gaussian distribution of forward process increments and the quadratic exponential form of their characteristic functions were the main reason for the spectral convergence of index 2 of the truncation error in Section 2.4. The space discretization error though is unchanged with first-order due to the second-order differentiability of the BSDE coefficients. However, the Fourier interpolation produces a space discretization error with a higher order when the driver f and the terminal function g have the required smoothness. In general, if $f\in {\mathcal{C}}_b^{m+1}$ and $g\in {\mathcal{C}}_b^{m+1}$, we can expect a space discretization error of order m which is the convergence order of the underlying Fourier interpolation.

We now turn to the global space discretization error defined as in Equation (63). The next theorem gives its error bound.

Theorem 8.

Suppose the conditions of Theorem 7 are satisfied. If the discretization is such that

$$\underset{i}{sup}\left\{\frac{C_0\Delta x}{\pi {\Delta }_i^{p_0}}\right\}\le 1$$

(135)

then the Fourier interpolation method is stable and yields a global discretization error $E_{l,\Delta x}$ of the form

$$E_{l,\Delta x}=\mathcal{O}(\Delta x)+\mathcal{O}\left(e^{-K{\left|\pi \right|}^{-s_0}{l}^{q_0}}\right)$$

(136)

where $K>0$ for any explicit q-stage Runge–Kutta scheme.

Proof. 

From the definition of the global space discretization error, we may write

$$\begin{array}{cc}\hfill e_{ik}& \le E_{n-i,k}+\left|{\mathbf{u}}_{n-i,k}-u_{n-i,k}\right|\hfill \end{array}$$

(137)

$$\begin{array}{cc}\hfill {\dot{e}}_{ik}& \le E_{n-i,k}+\left|{\dot{\mathbf{u}}}_{n-i,k}-{\dot{u}}_{n-i,k}\right|.\hfill \end{array}$$

(138)

Assume that the boundary values of the function ${\tilde{u}}_{i+1}$ and the sequence ${\tilde{u}}_{i+1,s}$ are matched on the alternative grid so that we do not have to treat the alternative transform. Under an explicit $q-$stage Runge–Kutta scheme, we have

$$\begin{array}{cc}\hfill & \left|{\dot{\mathbf{u}}}_{i,j,k}-{\dot{u}}_{i,j,k}\right|=\left|{\mathfrak{D}}^{-1}{\left[{\left\{{\mathsf{\Phi }}_{i,j}({\nu }_{i+1,m},x_{ik})\mathbb{D}{[{\tilde{u}}_{i+1}-{\tilde{u}}_{i+1,s}]}_m\right\}}_{m=0}^{N_{i+1}N-1}\right]}_{k+\frac{N}{2}}\right|\hfill \\ \hfill & \le \frac{{\sum }_{m=0}^{N_{i+1}N-1}\left|{\mathsf{\Phi }}_{i,j}({\nu }_{i+1,m},x_{ik})\right|}{N_{i+1}N}\underset{k}{sup}\left|{\tilde{u}}_{i+1}(x_{ik},{\beta }_{1,j})-{\tilde{u}}_{i+1,k}\right|\hfill \\ \hfill & \le \frac{\Delta x}{2\pi }\left({\int }_{\mathbb{R}}\left|{\mathsf{\Phi }}_{i,j}(\nu ,x_{i,k})\right|d\nu \right)\underset{k}{sup}\left|{\tilde{u}}_{i+1}(x_{ik},{\beta }_{1,j})-{\tilde{u}}_{i+1,k}\right|\hfill \\ \hfill & \le \frac{C_0\Delta x}{2\pi {\Delta }_i^{p_0}}\underset{k}{sup}\left|{\tilde{u}}_{i+1}(x_{ik},{\beta }_{1,j})-{\tilde{u}}_{i+1,k}\right|(\mathrm{using}\mathrm{Assumption}3)\hfill \\ \hfill & \le \frac{C_0\Delta x}{2\pi {\Delta }_i^{p_0}}(1+{\Delta }_{i}K)\underset{k}{sup}{e}_{n-i-1,k}+\frac{C_0\Delta x}{2\pi {\Delta }_i^{p_0}}{\Delta }_{i}K\underset{k}{sup}{\dot{e}}_{n-i-1,k}\hfill \\ \hfill & \hfill (\mathrm{since}f\mathrm{is}\mathrm{Lipschitz}\mathrm{and}{\beta }_{1,j}\mathrm{is}\mathrm{bounded})\end{array}$$

$$\begin{array}{ll}& \le \frac{C_0\Delta x}{2\pi {\Delta }_i^{p_0}}(1+{\Delta }_{i}K){\displaystyle \underset{k}{sup}}{e}_{n-i-1,k}+\frac{C_0\Delta x}{2\pi {\Delta }_i^{p_0}}(1+{\Delta }_{i}K){\displaystyle \underset{k}{sup}}{\dot{e}}_{n-i-1,k}.\hfill \end{array}$$

(139)

Similarly, we get

$$\begin{array}{cc}\hfill & \left|{\mathbf{u}}_{i,2,k}-u_{i,2,k}\right|\le \left|{\mathfrak{D}}^{-1}{\left[{\left\{{\varphi }_i({\nu }_{i+1,m},x_{ik})\mathbb{D}{[{\tilde{u}}_{i+1}-{\tilde{u}}_{i+1,s}]}_m\right\}}_{m=0}^{N_{i+1}N-1}\right]}_{k+\frac{N}{2}}\right|\hfill \\ \hfill & \le \frac{\Delta x}{2\pi }\left({\int }_{\mathbb{R}}\left|{\varphi }_i(\nu ,x_{i,k})\right|d\nu \right)\underset{k}{sup}\left|{\tilde{u}}_{i+1}(x_{ik},{\alpha }_{1,2})-{\tilde{u}}_{i+1,k}\right|\hfill \\ \hfill & \le \frac{C_0\Delta x}{2\pi {\Delta }_i^{p_0}}\underset{k}{sup}\left|{\tilde{u}}_{i+1}(x_{ik},{\alpha }_{1,2})-{\tilde{u}}_{i+1,k}\right|(\mathrm{using}\mathrm{Assumption}3)\hfill \\ \hfill & \le \frac{C_0\Delta x}{2\pi {\Delta }_i^{p_0}}(1+{\Delta }_{i}K)\underset{k}{sup}{e}_{n-i-1,k}+\frac{C_0\Delta x}{2\pi {\Delta }_i^{p_0}}(1+{\Delta }_{i}K)\underset{k}{sup}{\dot{e}}_{n-i-1,k}\hfill \end{array}$$

so that we get

$$\begin{array}{cc}\hfill \left|{\mathbf{u}}_{i,j,k}-u_{i,j,k}\right|& \le \frac{C_0\Delta x}{2\pi {\Delta }_i^{p_0}}(1+{\Delta }_{i}K)\left[\underset{k}{sup}{e}_{n-i-1,k}+2\underset{k}{sup}{\dot{e}}_{n-i-1,k}\right]\hfill \end{array}$$

(140)

recursively for $1<j\le q+1$ using the Lipschitz property of the driver f and the boundedness of the Runge–Kutta coefficients. Equations (137) and (138) combined with Equations (140) and (139) lead to

$$\begin{array}{cc}\hfill \underset{k}{sup}{e}_{i,k}+\underset{k}{sup}{\dot{e}}_{i,k}& \le 2\underset{i,k}{sup}{E}_{ik}+\frac{C_0\Delta x}{\pi {\Delta }_i^{p_0}}(1+{\Delta }_{n-i}K)\left(\underset{k}{sup}{e}_{i-1,k}+\underset{k}{sup}{\dot{e}}_{i-1,k}\right)\hfill \\ \hfill & \le 2\underset{i,k}{sup}{E}_{ik}+\zeta (1+{\Delta }_{n-i}K)\left(\underset{k}{sup}{e}_{i-1,k}+\underset{k}{sup}{\dot{e}}_{i-1,k}\right)\hfill \end{array}$$

where

$$\underset{i}{sup}\left\{\frac{C_0\Delta x}{\pi {\Delta }_i^{p_0}}\right\}\le \zeta \le 1.$$

Gronwall’s Lemma then yields

$$\underset{k}{sup}{e}_{i,k}+\underset{k}{sup}{\dot{e}}_{i,k}\le 2e^{TK}\underset{i,k}{sup}{E}_{ik}$$

(141)

so that the scheme is stable. The result of Equation (136) follows by taking the supremum on the left hand side of Equation (141) over time steps and applying Theorem 7.   □

In this general case, the global discretization error maintains the structure of the local discretization error under a stability condition. Equation (135) indicates that the space discretization has to be relatively as fine as the time discretization to ensure stability. Hence, stability can always be reached for any time discretization by refining the space discretization. However, the structure of the characteristic functions ${\varphi }_i$ and ${\mathsf{\Phi }}_{ij}$ determines the relative refinement needed for the space discretization.

4. Numerical Results

We test the convergence properties of the Fourier interpolation method on Runge–Kutta schemes with a problem of commodity derivative pricing under a model proposed by Lucia and Schwartz (2002). We shall test the method’s convergence and behaviour on smooth and unbounded FBSDE coefficients.

The commodity spot price X is defined by

$$X_t=e^{S\left(t\right)+V_t}$$

(142)

where the deterministic function $S:{\mathbb{R}}^{+}\to \mathbb{R}$ represents the seasonality component of the commodity and V is the price diffusion following an Ornstein-Uhlenbeck process according to the Vašíček (1977) model

$$d{V}_t=-\kappa V_{t}dt+\sigma d{W}_t.$$

(143)

As indicated by Lucia and Schwartz (2002), the commodity spot price X satisfies the stochastic differential equation

$$d{X}_t=\kappa (\theta \left(t\right)-ln{X}_t)X_{t}dt+\sigma X_{t}d{W}_t$$

(144)

where

$$\theta \left(t\right)=\frac{1}{\kappa }\left(\frac{{\sigma }^2}{2}+\frac{dS}{dt}\left(t\right)\right)+S\left(t\right).$$

(145)

We consider the commodity price as our forward process through equation (144).

When the risk-free rate r and the market price of risk $\lambda$ are both constant, the forward (or future) price $F_{t,T}:=Y_t=u(t,X_t)$ with maturity $T>0$ at time $t<T$ is given by

$$\begin{array}{cc}\hfill Y_t& ={\mathbf{E}}_t^{\mathbf{Q}}\left[X_T\right]\hfill \\ \hfill & =exp\left(S\left(T\right)+(ln{X}_t-S\left(t\right))e^{-\kappa (T-t)}-\frac{\sigma \lambda }{\kappa }h(T-t,\kappa )+\frac{{\sigma }^2}{4\kappa }h(T-t,2\kappa )\right)\hfill \end{array}$$

(146)

with

$$h(\tau ,\kappa )=1-e^{-\kappa \tau }$$

(147)

where the expectation is taken under the equivalent risk measure $\mathbf{Q}$. It can be shown that the forward price solves a BSDE with linear driver

$$f(t,x,y,z)=-\lambda z$$

(148)

and terminal condition

$$g\left(x\right)=x.$$

(149)

Options on forward contracts can also be represented in form of BSDEs in this spot price model but we limit our analysis to forward price estimation. From Equation (146) the control process (or equivalently the forward price delta) is given by

$$\begin{array}{cc}\hfill Z_t& =\sigma X_t\nabla u(t,X_t)\hfill \\ \hfill & =\sigma e^{-\kappa (T-t)}u(t,X_t).\hfill \end{array}$$

(150)

That is, we wish to use the results of Section 2 and Section 3 to numerically solve the FBSDE

$$\begin{array}{cc}\hfill d{X}_t& =\kappa (\theta \left(t\right)-ln{X}_t)X_{t}dt+\sigma X_{t}d{W}_t\hfill \\ \hfill -d{Y}_t& =-\lambda Z_{t}dt-Z_{t}d{W}_t\hfill \\ \hfill X_0& =x_0,Y_T=X_T\hfill \end{array}$$

and compare our numerical solutions under different discretization schemes to the explicit solution given by Equations (146) and (150).

The adjustment speed of the diffusion process is $\kappa =1.5$ and the volatility of the diffusion is set to be $\sigma =0.065$. The seasonality component is given by

$$S\left(t\right)=ln\overline{P}+0.05sin\left(2\pi t\right)$$

(151)

and the initial spot price by

$$X_0=\overline{P}{e}^{V_0}=0.95\overline{P}$$

(152)

where we normalize the real value3 of the commodity $\overline{P}=1$. Furthermore, the maturity of the forward contract is $T=0.25$ and we suppose a market price of risk of $\lambda =0.25$.

The FBSDE is solved on an alternative grid centred at $X_0$ with a uniform time mesh. For a given number of time steps n and the initial number $N_0=1$ of intervals, the length of an increment interval is set as

$$l=\frac{1.8}{N_0+n}$$

(153)

so that the truncated interval at time $t_n$ has length $1.8$. This restriction keeps the space nodes in the upper half plane knowing that the commodity price is a positive process. Moreover, the number of space steps on an increment interval is $N=2$.

We numerically solve the BSDE with the explicit $1-$stage Runge–Kutta scheme of half-order and an explicit $2-$stage Runge–Kutta scheme of first-order. Under the explicit $1-$stage scheme, the commodity price is discretized with an Euler scheme whereas a Milstein scheme is used for the forward process X under the explicit $2-$stage Runge–Kutta scheme. In addition, we use an explicit $2-$stage Runge–Kutta scheme with tableau

0	0	0	0	0	0
$\frac{2}{3}$	$\frac{2}{3}$	0	0	$\frac{2}{3}$	0
1	$\frac{1}{4}$	$\frac{3}{4}$	0	1	0

Under both FBSDE discretizations, we compute two different types of error. The first error $E_{True}$ evaluates the maximal absolute error of the numerical solution with respect to the true solution

$$\begin{array}{cc}\hfill E_{True}& =\underset{0\le i<n}{max}\underset{0\le k\le N{N}_i}{max}\left|u(t_i,x_{ik})-u_{ik}\right|+\underset{0\le i<n}{max}\underset{0\le k\le N{N}_i}{max}\left|\dot{u}(t_i,x_{ik})-{\dot{u}}_{ik}\right|\hfill \end{array}$$

(154)

where

$$\dot{u}(t,x)=\sigma x\nabla u(t,x)=\sigma e^{-\kappa (T-t)}u(t,x).$$

(155)

The second error $E_{Sim}$ is a simulation error. Given the numerical solution ${\left\{X_{t_i,j}^{\pi }\right\}}_{j=1}^m$, $i=0,1,\dots ,n-1$ with $m>0$ simulated paths for the forward process, we compute the numerical solution ${\left\{(y_{t_i,j},z_{t_i,j})\right\}}_{j=1}^m$ of the backward processes by linearly interpolating the simulated paths through the BSDE numerical solutions ${\left\{u_{ik}\right\}}_{k=0}^{N_{i}N}$ and ${\left\{{\dot{u}}_{ik}\right\}}_{k=0}^{N_{i}N}$ at each time step $t_i$. The error $E_{Sim}$ can be written as

$$\begin{array}{cc}\hfill E_{Sim}& =\frac{1}{m}\sum _{j=1}^m\left[\underset{0\le i<n}{max}\left|u(t_i,X_{t_i,j}^{\pi })-y_{t_i,j}\right|+{\left(\sum _{i=0}^{n-1}{\Delta }_i{(\dot{u}(t_i,X_{t_i,j}^{\pi })-z_{t_i,j})}^2\right)}^{\frac{1}{2}}\right].\hfill \end{array}$$

We systematically use $m=1000$ paths. Even if the errors $E_{True}$ and $E_{Sim}$ may be of the same order, they are interpreted differently. The error $E_{True}$ gives the behaviour of the maximal approximation error on the grid whereas $E_{Sim}$ gives the behaviour of the error on the relevant part of grid when solving the FBSDE numerically. Figure 2 displays the errors under the explicit $1-$stage Runge–Kutta scheme with $n\in \{5,10,20,50,100\}$ and Figure 3 shows the errors under the explicit $2-$stage scheme.

The error graphs of Figure 2 and Figure 3 look almost identical and confirm that the $2-$stage scheme is of first order and the $1-$stage scheme of (at least) half-order. The extra-efficiency of the $1-$stage scheme may be attributed in this particular case to the simplicity of the driver f and the terminal condition g.

In Figure 4, we present the absolute errors along the simulated paths for the BSDE solution. One notices that the maximal errors occur at the initial time $t_0=0$ for the forward price ($Y_t$) and at maturity $T=0.25$ for the control process ($Z_t$). Nonetheless, the simulation errors are of the same order (${10}^{-4}$) for both processes. This information is confirmed by the contour plot of Figure 5 not only along the simulated paths but on the entire grid.

Moreover, the contour plot gives indication on the source of errors. Indeed, Figure 5 shows that the maximal errors mainly occur for the upper space node values on the alternative grid and they decrease for lower space node values. This is due to the unbounded nature of the spot price process coefficients. Since the volatility of the spot price is a positive and increasing function of the spot price4, higher spot price values lead to higher local volatility. Hence, the fixed length of increment interval l may not be sufficiently large to ensure accuracy for higher space node values. In general, the phenomenon is amplified with the magnitude of the forward process coefficients as illustrated in the contour plot of Figure 6 where we choose a higher value for the volatility $\sigma$ and keep the other parameters unchanged. Similar results can be obtained by selecting a higher value for the speed of adjustment $\kappa$ as shown in Figure 7.

We end this section with an efficiency study of our schemes. Using the parameters initially given, the BSDE is solved on a uniform time grid with $n\in \{10,20,40,50,60,80,100\}$ time steps and $N\in \{2,2^2,2^3,2^4\}$ space steps and value the computation time. Figure 8 displays the results. First note that since the Fourier interpolation method performs matrix multiplications, it is much slower than the convolution method of Hyndman and Oyono Ngou (2017).

As shown in Figure 8, the computation time of Fourier interpolation method increases with the number of time steps leading to a trade-off between computation speed and accuracy. The exponential nature of the curves suggests that preference has to be given to the coarsest time discretization providing a satisfactory level of accuracy. Similarly, the computation time also increase drastically with the number N of space steps. Coarse space grids insuring accuracy are hence also preferable. Since a total number of $2q$ conditional expectations are computed under a q-stage Runge–Kutta scheme, we can expect the 1-stage scheme to run twice as fast as the 2-stage scheme. This is confirmed on Figure 8, especially when looking at the computation times for $n=100$.

5. Conclusions

In order to solve the problem of extrapolation errors in the initial implementation of the convolution method, we proposed an alternative space discretization. The new tree-like space grid naturally allows the usage of the FFT algorithm when computing the conditional expectation included in the underlying explicit Euler scheme. The error analysis shows that both the alternative grid and the (alternative) transform suit the periodic nature of the FFT algorithm and help in producing a stable, consistent and globally convergent numerical procedure for the FBSDE approximate solutions. The second part of the paper deals with the implementation of the Fourier interpolation method with higher order time discretizations of FBSDEs. When the forward process increments admit conditional characteristic functions satisfying certain regularity conditions, it was shown that the method is also consistent, conditionally stable and globally convergent under Runge–Kutta schemes for FBSDEs. A challenging area of research is the implementation of the methods of this paper in multidimensional and jump cases.

References [custom]