Abstract

A Multi-step Scheme based on Cubic Spline for solving Backward Stochastic Differential Equations

Long Teng^1,^*^* * Corresponding author (teng@math.uni-wuppertal.de), Aleksandr Lapitckii, Michael GÜNTHER¹

$\mbox{}^{{1}}$ Lehrstuhl für Angewandte Mathematik und Numerische Analysis,

Fakultät für Mathematik und Naturwissenschaften,

Bergische Universität Wuppertal, Gaußstr. 20, 42119 Wuppertal, Germany

Abstract

In this work we study a multi-step scheme on time-space grids proposed by W. Zhao et al. [Zhao et al., 2010] for solving backward stochastic differential equations, where Lagrange interpolating polynomials are used to approximate the time-integrands with given values of these integrands at chosen multiple time levels. For a better stability and the admission of more time levels we investigate the application of spline instead of Lagrange interpolating polynomials to approximate the time-integrands. The resulting scheme is a semi-discretization in the time direction involving conditional expectations, which can be numerically solved by using the Gaussian quadrature rules and polynomial interpolations on the spatial grids. Several numerical examples including applications in finance are presented to demonstrate the high accuracy and stability of our new multi-step scheme.

Keywords backward stochastic differential equations, multi-step scheme, cubic splines, time-space grid, Gauss-Hermite quadrature rule

1 Introduction

Recently, the forward-backward stochastic differential equation (FBSDE) becomes an important tool for formulating many problems in, e.g., mathematical finance and stochastic control. The BSDE exhibits usually no analytical solution, see e.g., [Karoui et al., 1997a]. Their numerical solutions have thus been extensively studied by many researchers. The general form of (decoupled) FBSDEs reads

\left\{\begin{array}[]{l}\,\,\,dX_{t}=a(t,X_{t})\,dt+b(t,X_{t})\,dW_{t},\quad X_{0}=x_{0},\\ -dY_{t}=f(t,X_{t},Y_{t},Z_{t})\,dt-Z_{t}\,dW_{t},\\ \quad Y_{T}=\xi=g(X_{T}),\end{array}\right.

(1)

where $X_{t},a\in\mathbb{R}^{n},$ $b$ is a $n\times d$ matrix, $W_{t}=(W^{1}_{t},\cdots,W^{d}_{t})^{T}$ is a $d$ -dimensional Brownian motion (all Brownian motions are independent with each other), $f(t,X_{t},Y_{t},Z_{t}):[0,T]\times\mathbb{R}^{n}\times\mathbb{R}^{m}\times\mathbb{R}^{m\times d}\to\mathbb{R}^{m}$ is the driver function and $\xi$ is the square-integrable terminal condition. We see that the terminal condition $Y_{T}$ depends on the final value of a forward stochastic differential equation (SDE).

For $a=0~\mbox{and}~b=1,$ namely $X_{t}=W_{t},$ one obtains a backward stochastic differential equation (BSDE) of the form

\left\{\begin{array}[]{l}-dY_{t}=f(t,Y_{t},Z_{t})\,dt-Z_{t}\,dW_{t},\\ \quad Y_{T}=\xi=g(W_{T}),\end{array}\right.

(2)

where $Y_{t}\in\mathbb{R}^{m}$ and $f:[0,T]\times\mathbb{R}^{m}\times\mathbb{R}^{m\times d}\to\mathbb{R}^{m}.$ In the sequel of this paper, we investigate the numerical scheme for solving (2). Note that the developed schemes can be applied also for solving (1), where the general Markovian diffusion $X_{t}$ can be approximated, e.g., by using the Euler-Scheme.

The existence and uniqueness of solution of (2) assuming the Lipschitz conditions on $f,a(t,X_{t}),b(t,X_{t})~\mbox{and}~g$ are proven by Pardoux and Peng [Pardoux and Peng, 1990, Pardoux and Peng, 1992]. The uniqueness of solution is extended under more general assumptions for $f$ in [Lepeltier and Martin, 1997], but only in the one-dimensional case.

In recent years, many numerical methods have been proposed for the FBSDEs and BSDEs. Peng [Peng, 1991] obtained a direct relation between FBSDEs and partial differential equations (PDEs), see also [Karoui et al., 1997b]. Based on this relation, several numerical schemes are proposed, e.g., [Douglas et al., 1996, Ma et al., 1994, Milsetin and Tretyakov, 2006]. As probabilistic methods, (least-squares) Monte-Carlo approaches are investigated in [Bender and Steiner, 2012, Bouchard and Touzi, 2004, Gobet et al., 2005, Lemor et al., 2006, Zhao et al., 2006], and tree-based approaches in [Crisan and Manolarakis, 2010, Teng, 2018]. For numerical approximation and analysis we refer to [Bally, 1997, Bender and Zhang, 2008, Ma et al., 2009, Ma and Zhang, 2005, Zhang, 2004, Zhao et al., 2010]. And many others, e.g., some numerical methods for BSDEs applying binomial tree are investigated in [Ma et al., 2002]. The approach based on the Fourier method for BSDEs is developed in [Ruijter and Oosterlee, 2015].

In [Zhao et al., 2010], a multi-step scheme is achieved by using Lagrange interpolating polynomials. However, the number of multiple time levels is restricted, the stability condition cannot be satisfied for a high number of time steps. This is actually to be expected due to Runge’s phenomenon. For this reason, we study in this work a stable multi-step scheme by using the cubic spline polynomials, for numerically solving BSDEs on the time-space grids. More precisely, we use the cubic spline polynomials to approximate the integrands, which are conditional mathematical expectations derived from the original BSDEs. For this, we need to know values of integrands at multiple time levels, which can be numerically evaluated, e.g., using the Gauss-Hermite quadrature and polynomial interpolations on the spatial grids. We will study the convergence and the error estimates for the proposed multi-step scheme.

In the next section, we start with notation and definitions and derive in Section 3 the reference equations for our multi-step scheme for the BSDEs. In Section 4, we introduce the multi-step scheme for their discretizations. Section 5 is devoted to error estimates. In Section 6, several numerical experiments on different types of (F)BSDEs including financial applications are provided to show the high accuracy and stability. Finally, Section 7 concludes this work.

2 Preliminaries

Throughout the paper, we assume that $(\Omega,\mathcal{F},P;\{\mathcal{F}_{t}\}_{0\leq t\leq T})$ is a complete, filtered probability space. In this space, a standard $d$ -dimensional Brownian motion $W_{t}$ with a finite terminal time $T$ is defined, which generates the filtration $\{\mathcal{F}_{t}\}_{0\leq t\leq T},$ i.e., $\mathcal{F}_{t}=\sigma\{X_{s},0\leq s\leq t\}$ for FBSDEs or $\mathcal{F}_{t}=\sigma\{W_{s},0\leq s\leq t\}$ for BSDEs. And the usual hypotheses should be satisfied. We denote the set of all $\mathcal{F}_{t}$ -adapted and square integrable processes in $\mathbb{R}^{d}$ with $L^{2}=L^{2}(0,T;\mathbb{R}^{d}).$ A pair of process $(Y_{t},Z_{t}):[0,T]\times\Omega\to\mathbb{R}^{m}\times\mathbb{R}^{m\times d}$ is the solution of the BSDE (2) if it is $\mathcal{F}_{t}$ -adapted and square integrable and satisfies (2) as

Y_{t}=\xi+\int_{t}^{T}f(s,Y_{s},Z_{s})\,ds-\int_{t}^{T}Z_{s}\,dW_{s},\quad t\in[0,T],

(3)

where $f(t,Y_{s},Z_{s}):[0,T]\times\mathbb{R}^{m}\times\mathbb{R}^{m\times d}\to\mathbb{R}^{m}$ is $\mathcal{F}_{t}$ adapted, $g:\mathbb{R}^{d}\to\mathbb{R}^{m}.$ As mentioned above, these solutions exist uniquely under Lipschitz conditions.

Suppose that the terminal value $Y_{T}$ is of the form $g(W^{t,x}_{T}),$ where $W^{t,x}_{T}$ denotes the value of $W_{T}$ starting from $x$ at time $t.$ Then the solution $(Y^{t,x}_{t},Z^{t,x}_{t})$ of BSDEs (2) can be represented [Karoui et al., 1997b, Ma and Zhang, 2005, Pardoux and Peng, 1992, Peng, 1991] as

Y^{t,x}_{t}=u(t,x),\quad Z^{t,x}_{t}=\nabla u(t,x)\quad\forall t\in[0,T),

(4)

which is the solution of the semilinear parabolic PDE of the form

\frac{\partial u}{\partial t}+\frac{1}{2}\sum_{i}^{d}\partial^{2}_{i,i}u+f(t,u,\nabla u)=0

(5)

with the terminal condition $u(T,x)=g(x).$ In turn, suppose $(Y,Z)$ is the solution of BSDEs, $u(t,x)=Y^{t,x}_{t}$ is a viscosity solution to the PDE.

3 Reference equations for the multi-step scheme

In this section we drive the reference equations for the multi-step scheme by using the cubic spline polynomials.

3.1 The one-dimensional reference equations

We start with the one-dimensional processes, namely $m=n=d=1.$ We introduce the uniform time partition for the time interval $[0,T]$

\Delta_{t}=\{t_{i}|t_{i}\in[0,T],i=0,1,\cdots,N_{T},t_{i}<t_{i+1},t_{0}=0,t_{N_{T}}=T\}.

(6)

Let $\Delta t:=h=\frac{T}{N_{T}}$ be the time step, and thus $t_{i}=t_{0}+ih,~\mbox{for}~i=0,1,\cdots,N_{T}.$ Then one needs to discretize the backward process (3), namely

Y_{t}=\xi+\int_{t}^{T}f(s,\mathbb{V}_{s})\,ds-\int_{t}^{T}Z_{s}\,dW_{s},

(7)

where $\xi=g(W_{T}),\mathbb{V}_{s}=(Y_{s},Z_{s}).$ Let $(Y_{t},Z_{t})$ be the adapted solution of (7), we thus have

Y_{i}=Y_{i+k}+\int_{t_{i}}^{t_{i+k}}f(s,\mathbb{V}_{s})\,ds-\int_{t_{i}}^{t_{i+k}}Z_{s}\,dW_{s},\quad t\in[0,T),

(8)

where $1\leq k\leq K_{y}\leq N_{T}$ with two given positive integers $k~\mbox{and}~K_{y}.$ To obtain the adaptability of the solution $(Y_{t},Z_{t}),$ we use conditional expectations $E_{i}[\cdot](=E[\cdot|\mathcal{F}_{t_{i}}]).$ We start finding the reference equation for $Y.$ We take the conditional expectations $E_{i}[\cdot]$ on the both sides of (8) to obtain

Y_{i}=E_{i}[Y_{i+k}]+\int_{t_{i}}^{t_{i+k}}E_{i}[f(s,\mathbb{V}_{s})]\,ds.

(9)

We see that the integrand on the right-hand side of (9) is deterministic of time $s.$ When the values of $\mathbb{V}_{s},$ $(y_{t},z_{t})$ are available on the time levels $t_{i+1},t_{i+2},\cdots,t_{i+K_{y}},$ an approximation of the integrand in (9) can be found. In this work we choose the cubic spline interpolant $\tilde{S}_{K_{y},t_{i}}(s)$ based on the support values $(t_{i+j},E_{i}[f(t_{i+j},Y_{i+j},Z_{i+j})]),j=0,\cdots,K_{y},$ namely we have

\int_{t_{i}}^{t_{i+k}}E_{i}[f(s,\mathbb{V}_{s})]\,ds=\int_{t_{i}}^{t_{i+k}}\tilde{S}_{K_{y},t_{i}}(s)\,ds+R_{y}^{i}

(10)

with the residual

R_{y}^{i}=\int_{t_{i}}^{t_{i+k}}\left(E_{i}[f(s,\mathbb{V}_{s})]-\tilde{S}_{K_{y},t_{i}}(s)\right)\,ds.

(11)

Then we can calculate

\int_{t_{i}}^{t_{i+k}}\tilde{S}_{K_{y},t_{i}}(s)\,ds=\int_{t_{i}}^{t_{i+k}}\sum_{j=0}^{K_{y}-1}\tilde{s}^{y}_{t_{i},j}(s)\,ds=\sum_{j=0}^{K_{y}-1}\int_{t_{i}}^{t_{i+k}}\tilde{s}^{y}_{t_{i},j}(s)\,ds

(12)

with

\tilde{s}^{y}_{t_{i},j}(s)=a^{y}_{j}+b^{y}_{j}(s-t_{i+j})+c^{y}_{j}(s-t_{i+j})^{2}+d^{y}_{j}(s-t_{i+j})^{3},

(13)

where $s\in[t_{i+j},t_{i+j+1}],\,j=0,\cdots,K_{y}-1.$ We straightforwardly calculate

\begin{split}\int_{t_{i}}^{t_{i+k}}\tilde{s}^{y}_{t_{i},j}\,ds&=\int_{t_{i+j}}^{t_{i+j+1}}\tilde{s}^{y}_{t_{i},j}(s)\,ds\\ &=a^{y}_{j}h+\frac{b^{y}_{j}h^{2}}{2}+\frac{c^{y}_{j}h^{3}}{3}+\frac{d^{y}_{j}h^{4}}{4}.\end{split}

(14)

Note that $j$ satisfying $k-1<j\leq K_{y}-1$ results an integral with zero value when $k<K_{y}.$ And the coefficients $a^{y}_{j},b^{y}_{j},c^{y}_{j}~\mbox{and}~d^{y}_{j}$ are obtained with the support points $(t_{i+j},E_{i}[f(t_{i+j},Y_{i+j},Z_{i+j})]),j=0,\cdots,K_{y}$ as

\begin{cases}\tilde{S}_{K_{y},t_{i}}(t_{i+j})=E_{i}[f(t_{i+j},Y_{{i+j}},Z_{{i+j}})]\quad&j=0,...,K_{y}\\ \tilde{s}^{y}_{t_{i},j}(t_{i+j})=\tilde{s}^{y}_{t_{i},j+1}(t_{i+j})\quad&j=0,1,...,K_{y}-2\\ \tilde{s}^{{}^{\prime}y}_{t_{i},j}(t_{i+j})=\tilde{s}^{{}^{\prime}y}_{t_{i},j+1}(t_{i+j})\quad&j=0,1,...,K_{y}-2\\ \tilde{s}^{{}^{\prime\prime}y}_{t_{i},j}(t_{i+j})=\tilde{s}^{{}^{\prime\prime}y}_{t_{i},j+1}(t_{i+j})\quad&j=0,1,...,K_{y}-2.\end{cases}

(15)

Obviously, we need two boundary conditions to solve the system above. Since the values of derivatives of $E_{i}[f(t_{i+j},Y_{i+j},Z_{i+j})]$ are unknown, we could thus choose e.g., the natural boundary conditions or Not-a-Knot conditions depending on the value of $K_{y}.$ Combining (9), (10), (12) and (14) we obtain the reference equation for $Y_{i}$ (based on those support points) as:

Y_{i}=E_{i}[Y_{i+k}]+\sum_{j=0}^{Ky-1}\left[a^{y}_{j}h+\frac{b^{y}_{j}h^{2}}{2}+\frac{c^{y}_{j}h^{3}}{3}+\frac{d^{y}_{j}h^{4}}{4}\right]+R_{y}^{i},

(16)

where the coefficients $a_{j}^{y},b_{j}^{y},c_{j}^{y}~\mbox{and}~d_{j}^{y}$ will be obtained by solving (15) together with appropriate boundary conditions and depend on $Y_{i}.$ Therefore, (16) is an implicit scheme.

We now start with the reference equation for $Z.$ By multiplying both sides of the equation (8) by $\Delta W_{i+1}:=W_{t_{i+1}}-W_{t_{i}}$ and taking the conditional expectations $E_{i}[\cdot]$ on both sides of the derived equation we obtain

-E_{i}[Y_{i+l}\Delta W_{i+l}]=\int_{t_{i}}^{t_{i+l}}E_{i}[f(s,\mathbb{V}_{s})\Delta W_{s}]\,ds-\int_{t_{i}}^{t_{i+l}}E_{i}[Z_{s}]\,ds,

(17)

where the Itô isometry and Fubini’s theorem are used, $\Delta W_{s}=W_{s}-W_{t_{i}}$ and the given integers $l~\mbox{and}~K_{z}$ satisfy $1\leq l\leq K_{z}.$ Similarly, we derive the reference equation of $Z$ also based on the support points $(t_{i+j},E_{i}[f(t_{i+j},y_{i+j},z_{i+j})\Delta w_{i+j}])$ and $((t_{i+j},E_{i}[z_{i+j}])$ , $j=0,\cdots,K_{z}.$ Then, we again use the cubic spline polynomials to approximate the time deterministic integers and obtain

\begin{split}\int_{t_{i}}^{t_{i+l}}E_{i}[f(t_{s},Y_{s},Z_{s})\Delta w_{s}]\,ds&=\int_{t_{i}}^{t_{i+l}}\tilde{S}_{K_{z_{1}},t_{i}}(s)\,ds+R_{z_{1}}^{i}\\ &=\sum_{j=0}^{K_{z}-1}\int_{t_{i}}^{t_{i+l}}\tilde{s}^{z_{1}}_{t_{i},j}(s)\,ds+R_{z_{1}}^{i}\end{split}

(18)

with

	$\displaystyle R_{z_{1}}^{i}$	$\displaystyle=\int_{t_{i}}^{t_{i+l}}\left(E_{i}[f(t_{s},Y_{s},Z_{s})\Delta w_{s}]-\tilde{S}_{K_{z_{1}},t_{i}}(s)\right)\,ds,$		(19)
	$\displaystyle\tilde{s}^{z_{1}}_{t_{i},j}(s)$	$\displaystyle=a^{z_{1}}_{j}+b^{z_{1}}_{j}(s-t_{i+j})+c^{z_{1}}_{j}(s-t_{i+j})^{2}+d^{z_{1}}_{j}(s-t_{i+j})^{3}$		(20)

$~\mbox{for}~s\in[t_{i+j},t_{i+j+1}],\,j=0,\cdots,K_{z}-1,$ and

\begin{split}\int_{t_{i}}^{t_{i+l}}E_{i}[Z_{s}]\,ds&=\int_{t_{i}}^{t_{i+l}}\tilde{S}_{K_{z_{2}},t_{i}}(s)\,ds+R_{z_{2}}^{i}\\ &=\sum_{j=0}^{K_{z}-1}\int_{t_{i}}^{t_{i+l}}\tilde{s}^{z_{2}}_{t_{i},j}(s)\,ds+R_{z_{2}}^{i}\end{split}

(21)

with

	$\displaystyle R_{z_{2}}^{i}$	$\displaystyle=\int_{t_{i}}^{t_{i+l}}\left(E_{i}[Z_{s}]-\tilde{S}_{K_{z_{2}},t_{i}}(s)\right)\,ds,$		(22)
	$\displaystyle\tilde{s}^{z_{2}}_{t_{i},j}(s)$	$\displaystyle=a^{z_{2}}_{j}+b^{z_{2}}_{j}(s-t_{i+j})+c^{z_{2}}_{j}(s-t_{i+j})^{2}+d^{z_{2}}_{j}(s-t_{i+j})^{3}$		(23)

$~\mbox{for}~s\in[t_{i+j},t_{i+j+1}],\,j=0,\cdots,K_{z}-1$ and we let

R_{z}^{i}:=R_{z_{1}}^{i}+R_{z_{2}}^{i}.

(24)

Furthermore, using the relation (4) and integration by parts it can be verified that

E_{i}[Y_{i+l}\Delta W_{i+l}]=lhE_{i}[Z_{i+1}].

(25)

Integrating (20), (23) and combining (17), (18), (21) and (25) we obtain the reference equation for $Z_{i}$ as:

\begin{split}0=lhE_{i}[Z_{i+l}]&+\sum_{j=0}^{Kz-1}\left[a^{z_{1}}_{j}h+\frac{b^{z_{1}}_{j}h^{2}}{2}+\frac{c^{z_{1}}_{j}h^{3}}{3}+\frac{d^{z_{1}}_{j}h^{4}}{4}\right]\\ &-\sum_{j=0}^{Kz-1}\left[a^{z_{2}}_{j}h+\frac{b^{z_{2}}_{j}h^{2}}{2}+\frac{c^{z_{2}}_{j}h^{3}}{3}+\frac{d^{z_{2}}_{j}h^{4}}{4}\right]+R_{z}^{i},\end{split}

(26)

where the coefficients $a^{z_{1}}_{j},b^{z_{1}}_{j},c^{z_{1}}_{j},d^{z_{1}}_{j}$ are solutions of

\begin{cases}\tilde{S}_{K_{z},t_{i}}(t_{i+j})=E_{i}[f(t_{i+j},Y_{i+j},Z_{i+j})\Delta W_{i+j}]\quad&j=0,...,K_{z}\\ \tilde{s}^{{z_{1}}}_{t_{i},j}(t_{i+j})=\tilde{s}^{{z_{1}}}_{t_{i},j+1}(t_{i+j})\quad&j=0,...,K_{z}-2\\ \tilde{s}^{{}^{\prime}{z_{1}}}_{t_{i},j}(t_{i+j})=\tilde{s}^{{}^{\prime}{z_{1}}}_{t_{i},j+1}(t_{i+j})\quad&j=0,...,K_{z}-2\\ \tilde{s}^{{}^{\prime\prime}{z_{1}}}_{t_{i},j}(t_{i+j})=\tilde{s}^{{}^{\prime\prime}{z_{1}}}_{t_{i},j+1}(t_{i+j})\quad&j=0,...,K_{z}-2\end{cases}

(27)

with the appropriate boundary conditions, and the coefficients $a^{z_{2}}_{j},b^{z_{2}}_{j},c^{z_{2}}_{j},d^{z_{2}}_{j}$ are solutions of

\begin{cases}\tilde{S}_{K_{z},t_{i}}(t_{i+j})=E_{i}[Z_{{i+j}}]\quad&j=0,...,K_{z}\\ \tilde{s}^{{z_{2}}}_{t_{i},j}(t_{i+j})=\tilde{s}^{z_{2}}_{t_{i},j+1}(t_{i+j})\quad&j=0,...,K_{z}-2\\ \tilde{s}^{{}^{\prime}{z_{2}}}_{t_{i},j}(t_{i+j})=\tilde{s}^{{}^{\prime}{z_{2}}}_{t_{i},j+1}(t_{i+j})\quad&j=0,...,K_{z}-2\\ \tilde{s}^{{}^{\prime\prime}{z_{2}}}_{t_{i},j}(t_{i+j})=\tilde{s}^{{}^{\prime\prime}{z_{2}}}_{t_{i},j+1}(t_{i+j})\quad&j=0,...,K_{z}-2\end{cases}

(28)

with the appropriate boundary conditions, respectively.

3.2 The high-dimensional reference equations

In this section, we give the reference equations for the high-dimensional case. With the aid of (16) we can straightforwardly write the reference equation for $y_{i}$ in component-wise as

Y^{\tilde{m}}_{i}=E_{i}[Y^{\tilde{m}}_{i+k}]+\sum_{j=0}^{Ky-1}\left[a^{y_{j},\tilde{m}}h+\frac{b^{y_{j},\tilde{m}}h^{2}}{2}+\frac{c^{y_{j},\tilde{m}}h^{3}}{3}+\frac{d^{y_{j},\tilde{m}}h^{4}}{4}\right]+R_{y}^{i,\tilde{m}},

(29)

with

\begin{cases}\tilde{S}^{\tilde{m}}_{K_{y},t_{i}}(t_{i+j})=\mathbb{E}_{i}[f^{\tilde{m}}(t_{i+j},Y_{{i+j}},Z_{{i+j}})]\quad&j=0,...,K_{y}\\ \tilde{s}^{y,\tilde{m}}_{t_{i},j}(t_{i+j})=\tilde{s}^{y,\tilde{m}}_{t_{i},j+1}(t_{i+j})\quad&j=0,1,...,K_{y}-2\\ \tilde{s}^{{}^{\prime}y,\tilde{m}}_{t_{i},j}(t_{i+j})=\tilde{s}^{{}^{\prime}y,\tilde{m}}_{t_{i},j+1}(t_{i+j})\quad&j=0,1,...,K_{y}-2\\ \tilde{s}^{{}^{\prime\prime}y,\tilde{m}}_{t_{i},j}(t_{i+j})=\tilde{s}^{{}^{\prime\prime}y,\tilde{m}}_{t_{i},j+1}(t_{i+j})\quad&j=0,1,...,K_{y}-2,\end{cases}

(30)

where $f^{\tilde{m}}$ is the $\tilde{m}$ -th component of the vector $f$ for $\tilde{m}=1,2,\cdots,m.$ The coefficients $a_{j}^{y,\tilde{m}},b_{j}^{y,\tilde{m}},c_{j}^{y,\tilde{m}}~\mbox{and}~d_{j}^{y,\tilde{m}}$ will be obtained by solving the $\tilde{m}$ -th system (30) together with appropriate boundary conditions. The $\tilde{m}$ -th component residual reads

R_{y}^{i,\tilde{m}}=\int_{t_{i}}^{t_{i+k}}\left(E_{i}[f^{\tilde{m}}(s,Y_{s},Z_{s})]-\tilde{S}^{\tilde{m}}_{K_{y},t_{i}}(s)\right)\,ds.

(31)

Similarly, the reference equation for $Z_{i}$ can be formulated as follows:

\begin{split}0=lhE_{i}[Z^{\tilde{m},\tilde{d}}_{i+l}]&+\sum_{j=0}^{Kz-1}\left[a^{z_{1},\tilde{m},\tilde{d}}_{j}h+\frac{b^{z_{1},\tilde{m},\tilde{d}}_{j}h^{2}}{2}+\frac{c^{z_{1},\tilde{m},\tilde{d}}_{j}h^{3}}{3}+\frac{d^{z_{1},\tilde{m},\tilde{d}}_{j}h^{4}}{4}\right]\\ &-\sum_{j=0}^{Kz-1}\left[a^{z_{2},\tilde{m},\tilde{d}}_{j}h+\frac{b^{z_{2},\tilde{m},\tilde{d}}_{j}h^{2}}{2}+\frac{c^{z_{2},\tilde{m},\tilde{d}}_{j}h^{3}}{3}+\frac{d^{z_{2},\tilde{m},\tilde{d}}_{j}h^{4}}{4}\right]+R_{z}^{i,{\tilde{m},\tilde{d}}},\end{split}

(32)

where the coefficients $a^{z_{1},\tilde{m},\tilde{d}}_{j},b^{z_{1},\tilde{m},\tilde{d}}_{j},c^{z_{1},\tilde{m},\tilde{d}}_{j},d^{z_{1},\tilde{m},\tilde{d}}_{j}$ are solutions of

\begin{cases}\tilde{S}^{\tilde{m},\tilde{d}}_{K_{z},t_{i}}(t_{i+j})=E_{i}[f^{\tilde{m}}(t_{i+j},Y_{i+j},Z_{i+j})\Delta W^{\tilde{d}}_{i+j}]\quad&j=0,...,K_{z}\\ \tilde{s}^{{z_{1}},\tilde{m},\tilde{d}}_{t_{i},j}(t_{i+j})=\tilde{s}^{{z_{1}},\tilde{m},\tilde{d}}_{t_{i},j+1}(t_{i+j})\quad&j=0,...,K_{z}-2\\ \tilde{s}^{{}^{\prime}{z_{1}},\tilde{m},\tilde{d}}_{t_{i},j}(t_{i+j})=\tilde{s}^{{}^{\prime}{z_{1}},\tilde{m},\tilde{d}}_{t_{i},j+1}(t_{i+j})\quad&j=0,...,K_{z}-2\\ \tilde{s}^{{}^{\prime\prime}{z_{1}},\tilde{m},\tilde{d}}_{t_{i},j}(t_{i+j})=\tilde{s}^{{}^{\prime\prime}{z_{1}},\tilde{m},\tilde{d}}_{t_{i},j+1}(t_{i+j})\quad&j=0,...,K_{z}-2\end{cases}

(33)

with the appropriate boundary conditions, and the coefficients $a^{z_{2},\tilde{m},\tilde{d}}_{j},b^{z_{2},\tilde{m},\tilde{d}}_{j},c^{z_{2},\tilde{m},\tilde{d}}_{j},d^{z_{2},\tilde{m},\tilde{d}}_{j}$ are solutions of

\begin{cases}\tilde{S}^{\tilde{m},\tilde{d}}_{K_{z},t_{i}}(t_{i+j})=E_{i}[Z^{\tilde{m},\tilde{d}}_{{i+j}}]\quad&j=0,...,K_{z}\\ \tilde{s}^{{z_{2}},\tilde{m},\tilde{d}}_{t_{i},j}(t_{i+j})=\tilde{s}^{z_{2},\tilde{m},\tilde{d}}_{t_{i},j+1}(t_{i+j})\quad&j=0,...,K_{z}-2\\ \tilde{s}^{{}^{\prime}{z_{2}},\tilde{m},\tilde{d}}_{t_{i},j}(t_{i+j})=\tilde{s}^{{}^{\prime}{z_{2}},\tilde{m},\tilde{d}}_{t_{i},j+1}(t_{i+j})\quad&j=0,...,K_{z}-2\\ \tilde{s}^{{}^{\prime\prime}{z_{2}},\tilde{m},\tilde{d}}_{t_{i},j}(t_{i+j})=\tilde{s}^{{}^{\prime\prime}{z_{2}},\tilde{m},\tilde{d}}_{t_{i},j+1}(t_{i+j})\quad&j=0,...,K_{z}-2.\end{cases}

(34)

The corresponding residual reads

R_{z}^{i,{\tilde{m},\tilde{d}}}=R_{z_{1}}^{i,{\tilde{m},\tilde{d}}}+R_{z_{2}}^{i,{\tilde{m},\tilde{d}}}

(35)

with

R_{z_{1}}^{i,{\tilde{m},\tilde{d}}}=\int_{t_{i}}^{t_{i+l}}\left(E_{i}[f^{{\tilde{m}}}(t_{s},Y_{s},Z_{s})\Delta W^{{\tilde{d}}}_{s}]-\tilde{S}^{{\tilde{m},\tilde{d}}}_{K_{z_{1}},t_{i}}(s)\right)\,ds,

(36)

R_{z_{2}}^{i,{\tilde{m},\tilde{d}}}=\int_{t_{i}}^{t_{i+l}}\left(E_{i}[Z^{\tilde{m},\tilde{d}}_{s}]-\tilde{S}^{\tilde{m},\tilde{d}}_{K_{z_{2}},t_{i}}(s)\right)\,ds,

(37)

where $\tilde{m}=1,2,\cdots,m$ and $\tilde{d}=1,2,\cdots,d.$ Note that, by removing superscripts $\tilde{m}$ and $\tilde{d},$ we can write (29) and (32) in matrix form.

3.3 The cubic spline coefficients

As mentioned before, due to the lack of derivative values of the integrands, we should choose some cubic spline which does not need those derivative values. Furthermore, it will be shown in the next section that (29) is stable for any positive $k$ and $K_{y},$ we thus fix $k=K_{y}.$ However, (32) is only stable for any positive $K_{z}$ and $l=1.$ Therefore, in the sequel of this paper we fix $k=K_{y}$ and $l=1.$

For the reference equation (15), we calculate cubic spline coefficients for different values of $K_{y}$ as follows. For notational simplicity, we let $g_{i+j}=E_{i}[f(t_{i+j},Y_{i+j},Z_{i+j})]$ for $j=0,\cdots,K_{y}.$

•

$K_{y}=1:$ there are only two points available. One can just construct a straight line and obtain $a^{y}_{0}=g_{i},b^{y}_{0}=\frac{g_{i+1}-g_{i}}{h},c^{y}_{0}=0,d^{y}_{0}=0.$ Now, we can rewrite (16) as

	$\displaystyle Y_{i}$	$\displaystyle=E_{i}[Y_{i+K_{y}}]+\frac{h}{2}g_{i}+\frac{h}{2}g_{i+1}+R_{y}^{i}$		(38)
		$\displaystyle:=E_{i}[Y_{i+K_{y}}]+hK_{y}\sum_{j=0}^{Ky}\gamma_{K_{y},j}^{K_{y}}E_{i}[f(t_{i+j},Y_{i+j},Z_{i+j})]+R_{y}^{i},$		(39)

where $\gamma_{K_{y},0}^{K_{y}}=\gamma_{K_{y},1}^{K_{y}}=\frac{1}{2}.$

•

$K_{y}=2:$ we can already construct e.g., a natural cubic spline based on three points. The corresponding coefficients can be calculated as follows.

For $\tilde{s}^{y}_{t_{i},0}(s),s\in[t_{i},t_{i+1}]:$

	$\displaystyle a_{0}$	$\displaystyle=g_{i},b_{0}=-(5g_{i}-6g_{i+1}+g_{i+2})/4h$
	$\displaystyle c_{0}$	$\displaystyle=0,d_{0}=(g_{i}-2g_{i+1}+g_{i+2})/4h^{3}$

For $\tilde{s}^{y}_{t_{i},1}(s),s\in[t_{i+1},t_{i+2}]:$

	$\displaystyle a_{1}$	$\displaystyle=g_{i+1},b_{1}=-(g_{i}-g_{i+2})/2h$
	$\displaystyle c_{1}$	$\displaystyle=(3g_{i}-6g_{i+1}+3g_{i+2})/4h^{2},d_{1}=-(g_{i}-2g_{i+1}+g_{i+2})/4h^{3}$

Thus, (16) can be rewritten as

	$\displaystyle Y_{i}$	$\displaystyle=E_{i}[Y_{i+K_{y}}]+\frac{3h}{8}g_{i}+\frac{10h}{8}g_{i+1}+\frac{3h}{8}g_{i+2}+R_{y}^{i}$		(40)
		$\displaystyle:=E_{i}[Y_{i+K_{y}}]+hK_{y}\sum_{j=0}^{Ky}\gamma_{K_{y},j}^{K_{y}}E_{i}[f(t_{i+j},Y_{i+j},Z_{i+j})]+R_{y}^{i},$		(41)

where $\gamma_{K_{y},0}^{K_{y}}=\gamma_{K_{y},2}^{K_{y}}=\frac{3}{16},\gamma_{K_{y},1}^{K_{y}}=\frac{5}{8}.$

Moreover, for the cubic spline we set the second derivatives of cubic interpolants at boundaries to be zero. Instead of this, one can also choose a second order polynomial for the whole interval, namely $(t_{i},t_{i+2}).$ In this way we obtain the polynomial $p_{i}(s)$ as

g_{i}(s-t_{i})-\left(\frac{3}{2}g_{i}-2g_{i+1}+\frac{1}{2}g_{i+2}\right)(s-t_{i})/h+\left(\frac{1}{2}g_{i}-g_{i+1}+\frac{1}{2}g_{i+2}\right)(s-t_{i})^{2}/h^{2}

(42)

and its integration as

\int_{t_{i}}^{t_{i+2}}p_{i}(s)ds=h\frac{g_{i}+4g_{i+1}+g_{i+2}}{3}.

(43)

By using the second order polynomial we rewrite (16) as

	$\displaystyle Y_{i}$	$\displaystyle=E_{i}[Y_{i+K_{y}}]+\frac{h}{3}g_{i}+\frac{4h}{3}g_{i+1}+\frac{h}{3}g_{i+2}+R_{y}^{i}$
		$\displaystyle:=E_{i}[Y_{i+K_{y}}]+hK_{y}\sum_{j=0}^{Ky}\gamma_{K_{y},j}^{K_{y}}E_{i}[f(t_{i+j},Y_{i+j},Z_{i+j})]+R_{y}^{i},$		(44)

where $\gamma_{K_{y},0}^{K_{y}}=\gamma_{K_{y},2}^{K_{y}}=\frac{1}{6},\gamma_{K_{y},1}^{K_{y}}=\frac{2}{3}.$

•

$K_{y}=3:$ for $K_{y}\geq 3$ we will use the Not-a-knot cubic spline and calculate the corresponding coefficients as follows.

For $\tilde{s}^{y}_{t_{i},0}(s),s\in[t_{i},t_{i+1}]:$

	$\displaystyle a_{0}$	$\displaystyle=g_{i},b_{0}=-(11g_{i}-18g_{i+1}+9g_{i+2}-2g_{i+3})/6h$
	$\displaystyle c_{0}$	$\displaystyle=(2g_{i}-5g_{i+1}+4g_{i+2}-g_{i+3})/2h^{2},d_{0}=-(g_{i}-3g_{i+1}+3g_{i+2}-g_{i+3})/6h^{3}$

For $\tilde{s}^{y}_{t_{i},1}(s),s\in[t_{i+1},t_{i+2}]:$

	$\displaystyle a_{1}$	$\displaystyle=g_{i+1},b_{1}=-(2g_{i}+3g_{i+1}-6g_{i+2}+g_{i+3})/6h$
	$\displaystyle c_{1}$	$\displaystyle=(g_{i}-2g_{i+1}+g_{i+2})/2h^{2},d_{1}=-(g_{i}-3g_{i+1}+3g_{i+2}-g_{i+3})/6h^{3}$

For $\tilde{s}^{y}_{t_{i},2}(s),s\in[t_{i+2},t_{i+3}]:$

	$\displaystyle a_{2}$	$\displaystyle=g_{i+2},b_{2}=(g_{i}-6g_{i+1}+3g_{i+2}+2g_{i+3})/6h$
	$\displaystyle c_{2}$	$\displaystyle=(g_{i}-2g_{i+1}+g_{i+3})/2h^{2},d_{2}=-(g_{i}-3g_{i+1}+3g_{i+2}-g_{i+3})/6h^{3}$

Thus, (16) can be rewritten as

	$\displaystyle Y_{i}$	$\displaystyle=E_{i}[Y_{i+K_{y}}]+\frac{3h}{8}g_{i}+\frac{9h}{8}g_{i+1}+\frac{9h}{8}g_{i+2}+\frac{3h}{8}g_{i+2}+R_{y}^{i}$		(45)
		$\displaystyle:=E_{i}[Y_{i+K_{y}}]+hK_{y}\sum_{j=0}^{Ky}\gamma_{K_{y},j}^{K_{y}}E_{i}[f(t_{i+j},Y_{i+j},Z_{i+j})]+R_{y}^{i},$		(46)

where $\gamma_{K_{y},0}^{K_{y}}=\gamma_{K_{y},3}^{K_{y}}=\frac{1}{8},\gamma_{K_{y},1}^{K_{y}}=\gamma_{K_{y},2}^{K_{y}}=\frac{3}{8}.$

In an analogous way we can also find coefficients for $K_{y}\geq 3,$ and report them for $1\leq K_{y}\leq 6$ in Table 1.

$K_{y}$	$\gamma_{K_{y},j}^{K_{y}}$
	$j=0$	$j=1$	$j=2$	$j=3$	$j=4$	$j=5$	$j=6$
1	$\frac{1}{2}$	$\frac{1}{2}$
2 (Second Order Polynomial)	$\frac{1}{6}$	$\frac{2}{3}$	$\frac{1}{6}$
2 (Natural Cubic Spline )	$\frac{3}{16}$	$\frac{5}{8}$	$\frac{3}{16}$
3	$\frac{1}{8}$	$\frac{3}{8}$	$\frac{3}{8}$	$\frac{1}{8}$
4	$\frac{1}{12}$	$\frac{1}{3}$	$\frac{1}{6}$	$\frac{1}{3}$	$\frac{1}{12}$
5	$\frac{41}{600}$	$\frac{19}{75}$	$\frac{107}{600}$	$\frac{107}{600}$	$\frac{19}{75}$	$\frac{41}{600}$
6	$\frac{19}{336}$	$\frac{3}{14}$	$\frac{15}{112}$	$\frac{4}{21}$	$\frac{15}{112}$	$\frac{3}{14}$	$\frac{19}{336}$

Table 1: The coefficients

[\gamma_{K_{y},j}^{K_{y}}]_{j=0}^{K_{y}}

for

K_{y}=1,2,\cdots,6.

We substitute $l=1$ into (26) and thus obtain

\begin{split}0=hE_{i}[Z_{i+1}]&+\sum_{j=0}^{Kz-1}\left[a^{z_{1}}_{j}h+\frac{b^{z_{1}}_{j}h^{2}}{2}+\frac{c^{z_{1}}_{j}h^{3}}{3}+\frac{d^{z_{1}}_{j}h^{4}}{4}\right]\\ &-\sum_{j=0}^{Kz-1}\left[a^{z_{2}}_{j}h+\frac{b^{z_{2}}_{j}h^{2}}{2}+\frac{c^{z_{2}}_{j}h^{3}}{3}+\frac{d^{z_{2}}_{j}h^{4}}{4}\right]+R_{z}^{i}.\end{split}

(47)

Note that both the sum terms in the latter equation have the same structure, they will have the same coefficients. We use $g_{i+j}$ for $E_{i}[f(t_{i+j},Y_{i+j},Z_{i+j})\Delta W_{i+j}]$ and $\tilde{g}_{i+j}$ for $E_{i}[Z_{i+j}]$ for $j=0,1,\cdots,K_{z}.$ Similar to the way of calculating the coefficients for the reference equation of $Y_{i},$ in the following we calculate the coefficients for (47).

•

$K_{z}=1:$ we construct straight lines $a^{z_{1}}_{0}=g_{i},b^{z_{1}}_{0}=\frac{g_{i+1}-g_{i}}{h},c^{z_{1}}_{0}=0,d^{z_{1}}_{0}=0$ and $a^{z_{2}}_{0}=\tilde{g}_{i},b^{z_{2}}_{0}=\frac{\tilde{g}_{i+1}-\tilde{g}_{i}}{h},c^{z_{2}}_{0}=0,d^{z_{2}}_{0}=0$ Now, we can rewrite (47) as

		$\displaystyle 0=hE_{i}[Z_{i+1}]+\frac{h}{2}g_{i}+\frac{h}{2}g_{i+1}-\frac{h}{2}\tilde{g}_{i}-\frac{h}{2}\tilde{g}_{i+1}+R_{z}^{i}$		(48)
		$\displaystyle:=hE_{i}[Z_{i+1}]+h\sum_{j=0}^{Kz}\gamma_{K_{z},j}^{1}E_{i}[f(t_{i+j},Y_{i+j},Z_{i+j})\Delta W_{i+j}]-h\sum_{j=0}^{Kz}\gamma_{K_{z},j}^{1}E_{i}[Z_{i+j}]+R_{z}^{i},$		(49)

where $\gamma_{K_{z},0}^{1}=\gamma_{K_{z},1}^{1}=\frac{1}{2}.$

•

$K_{z}=2:$ due to $l=1$ we only need to consider the interval $[t_{i},t_{i+1}].$
Using natural cubic splines: $\tilde{s}^{z_{1}}_{t_{i},0}(s),\tilde{s}^{z_{2}}_{t_{i},0}(s),s\in[t_{i},t_{i+1}]:$

	$\displaystyle a^{z_{1}}_{0}$	$\displaystyle=g_{i},b^{z_{1}}_{0}=-(5g_{i}-6g_{i+1}+g_{i+2})/4h,c^{z_{1}}_{0}=0,d^{z_{1}}_{0}=(g_{i}-2g_{i+1}+g_{i+2})/4h^{3}$
	$\displaystyle a^{z_{2}}_{0}$	$\displaystyle=\tilde{g}_{i},b^{z_{2}}_{0}=-(5\tilde{g}_{i}-6\tilde{g}_{i+1}+\tilde{g}_{i+2})/4h,c^{z_{2}}_{0}=0,d^{z_{2}}_{0}=(\tilde{g}_{i}-2\tilde{g}_{i+1}+\tilde{g}_{i+2})/4h^{3}$

Thus, (47) can be rewritten as

		$\displaystyle 0=hE_{i}[Z_{i+1}]+\frac{7h}{16}g_{i}+\frac{10h}{16}g_{i+1}-\frac{h}{16}g_{i+2}-(\frac{7h}{16}\tilde{g}_{i}+\frac{10h}{16}\tilde{g}_{i+1}-\frac{h}{16}\tilde{g}_{i+2})+R_{z}^{i}$		(50)
		$\displaystyle:=hE_{i}[Z_{i+1}]+h\sum_{j=0}^{Kz}\gamma_{K_{z},j}^{1}E_{i}[f(t_{i+j},Y_{i+j},Z_{i+j})\Delta W_{i+j}]-h\sum_{j=0}^{Kz}\gamma_{K_{z},j}^{1}E_{i}[Z_{i+j}]+R_{z}^{i},$		(51)

where $\gamma_{K_{z},0}^{1}=\frac{7}{16},\gamma_{K_{z},1}^{1}=\frac{5}{8},\gamma_{K_{z},2}^{1}=-\frac{1}{16}.$

Using the second order polynomials we obtain

\begin{split}p_{i}(s)&=g_{i}(s-t_{i})-\left(\frac{3}{2}g_{i}-2g_{i+1}+\frac{1}{2}g_{i+2}\right)(s-t_{i})/h\\ &+\left(\frac{1}{2}g_{i}-g_{i+1}+\frac{1}{2}g_{i+2}\right)(s-t_{i})^{2}/h^{2}\end{split}

(52)

\begin{split}\tilde{p}_{i}(s)&=\tilde{g}_{i}(s-t_{i})-\left(\frac{3}{2}\tilde{g}_{i}-2\tilde{g}_{i+1}+\frac{1}{2}\tilde{g}_{i+2}\right)(s-t_{i})/h\\ &+\left(\frac{1}{2}\tilde{g}_{i}-\tilde{g}_{i+1}+\frac{1}{2}\tilde{g}_{i+2}\right)(s-t_{i})^{2}/h^{2}\end{split}

(53)

whose integrations are given by

\int_{t_{i}}^{t_{i+1}}p_{i}(s)ds=h\frac{5g_{i}+8g_{i+1}-g_{i+2}}{12}.

(54)

\int_{t_{i}}^{t_{i+1}}\tilde{p}_{i}(s)ds=h\frac{5\tilde{g}_{i}+8\tilde{g}_{i+1}-\tilde{g}_{i+2}}{12}.

(55)

By using the second order polynomial we rewrite (16) as

		$\displaystyle 0=hE_{i}[Z_{i+1}]+\frac{5h}{12}g_{i}+\frac{2h}{3}g_{i+1}-\frac{h}{12}g_{i+2}-(\frac{5h}{12}\tilde{g}_{i}+\frac{2h}{3}\tilde{g}_{i+1}-\frac{h}{12}\tilde{g}_{i+2})+R_{z}^{i}$		(56)
		$\displaystyle:=hE_{i}[Z_{i+1}]+h\sum_{j=0}^{Kz}\gamma_{K_{z},j}^{1}E_{i}[f(t_{i+j},Y_{i+j},Z_{i+j})\Delta W_{i+j}]-h\sum_{j=0}^{Kz}\gamma_{K_{z},j}^{1}E_{i}[Z_{i+j}]+R_{z}^{i},$		(57)

where $\gamma_{K_{z},0}^{1}=\frac{5}{12},\gamma_{K_{z},1}^{1}=\frac{2}{3},\gamma_{K_{z},2}^{1}=-\frac{1}{12}.$

•

$K_{z}=3:$ for $K_{z}\geq 3$ we will use the Not-a-knot cubic spline.

For $\tilde{s}^{z_{1}}_{t_{i},0}(s),\tilde{s}^{z_{2}}_{t_{i},0}(s),s\in[t_{i},t_{i+1}]:$

	$\displaystyle a^{z_{1}}_{0}$	$\displaystyle=g_{i},b^{z_{1}}_{0}=-(11g_{i}-18g_{i+1}+9g_{i+2}-2g_{i+3})/6h$
	$\displaystyle c^{z_{1}}_{0}$	$\displaystyle=(2g_{i}-5g_{i+1}+4g_{i+2}-g_{i+3})/2h^{2},d^{z_{1}}_{0}=-(g_{i}-3g_{i+1}+3g_{i+2}-g_{i+3})/6h^{3}$

	$\displaystyle a^{z_{2}}_{0}$	$\displaystyle=\tilde{g}_{i},b^{z_{2}}_{0}=-(11\tilde{g}_{i}-18\tilde{g}_{i+1}+9\tilde{g}_{i+2}-2\tilde{g}_{i+3})/6h$
	$\displaystyle c^{z_{2}}_{0}$	$\displaystyle=(2\tilde{g}_{i}-5\tilde{g}_{i+1}+4\tilde{g}_{i+2}-\tilde{g}_{i+3})/2h^{2},d^{z_{2}}_{0}=-(\tilde{g}_{i}-3\tilde{g}_{i+1}+3\tilde{g}_{i+2}-g_{i+3})/6h^{3}$

Thus, (47) can be rewritten as

		$\displaystyle 0=hE_{i}[Z_{i+1}]+\frac{3h}{8}g_{i}+\frac{19h}{24}g_{i+1}-\frac{5h}{24}g_{i+2}+\frac{h}{24}g_{i+3}$
		$\displaystyle-(\frac{3h}{8}\tilde{g}_{i}+\frac{19h}{24}\tilde{g}_{i+1}-\frac{5h}{24}\tilde{g}_{i+2}+\frac{h}{24}\tilde{g}_{i+3})+R_{z}^{i}$		(58)
		$\displaystyle:=hE_{i}[Z_{i+1}]+h\sum_{j=0}^{Kz}\gamma_{K_{z},j}^{1}E_{i}[f(t_{i+j},Y_{i+j},Z_{i+j})\Delta W_{i+j}]-h\sum_{j=0}^{Kz}\gamma_{K_{z},j}^{1}E_{i}[Z_{i+j}]+R_{z}^{i},$		(59)

where $\gamma_{K_{z},0}^{1}=\frac{3}{8},\gamma_{K_{z},1}^{1}=\frac{19}{24},\gamma_{K_{z},2}^{1}=-\frac{5}{24},\gamma_{K_{z},3}^{1}=\frac{1}{24}.$

The coefficients for $1\leq K_{z}\leq 6$ are reported in Table 2.

$K_{z}$	$\gamma_{K,j}^{1}$
	$j=0$	$j=1$	$j=2$	$j=3$	$j=4$	$j=5$	$j=6$
1	$\frac{1}{2}$	$\frac{1}{2}$
2 (Second Order Polynomial)	$\frac{5}{12}$	$\frac{2}{3}$	$-\frac{1}{12}$
2 (Natural Cubic Spline )	$\frac{7}{16}$	$\frac{5}{8}$	$-\frac{1}{16}$
3	$\frac{3}{8}$	$\frac{19}{24}$	$-\frac{5}{24}$	$\frac{1}{24}$
4	$\frac{35}{96}$	$\frac{5}{6}$	$-\frac{13}{48}$	$\frac{1}{12}$	$-\frac{1}{96}$
5	$\frac{131}{360}$	$\frac{151}{180}$	$-\frac{103}{360}$	$\frac{37}{360}$	$-\frac{1}{45}$	$\frac{1}{360}$
6	$\frac{163}{448}$	$\frac{47}{56}$	$-\frac{129}{448}$	$\frac{3}{28}$	$-\frac{37}{1344}$	$\frac{1}{168}$	$-\frac{1}{1344}$

Table 2: The coefficients

[\gamma_{K_{z},j}^{1}]_{j=0}^{K_{z}}

for

K_{z}=1,2,\cdots,6.

Note that $\Delta W_{t_{i}}=0$ and $E_{i}[Z_{i}]=Z_{i},$ based on the calculations above we can obtain the reference equations of the BSDEs as

	$\displaystyle Y_{i}$	$\displaystyle=E_{i}[Y_{i+K_{y}}]+hK_{y}\sum_{j=0}^{Ky}\gamma_{K_{y},j}^{K_{y}}E_{i}[f(t_{i+j},Y_{i+j},Z_{i+j})]+R_{y}^{i},$		(60)
	$\displaystyle Z_{i}$	$\displaystyle=\left(E_{i}[Z_{i+1}]+\sum_{j=1}^{Kz}\gamma_{K_{z},j}^{1}E_{i}[f(t_{i+j},Y_{i+j},Z_{i+j})\Delta W_{i+j}]-\sum_{j=1}^{Kz}\gamma_{K_{z},j}^{1}E_{i}[Z_{i+j}]\right)/\gamma_{K_{z},0}^{1}+R_{z}^{i},$		(61)

where $Y_{i}=\left(Y_{i}^{1},Y_{i}^{2},\cdots,Y_{i}^{m}\right)^{\top}$ , $Z_{i}=\left(Z^{\tilde{m},\tilde{d}}_{i}\right)_{m\times d},$ $\Delta W_{i+j}=(W_{i+j}^{1},W_{i+j}^{2},\cdots,W_{i+j}^{d})^{\top}-(W_{i}^{1},W_{i}^{2},\cdots,W_{i}^{d})^{\top},$ $R_{y}^{i}=\left(R_{y}^{i,1},R_{y}^{i,2},\cdots,R_{y}^{i,{m}}\right)^{\top}$ and $R_{z}^{i}=\left(R_{z}^{i,{\tilde{m},\tilde{d}}}\right)_{m\times d}.$ It is easy to see that (60) is implicit, and (61) is always explicit for solving $Z_{i}.$ One can show that estimates for the local error terms $R^{i}_{y}$ and $R^{i}_{z}$ (componentwise in (31) and (35)) are given by

|R_{y}^{i}|=\mathcal{O}(h^{5}),\quad|R_{z}^{i}|=\mathcal{O}(h^{5})

(62)

provided that the generator function $f$ and the terminal function $g$ are smooth. It is worth noting that $R_{z}^{i}$ will be divided by $h$ for solving $Z_{i},$ see e.g., (59), one might set $K_{z}=K_{y}+1$ in order to balance the local truncation errors.

4 A stable multistep discretization scheme

In this Section we present a stable multistep scheme fully discrete in time and space.

4.1 The Semi-discretization in time

We denote $Y^{i}=\left(Y^{1,i},Y^{2,i},\cdots,Y^{m,i}\right)^{\top}$ and $Z^{i}=\left(Z^{\tilde{m},\tilde{d},i}\right)_{m\times d}$ as the approximations to $Y_{i}$ and $Z_{i},$ namely at the time $t_{i}$ in the reference equations, respectively. Furthermore, we have $W_{i}=(W_{i}^{1},W_{i}^{2},\cdots,W_{i}^{d})^{\top},$ whereas all Brownian motions are independent with each other. Since $Z_{i}$ is needed for computing $Y_{i}$ in our scheme, we thus need to consider the larger step size between $K_{y}$ and $K_{z}.$ Therefore, we define the number of time steps as $K=\max\left\{K_{y},K_{z}\right\}.$ Suppose that the random variables $Y^{N_{T}-j}$ and $Z^{N_{T}-j}$ are given for $j=0,1,\cdots,K-1,$ then $Y^{i}$ and $Z^{i}$ can be found for $i=N_{T}-K,\cdots,0$ by

	$\displaystyle Y^{i}$	$\displaystyle=E_{i}[Y^{i+K_{y}}]+hK_{y}\sum_{j=0}^{Ky}\gamma_{K_{y},j}^{K_{y}}E_{i}[f(t_{i+j},Y^{i+j},Z^{i+j})],$		(63)
	$\displaystyle Z^{i}$	$\displaystyle=\left(E_{i}[Z^{i+1}]+\sum_{j=1}^{Kz}\gamma_{K_{z},j}^{1}E_{i}[f(t_{i+j},Y^{i+j},Z^{i+j})\Delta W^{\top}_{i+j}]-\sum_{j=1}^{Kz}\gamma_{K_{z},j}^{1}E_{i}[Z^{i+j}]\right)/\gamma_{K_{z},0}^{1},$		(64)

We follow the methodologies used in [Zhao et al., 2010] to check the stability. We set the generator function $f=0$ and take the expectation $E[\cdot]$ on both sides of (63)

E[Y^{i}]=E[Y^{i+k}].

(65)

Note that we have set $k=K_{y}$ in (63). We need to recall $k$ in (65) for a general stability analysis. (65) indicates that reference equation of $Y_{i}$ is stable for any integers $1\leq k\leq K_{y}\leq N_{T}.$ Furthermore, in (63) where $k=K_{y},$ we have checked that $\sum_{j=0}^{K_{y}}\gamma_{K_{y},j}^{K_{y}}=1$ for $1\leq K_{y}\leq N_{T}.$

In a similar way to above, (61) can be reformulated as

0=E[Z^{i+l}]-\sum_{j=1}^{Kz}\gamma_{K_{z},j}^{l}E[Z^{i+j}],

(66)

where $l$ is recalled substituting $1$ in (64). We see that (66) is a difference equation of $Z^{i},$ the characteristic polynomial of the backward difference equation (66) reads

p^{l}_{K_{z}}(\lambda)=\lambda^{K_{z}-l}-\sum_{j=1}^{Kz}\gamma_{K_{z},j}^{l}\lambda^{K_{z}-j}.

(67)

In order to have a stable reference equation of $Z^{i},$ the roots of (67) must satisfy the following condition:

•

The roots must be in the closed unit disc and the ones on the unit circle must be simple.

The values of $\gamma_{K_{z},j}^{1}$ have been given for $K_{z}=1,\cdots,6$ in Table 2. In the same way as we obtained those values one can calculate the values of $\gamma_{K_{z},l}^{j}$ for $1<l\leq K_{z}\leq N_{T}$ and obtain the corresponding roots of (67), see Table 3.

$K_{z}$	$l$	Roots $\lambda^{l}_{K_{z},j}$
		$j=1$	$j=2$	$j=3$	$j=4$	$j=5$	$j=6$
1	1	1
2	1	1	$-\frac{1}{5}\>\>\>(-\frac{1}{7}\text{ natural CS})$
2	2	1	$-5\>\>\>(-4.3333\text{ natural CS})$
3	1	1	$\frac{\sqrt{13}}{9}-\frac{2}{9}$	$-\frac{\sqrt{13}}{9}-\frac{2}{9}$
	2	1	$0$	$-5$
	3	1	$\sqrt{3}i-2$	$-\sqrt{3}i-2$
4	1	1	$-0.82662$	$0.14188-0.12014i$	$0.14188+0.12014i$
	2	1	$0$	$0$	$-5$
	3	1	$-0.01244$	$-2.31196+1.40033i$	$-2.31196-1.40033i$
	4	1	$-3.93114$	$-0.53442+1.5851i$	$-0.53442-1.5851i$
5	1	1	$-0.89193$	$0.20080$	$0.06693+0.19529i$	$0.06693-0.19529i$
	2	1	$0$	$0$	$0$	$-5$
	3	1	$-0.07259$	$0.04667$	$-2.34069-1.31158i$	$-2.34069+1.31158i$
	4	1	$-3.64370$	$-0.00620$	$-0.57668-1.60195i$	$-0.57668+1.60195i$
	5	1	$-2.45215+0.06565i$	$-2.45215-0.06565i$	$-0.09849-1.50203i$	$-0.09849+1.50203i$
6	1	1	$-0.91034$	$-0.01033-0.22612i$	$-0.01033+0.22612i$	$0.18636-0.09543i$	$0.18636+0.09543i$
	2	1	$0$	$0$	$0$	$0$	$-5$
	3	1	$-0.13432$	$-2.34031+1.29934i$	$-2.34031-1.29934i$	$0.05126+0.06452i$	$0.05126-0.06452i$
	4	1	$-3.61188$	$-0.04794$	$0.03504$	$-0.58234-1.59752i$	$-0.58234+1.59752i$
	5	1	$-3.00560$	$-1.94659$	$-0.00538$	$0.09695-1.51077i$	$0.09695+1.51077i$
	6	1	$-3.38909$	$-1.14732+1.07617i$	$-1.14732-1.07617i$	$0.44714+1.33772i$	$0.44714-1.33772i$

Table 3: The roots of (67) for

K_{z}=1,2,\cdots,6

and

l=1,\cdots,K_{z}

Note that, for $K_{y}=1,2,3$ and $K_{z}=1,2,3,$ our reference equations (with second order polynomial for $K=2$ ) coincide with the reference equations proposed in [Zhao et al., 2010], where Lagrange interpolating polynomials are employed. However, in [Zhao et al., 2010], the reference equation of $Y^{i}$ is stable only when $K_{y}=1,2,3,4,5,6,7,9;$ and the reference equation of $Z^{i}$ is stable only when $K_{z}=1,2,3.$ As mentioned already, our both reference equations are generally stable, namely for all $K_{y}\geq 1$ and $K_{z}\geq 1.$ This is to say that our method allows for considering more multi-time levels.

4.2 Error analysis

Due to the nested conditional expectations we still are confronted with a problem to perform error analysis for the proposed multi-step scheme. In [Zhao et al., 2010], the authors have finished some error analysis for the multi-step semidiscrete scheme in one-dimensional case using the Lagrange interpolating polynomials under several assumptions. In this section, we adopt their results to our multi-step scheme. Throughout this section we assume that the functions $f$ and $g$ are bounded and smooth enough with bounded derivatives for a uniquely existing solution. Furthermore, suppose that $f$ does not involve the variable $Z_{t},$ i.e.,

Y_{t}=\xi+\int_{t}^{T}f(s,Y_{s})\,ds-\int_{t}^{T}Z_{s}\,dW_{s},

(68)

for which the reference equation read

	$\displaystyle Y_{i}$	$\displaystyle=E_{i}[Y_{i+K_{y}}]+hK_{y}\sum_{j=0}^{Ky}\gamma_{K_{y},j}^{K_{y}}E_{i}[f(t_{i+j},Y_{i+j})]+R_{y}^{i},$		(69)
	$\displaystyle Z_{i}$	$\displaystyle=\left(E_{i}[Z_{i+1}]+\sum_{j=1}^{Kz}\gamma_{K_{z},j}^{1}E_{i}[f(t_{i+j},Y_{i+j})\Delta W_{i+j}]-\sum_{j=1}^{Kz}\gamma_{K_{z},j}^{1}E_{i}[Z_{i+j}]\right)/\gamma_{K_{z},0}^{1}+R_{z}^{i}/h,$		(70)

where the local truncation errors $R_{y}^{i}$ and $R_{z}^{i}$ are defined in (11) and (24). And the corresponding multi-step scheme for $Y^{i}$ and $Z^{i}$ can be immediately written down from (63) and (64).

Lemma 4.1.

The local estimates of the local truncation errors in (69) and (70) satisfy

|R_{y}^{i}|\leq Ch^{\min\{K_{y}+2,\,5\}}\quad|R_{z}^{i}|\leq Ch^{\min\{K_{z}+2,\,5\}},

where $C>0$ is a constant depending on $T,f,g$ and the derivatives of $f,g.$

The proof can be done directly by combining the proof of Lemma 3.2 in [Zhao et al., 2009] and the fact that not-a-knot cubic spline is fourth-order accurate.

Theorem 4.2.

Suppose that the initial values satisfy

\left\{\begin{array}[]{l}\max_{N_{T}-K_{y}<i\leq N_{T}}E\left[\left|Y_{i}-Y^{i}\right|\right]=\mathcal{O}(h^{K_{y}+1}),~\mbox{for}~K_{y}=1,2,3\\ \max_{N_{T}-K_{y}<i\leq N_{T}}E\left[\left|Y_{i}-Y^{i}\right|\right]=\mathcal{O}(h^{4}),~\mbox{for}~K_{y}>3\end{array}\right.

for sufficiently small time step $h$ it can be shown that

\sup_{0\leq i\leq N_{T}}E\left[\left|Y_{i}-Y^{i}\right|\right]\leq Ch^{\min\{K_{y}+1,\,4\}},

(71)

where $C>0$ is a constant depending on $T,f,g$ and the derivatives of $f,g.$

The proof can be done directly by combining the proof of Theorem 1. in [Zhao et al., 2010] and the fact that not-a-knot cubic spline is fourth-order accurate.

Theorem 4.3.

Suppose that the initial values satisfy

\left\{\begin{array}[]{l}\max_{N_{T}-K_{z}<i\leq N_{T}}E\left[\left|Z_{i}-Z^{i}\right|\right]=\mathcal{O}(h^{K_{z}}),~\mbox{for}~K_{z}=1,2,3\\ \max_{N_{T}-K_{z}<i\leq N_{T}}E\left[\left|Z_{i}-Z^{i}\right|\right]=\mathcal{O}(h^{3})~\mbox{for}~K_{z}>3\end{array}\right.

and the condition on the initial values in Theorem 4.2 is fulfilled. For sufficiently small time step $h$ it can be shown that

\sup_{0\leq i\leq N_{T}}E\left[\left|Z_{i}-Z^{i}\right|\right]\leq Ch^{\min(K_{y}+1,\,K_{z},\,3)},

where $C>0$ is a constant depending on $T,f,g$ and the derivatives of $f,g.$

The proof can be done directly by combining the proof of Theorem 2. in [Zhao et al., 2010] and the fact that not-a-knot cubic spline is fourth-order accurate.

4.3 The fully discretized scheme

We have checked that (63) and (64) are stable in the time direction. To solve $(Y^{i},Z^{i})$ numerically, next we consider the space discretization. We define firstly the partion of the one-dimensional $(\tilde{d}=d=1)$ real axis as

\mathcal{R}^{\tilde{d}}=\left\{x_{\gamma}^{\tilde{d}}|x_{\gamma}^{\tilde{d}}\in\mathbb{R},\gamma\in\mathbb{Z},x_{\gamma}^{\tilde{d}}<x_{\gamma+1}^{\tilde{d}},\lim_{i\to+\infty}x_{\gamma}^{\tilde{d}}=+\infty,\lim_{i\to-\infty}x_{\gamma}^{\tilde{d}}=-\infty\right\}.

(72)

Thus, the partition of $d$ -dimensional space $\mathcal{R}^{d}$ reads

\mathcal{R}^{\tilde{d}}=\mathcal{R}^{1}\times\cdots\times\mathcal{R}^{\tilde{d}}\times\cdots\times\mathcal{R}^{d},

(73)

where $\tilde{d}=1,2,\cdots,d.$ For simplicity of notation we will use $x_{\Gamma}=(x^{1}_{\gamma_{1}},x^{2}_{\gamma_{2}},\cdots,x^{d}_{\gamma_{d}})^{\top}$ for $\Gamma=(\gamma_{1},\gamma_{2},\cdots,\gamma_{d})\in\mathbb{Z}^{d}.$ We use $y^{N_{T}-\lambda}_{\Gamma}$ and $z^{N_{T}-\lambda}_{\Gamma}$ to denote the values of random variables $Y^{N_{T}-\lambda}$ and $Z^{N_{T}-\lambda}$ at the points $x_{\Gamma}.$ Given these values for $\lambda=0,1,\cdots,K-1,$ we need to find $(y^{i}_{\Gamma},z^{i}_{\Gamma}),i=N_{T}-K,\cdots,0$ such that

	$\displaystyle y^{i}_{\Gamma}$	$\displaystyle=E^{x_{\Gamma}}_{i}[\hat{Y}^{i+K_{y}}]+hK_{y}\sum_{j=0}^{Ky}\gamma_{K_{y},j}^{K_{y}}E^{x_{\Gamma}}_{i}[f(t_{i+j},\hat{Y}^{i+j},\hat{Z}^{i+j})],$		(74)
	$\displaystyle z^{i}_{\Gamma}$	$\displaystyle=\left(E^{x_{\Gamma}}_{i}[\hat{Z}^{i+1}]+\sum_{j=1}^{Kz}\gamma_{K_{z},j}^{1}E^{x_{\Gamma}}_{i}[f(t_{i+j},\hat{Y}^{i+j},\hat{Z}^{i+j})\Delta W^{\top}_{i+j}]-\sum_{j=1}^{Kz}\gamma_{K_{z},j}^{1}E^{x_{\Gamma}}_{i}[\hat{Z}^{i+j}]\right)/\gamma_{K_{z},0}^{1},$		(75)

where $E^{x_{\Gamma}}_{i}[\cdot]$ denotes the conditional expectation under the $\sigma$ -field $\mathcal{F}_{t}^{x_{\Gamma}}$ generated by $\{W_{i}=x_{\Gamma}\}.$ Correspondingly, $\hat{Y}^{i+j}$ and $\hat{Z}^{i+j}$ denote the functions of increment of Brownian motion $Y^{i+j}(\Delta W_{i})$ and $Z^{i+j}(\Delta W_{i})$ with the fixed $\{W_{i}=x_{\Gamma}\}.$

To approximate the conditional expectations in (74) and (75) we employ the Gauss-Hermite quadrature formula. For example, we compute $E^{x_{\Gamma}}_{i}[\hat{Y}^{i+K_{y}}]$ as

$\displaystyle E^{x_{\Gamma}}_{i}[\hat{Y}^{i+K_{y}}]$	$\displaystyle=\frac{1}{(2K_{y}\pi h)^{d/2}}\int_{\mathbb{R}^{d}}\hat{Y}^{i+K_{y}}(s)\exp\left(-\frac{(s-x)^{\top}(s-x)}{2K_{y}h}\right)\,ds$	(76)
	$\displaystyle\approx\frac{1}{(2K_{y}\pi h)^{d/2}}\int_{\mathbb{R}^{d}}\hat{y}^{i+K_{y}}(s)\exp\left(-\frac{(s-x)^{\top}(s-x)}{2K_{y}h}\right)\,ds$	(77)
	$\displaystyle\approx\frac{1}{\pi^{\frac{d}{2}}}\sum_{\Lambda=1}^{L}\omega_{\Lambda}\hat{y}^{i+K_{y}}(x_{\Gamma}+\sqrt{2K_{y}h}a_{\Lambda})$	(78)
	$\displaystyle:=\hat{E}^{x_{\Gamma}}_{i}[\hat{Y}^{i+K_{y}}],$	(79)

where $\hat{y}^{i+K_{y}}(s)$ are interpolating values at the space points $s$ based on $y_{\Gamma}^{i+K_{y}}$ at a finite number of the space grid points $x_{\Gamma}$ near $s,\,\Lambda=(\lambda_{1},\lambda_{2},\cdots,\lambda_{d}),\,\omega_{\Lambda}=\prod_{\tilde{d}=1}^{d}\omega_{\lambda_{\tilde{d}}},\,a_{\Lambda}=(a_{\lambda_{1}},a_{\lambda_{2}},\cdots,a_{\lambda_{d}}),\,\sum_{\Lambda=1}^{L}=\sum_{\lambda_{1}=1,\cdots,\lambda_{d}=1}^{L,\cdots,L}.$ For the weights $\omega_{\Lambda}$ and the roots $a_{\Lambda}$ we refer to e.g., [Abramowitz and Stegun, 1972]. The approximations of the other conditional expectations in (74) and (75) can be done similarly. Finally, by considering these approximations we rewrite (74) and (75) as

	$\displaystyle y^{i}_{\Gamma}$	$\displaystyle=\hat{E}^{x_{\Gamma}}_{i}[\hat{Y}^{i+K_{y}}]+hK_{y}\sum_{j=0}^{Ky}\gamma_{K_{y},j}^{K_{y}}\hat{E}^{x_{\Gamma}}_{i}[f(t_{i+j},\hat{Y}^{i+j},\hat{Z}^{i+j})],$		(80)
	$\displaystyle z^{i}_{\Gamma}$	$\displaystyle=\left(\hat{E}^{x_{\Gamma}}_{i}[\hat{Z}^{i+1}]+\sum_{j=1}^{Kz}\gamma_{K_{z},j}^{1}\hat{E}^{x_{\Gamma}}_{i}[f(t_{i+j},\hat{Y}^{i+j},\hat{Z}^{i+j})\Delta W^{\top}_{i+j}]-\sum_{j=1}^{Kz}\gamma_{K_{z},j}^{1}\hat{E}^{x_{\Gamma}}_{i}[\hat{Z}^{i+j}]\right)/\gamma_{K_{z},0}^{1}.$		(81)

We observe that the computations at each space grid point are independent, which can be thus parallelized. Usually, only the values of $y^{N_{T}}_{\Gamma}$ and $z^{N_{T}}_{\Gamma}$ are known because of the terminal condition. However, for a $K$ -step scheme we need to know the support values of $y^{N_{T}-j}_{\Gamma}$ and $z^{N_{T}-j}_{\Gamma},j=0,\cdots,K-1.$ One can use the following two ways to deal with this problem: before running the multi-step scheme, we choose a quite smaller $h$ and run one-step scheme until $N_{T}-K;$ Alternatively, one can prepare these initial values “iteratively”, namely we compute $y^{N_{T}-1}_{\Gamma}$ and $z^{N_{T}-1}_{\Gamma}$ based on $y^{N_{T}}_{\Gamma}$ and $z^{N_{T}}_{\Gamma}$ with $K=1,$ and the compute $y^{N_{T}-2}_{\Gamma}$ and $z^{N_{T}-2}_{\Gamma}$ based on $y^{N_{T}}_{\Gamma},y^{N_{T}-1}_{\Gamma},z^{N_{T}}_{\Gamma},z^{N_{T}-1}_{\Gamma}$ with $K=2$ and so on. Notice that we are faced with a computational complexity problem for solving high-dimensional problem, since the number of the Gauss-Hermite quadrature points grows exponentially with the dimension $d.$

5 Numerical experiments

In this section we use some numerical examples to show the high effectiveness and accuracy of our scheme for solving the BSDEs. We choose the truncated domain for the Brownian motion to be $[-8,8]^{d},$ and the degree of the Hermite polynomial (see $L$ in (78) )to be $8.$ Note that, for $L=8,$ the quadrature error is so small that it cannot affect the convergence rate. We use the Newton-Raphson method to implicitly solve (80). For the interpolation method we apply cubic spline interpolation which is a fourth-order accurate, namely $(\Delta x)^{4}.$ In order to be able to estimate the convergence rate in time, we adjust the space step size $\Delta x$ according to the time step size $h$ such that $(\Delta x)^{4}=(h)^{q+1}$ with $q=\min\{K_{y}+1,K_{z}\}.$ In the general case (the generator $f$ depends on both $Y_{t}$ and $Z_{t}$ ), from Theorem 4.3 we know that $q$ is only limited to 3, since not-a-knot cubic spline is maximal fourth-rate accurate. This is to say that we always take $q=3$ when $\min\{K_{y}+1,K_{z}\}\geq 3.$ However, when the generator $f$ does not involve the component $Z_{t},$ the approximation for $Y_{t}$ can reach fourth-order accurate, see Theorem 4.2. For this case, $q$ is allowed to be $4$ when $\min\{K_{y}+1,K_{z}\}\geq 4.$

Generally, only $Y_{N_{T}}$ and $Z_{N_{T}}$ are known analytically. However, as mentioned before, for a $K$ -step scheme we need to know $y^{N_{T}-j}_{\Gamma}$ and $z^{N_{T}-j}_{\Gamma},j=1,\cdots,K-1$ as initial values as well. To obtain these initial values, we start with $K=1$ and choose a extremely small time step size $h.$ Because the largest number of steps in our experiments is $K=6,$ we start thus with $N_{T}=8.$ In our computation we have used parallel computing using Python’s multiprocessing module. Note that a GPU-based parallelism will be much more cost-effective, which is left as a future work.

As mentioned before, our algorithm coincides with the algorithm proposed in [Zhao et al., 2010] for $K_{y}=1,2,3$ and $K_{z}=1,2,3.$ In [Zhao et al., 2010], the authors have compared the multi-step scheme to the implicit Euler scheme [Zhao et al., 2009] and the $\theta$ -scheme [Zhao et al., 2006]. For these implicit Euler scheme and $\theta$ -scheme, they have considered both the Monte-Carlo method and the Gaussian quadrature for approximating the conditional expectations. Therefore, we will not do any comparison with other methods, for this we refer [Zhao et al., 2010]. In our numerical examples we will demonstrate higher effectiveness and accuracy of our scheme, which allows for more than $3$ -step scheme, namely $K>3.$

Example 1

The first example reads

\begin{cases}&-dY_{t}=-\frac{5}{8}Y_{t}\,dt-Z_{t}\,dW_{t},\\ &Y_{T}=\exp(W_{T}/2+T/2),\end{cases}

with the analytic solution

\begin{cases}&Y_{t}=\exp(W_{t}/2+t/2),\\ &Z_{t}=\exp(W_{t}/2+t/2)/2.\end{cases}

The exact solution of $(Y_{0},Z_{0})$ is thus $\left(1,\frac{1}{2}\right).$ Obviously, in this example, the generator $f$ does not depend on $Z_{t}.$ We thus choose $q=\min\{K_{y}+1,K_{z}\}<4$ and keep $q=4$ when $\min\{K_{y}+1,K_{z}\}\geq 4.$ This is to say that the value of $q$ is exactly the theoretical convergence order for the $Y$ -component solver. For the $Z$ -component, the theoretical convergence order of our scheme is $\min\{K_{y}+1,K_{z}\}$ but limited by $3$ due to Theorem 4.3. The corresponding numerical results and estimated convergence rates are reported in Table 4 and 5. For $K=1,\cdots,4,$ we have considered many combinations with the different values of $K_{y},K_{z}$ and the corresponding values of $q.$ The results of these combinations are also similar for $K\geq 5.$ Therefore, for $K=5,6$ we only report the results for $K_{y}=K_{z}=5,6$ which are sufficient to show the benefit from a higher number of multi-step.

$\|Y_{0}-y_{0}^{0}\|$
	$N_{T}=8$	$N_{T}=16$	$N_{T}=32$	$N_{T}=64$	$N_{T}=128$	CR
$K_{y}=1,K_{z}=1,q=1$	3.40e-04	8.90e-05	2.48e-05	7.37e-06	2.48e-06	1.78
$K_{y}=1,K_{z}=2,q=2$	3.19e-04	7.96e-05	2.00e-05	5.00e-06	1.25e-06	2.00
$K_{y}=2,K_{z}=1,q=1$	6.26e-06	2.81e-06	1.46e-06	7.16e-07	3.69e-07	1.01
$K_{y}=2,K_{z}=2,q=2$	8.79e-07	3.24e-07	4.57e-08	8.83e-09	4.28e-09	2.06
$K_{y}=2,K_{z}=3,q=3$	2.05e-07	1.16e-08	2.03e-09	1.38e-10	3.11e-11	3.18
$K_{y}=3,K_{z}=1,q=1$	7.33e-07	2.06e-07	1.75e-07	8.88e-08	5.67e-08	0.86
$K_{y}=3,K_{z}=2,q=2$	6.60e-07	8.09e-08	2.55e-08	8.35e-09	1.33e-09	2.12
$K_{y}=3,K_{z}=3,q=3$	2.30e-07	2.52e-08	1.79e-09	2.58e-10	2.29e-11	3.32
$K_{y}=3,K_{z}=4,q=4$	1.99e-07	1.77e-08	1.07e-09	7.05e-11	4.50e-12	3.88
$K_{y}=4,K_{z}=1,q=1$	3.23e-07	5.36e-07	2.54e-07	1.42e-07	6.64e-08	0.64
$K_{y}=4,K_{z}=2,q=2$	5.11e-07	8.37e-08	3.77e-08	1.55e-09	1.68e-09	2.23
$K_{y}=4,K_{z}=3,q=3$	1.54e-07	1.50e-08	9.64e-10	1.49e-10	8.19e-12	3.51
$K_{y}=4,K_{z}=4,q=4$	1.54e-07	9.29e-09	5.59e-10	3.40e-11	2.04e-12	4.05
$K_{y}=4,K_{z}=5,q=4$	1.54e-07	9.29e-09	5.59e-10	3.40e-11	2.04e-12	4.05
$K_{y}=5,K_{z}=5,q=4$	6.48e-08	7.06e-09	4.12e-10	2.54e-11	1.66e-12	3.86
$K_{y}=6,K_{z}=6,q=4$	6.60e-08	3.81e-09	3.21e-10	1.92e-11	1.32e-12	3.89

Table 4: Errors and convergence rates for Example 1,

T=1

$\|Z_{0}-z_{0}^{0}\|$
	$N_{T}=8$	$N_{T}=16$	$N_{T}=32$	$N_{T}=64$	$N_{T}=128$	CR
$K_{y}=1,K_{z}=1,q=1$	1.71e-02	8.52e-03	4.25e-03	2.12e-03	1.06e-03	1.00
$K_{y}=1,K_{z}=2,q=2$	8.50e-04	2.24e-04	5.76e-05	1.46e-05	3.67e-06	1.97
$K_{y}=2,K_{z}=1,q=1$	1.72e-02	8.54e-03	4.26e-03	2.12e-03	1.06e-03	1.00
$K_{y}=2,K_{z}=2,q=2$	7.89e-04	2.09e-04	5.37e-05	1.36e-05	3.42e-06	1.96
$K_{y}=2,K_{z}=3,q=3$	4.17e-05	6.02e-06	8.03e-07	1.04e-07	1.32e-08	2.91
$K_{y}=3,K_{z}=1,q=1$	1.72e-02	8.54e-03	4.26e-03	2.12e-03	1.06e-03	1.00
$K_{y}=3,K_{z}=2,q=2$	7.89e-04	2.09e-04	5.37e-05	1.36e-05	3.42e-06	1.96
$K_{y}=3,K_{z}=3,q=3$	4.16e-05	6.02e-06	8.03e-07	1.04e-07	1.32e-08	2.91
$K_{y}=3,K_{z}=4,q=4$	1.98e-05	3.24e-06	4.59e-07	6.10e-08	7.84e-09	2.83
$K_{y}=4,K_{z}=1,q=1$	1.72e-02	8.54e-03	4.26e-03	2.12e-03	1.06e-03	1.00
$K_{y}=4,K_{z}=2,q=2$	7.89e-04	2.09e-04	5.37e-05	1.36e-05	3.42e-06	1.96
$K_{y}=4,K_{z}=3,q=3$	4.17e-05	6.02e-06	8.03e-07	1.04e-07	1.32e-08	2.91
$K_{y}=4,K_{z}=4,q=4$	1.98e-05	3.25e-06	4.60e-07	6.10e-08	7.90e-09	2.83
$K_{y}=4,K_{z}=5,q=4$	1.67e-05	3.34e-06	5.00e-07	6.77e-08	1.30e-08	2.83
$K_{y}=5,K_{z}=5,q=4$	1.67e-05	3.34e-06	4.99e-07	6.77e-08	1.10e-08	2.68
$K_{y}=6,K_{z}=6,q=4$	1.29e-05	2.93e-06	4.61e-07	6.39e-08	1.60e-10	3.81

Table 5: Errors and convergence rates for Example 1,

T=1

By a columnwise comparison we see that the approximation errors reduce mostly with the increasing number of steps, $K_{y}$ and $K_{z}.$ We have obtained $10^{-8}$ for approximating $Y_{t}$ already with $N_{T}=8,$ namely $h=\frac{1}{8}.$ The estimated convergence rates¹¹ 1 Estimated by using linear squares fitting. (CR) for both of $Y_{t}$ and $Z_{t}$ are consistent with the theoretical results explained before, if we ignore the quadrature and interpolation errors which can cause a slightly smaller estimated convergence rate. In Table 5 we even observe a better CR than the theoretical result for $K_{y}=K_{z}=6.$ We display the plots of $\log_{2}\left(|Y_{0}-y_{0}^{0}|\right)$ and $\log_{2}\left(|Z_{0}-z_{0}^{0}|\right)$ with respect to $\log_{2}(N_{T})$ in Figure 1.

(a)

Y

-component

(b)

Z

-component

Figure 1: Plots of

\log_{2}\left(|Y_{0}-y_{0}^{0}|\right)

and

\log_{2}\left(|Z_{0}-z_{0}^{0}|\right)

with respect to

\log_{2}(N_{T})

for

K=1,\cdots 6

for Example 1.

For this example, we also run our algorithm separately (without computing the $Z$ -component) for solving the $Y$ -component with smaller space step size $\Delta x$ (higher value of $q$ ). For $K_{y}\geq 4,$ we compare the numerical solutions computed with $q=4,\cdots,K_{y}+1.$

$\|Y_{0}-y_{0}^{0}\|$
	$N_{T}=8$	$N_{T}=16$	$N_{T}=32$	$N_{T}=64$	$N_{T}=128$	CR
$K_{y}=4,q=4$	1.54e-07	9.29e-09	5.59e-10	3.40e-11	2.04e-12	4.05
$K_{y}=4,q=5$	1.53e-07	8.85e-09	5.30e-10	3.23e-11	2.00e-12	4.05
$K_{y}=5,q=4$	6.48e-08	7.06e-09	4.12e-10	2.54e-11	1.66e-12	3.86
$K_{y}=5,q=5$	6.24e-08	6.73e-09	4.03e-10	2.44e-11	1.63e-12	3.86
$K_{y}=5,q=6$	6.21e-08	6.71e-09	4.02e-10	2.44e-11	1.63e-12	3.86
$K_{y}=6,q=4$	6.60e-08	3.81e-09	3.21e-10	1.92e-11	1.32e-12	3.89
$K_{y}=6,q=5$	6.53e-08	3.62e-09	3.10e-10	1.87e-11	1.25e-12	3.89
$K_{y}=6,q=6$	6.50e-08	3.62e-09	3.09e-10	1.87e-11	1.25e-12	3.89
$K_{y}=6,q=7$	6.49e-08	3.62e-09	3.09e-10	1.87e-11	1.25e-12	3.89

Table 6: Errors and convergence rates for Example 1, where

y_{0}^{0}

is separately computed for different higher values of

q

and

T=1.

The reported results in Table 6 have shown clearly that there is almost no benefit to setting $q=K_{y}+1$ when $K_{y}+1>4,$ i.e., we only need to keep $q=4$ for $K_{y}+1>4.$ We emphasise again that the generator $f$ does not depends on $Z$ -component in this example. In general, this experiment clarifies that we should set $q=\min\{K_{y}+1,K_{z}\}<4$ and keep $q=3$ for $\min\{K_{y}+1,K_{z}\}\geq 4,$ the value of $q$ is thus the theoretical convergence order, see Theorem 4.3.

Example 2

For the second example we consider the nonlinear BSDE (taken from [Zhao et al., 2010])

\begin{cases}&-dY_{t}=\frac{1}{2}[\exp(t^{2})-4tY_{t}-3\exp(t^{2}-Y_{t}\exp(-t^{2}))+Z_{t}^{2}\exp(-t^{2})]\,dt-Z_{t}\,dW_{t},\\ &Y_{T}=\ln(\sin W_{T}+3)\exp(T^{2}),\end{cases}

with the analytic solution

\begin{cases}&Y_{t}=\ln\left(\sin{W_{t}}+3\right)\exp(t^{2}),\\ &Z_{t}=\exp(t^{2})\frac{\cos{W_{t}}}{\sin{W_{t}}+3}.\end{cases}

The exact solution of $(Y_{0},Z_{0})$ is then $\left(\ln(3),\frac{1}{3}\right).$ In this example, the generator $f$ is nonlinear and depends on $t,Y_{t}$ and $Z_{t}.$ Thus, from Theorem 4.3 we see that the theoretical convergence order of our scheme for solving both $Y$ and $Z$ is $\min\{K_{y}+1,K_{z}\}$ but limited by $3.$ As clarified before, the used values of $q$ in both Table 7, 8 are the values of corresponding theoretical convergence order.

$\|Y_{0}-y_{0}^{0}\|$
	$N_{T}=8$	$N_{T}=16$	$N_{T}=32$	$N_{T}=64$	$N_{T}=128$	CR
$K_{y}=1,K_{z}=1,q=1$	2.72e-02	9.69e-03	3.87e-03	1.70e-03	7.87e-04	1.27
$K_{y}=1,K_{z}=2,q=2$	1.40e-02	3.41e-03	8.43e-04	2.10e-04	5.22e-05	2.02
$K_{y}=2,K_{z}=1,q=1$	1.17e-02	5.79e-03	2.89e-03	1.45e-03	7.24e-04	1.00
$K_{y}=2,K_{z}=2,q=2$	1.38e-03	4.60e-04	1.27e-04	3.33e-05	8.47e-06	1.85
$K_{y}=2,K_{z}=3,q=3$	6.39e-04	8.51e-05	1.13e-05	1.48e-06	1.89e-07	2.93
$K_{y}=3,K_{z}=1,q=1$	1.05e-02	5.76e-03	2.87e-03	1.44e-03	7.22e-04	0.97
$K_{y}=3,K_{z}=2,q=2$	1.44e-03	4.55e-04	1.27e-04	3.32e-05	8.48e-06	1.86
$K_{y}=3,K_{z}=3,q=3$	5.34e-04	9.44e-05	1.19e-05	1.53e-06	1.92e-07	2.88
$K_{y}=3,K_{z}=4,q=3$	2.33e-04	5.17e-05	6.55e-06	8.89e-07	1.13e-07	2.79
$K_{y}=4,K_{z}=1,q=1$	1.19e-02	5.82e-03	2.89e-03	1.45e-03	7.23e-04	1.01
$K_{y}=4,K_{z}=2,q=2$	1.38e-03	4.63e-04	1.28e-04	3.33e-05	8.48e-06	1.85
$K_{y}=4,K_{z}=3,q=3$	6.60e-04	8.63e-05	1.14e-05	1.48e-06	1.89e-07	2.94
$K_{y}=4,K_{z}=4,q=3$	3.49e-04	4.29e-05	6.04e-06	8.31e-07	1.10e-07	2.90
$K_{y}=4,K_{z}=5,q=3$	3.33e-04	4.14e-05	6.18e-06	8.90e-07	1.21e-07	2.84
$K_{y}=5,K_{z}=5,q=3$	1.13e-04	3.59e-05	5.81e-06	8.67e-07	1.20e-07	2.51
$K_{y}=6,K_{z}=6,q=3$	8.55e-05	2.13e-05	4.75e-06	7.70e-07	1.11e-07	2.40

Table 7: Errors and convergence rates for Example 2,

T=1

$\|Z_{0}-z_{0}^{0}\|$
	$N_{T}=8$	$N_{T}=16$	$N_{T}=32$	$N_{T}=64$	$N_{T}=128$	CR
$K_{y}=1,K_{z}=1,q=1$	5.80e-02	2.86e-02	1.42e-02	7.05e-03	3.52e-03	1.01
$K_{y}=1,K_{z}=2,q=2$	9.45e-03	2.53e-03	6.54e-04	1.66e-04	4.20e-05	1.96
$K_{y}=2,K_{z}=1,q=1$	5.99e-02	2.91e-02	1.43e-02	7.09e-03	3.53e-03	1.02
$K_{y}=2,K_{z}=2,q=2$	7.45e-03	2.02e-03	5.28e-04	1.35e-04	3.41e-05	1.94
$K_{y}=2,K_{z}=3,q=3$	2.25e-03	3.52e-04	4.91e-05	6.49e-06	8.35e-07	2.86
$K_{y}=3,K_{z}=1,q=1$	5.99e-02	2.91e-02	1.43e-02	7.09e-03	3.53e-03	1.02
$K_{y}=3,K_{z}=2,q=2$	7.46e-03	2.02e-03	5.28e-04	1.35e-04	3.41e-05	1.95
$K_{y}=3,K_{z}=3,q=3$	2.23e-03	3.50e-04	4.90e-05	6.48e-06	8.34e-07	2.85
$K_{y}=3,K_{z}=4,q=3$	6.84e-04	1.53e-04	2.53e-05	3.63e-06	4.86e-07	2.63
$K_{y}=4,K_{z}=1,q=1$	5.99e-02	2.91e-02	1.43e-02	7.09e-03	3.53e-03	1.02
$K_{y}=4,K_{z}=2,q=2$	7.44e-03	2.02e-03	5.28e-04	1.35e-04	3.41e-05	1.94
$K_{y}=4,K_{z}=3,q=3$	2.26e-03	3.52e-04	4.91e-05	6.49e-06	8.35e-07	2.86
$K_{y}=4,K_{z}=4,q=3$	7.10e-04	1.55e-04	2.54e-05	3.64e-06	4.86e-07	2.64
$K_{y}=4,K_{z}=5,q=3$	5.94e-04	1.53e-04	2.69e-05	3.97e-06	5.40e-07	2.55
$K_{y}=5,K_{z}=5,q=3$	5.86e-04	1.53e-04	2.69e-05	3.97e-06	5.40e-07	2.54
$K_{y}=6,K_{z}=6,q=3$	4.03e-04	1.22e-04	2.33e-05	3.63e-06	5.08e-07	2.43

Table 8: Errors and convergence rates for Example 2,

T=1

The given numerical results show that the proposed multi-step scheme works also well for a general nonlinear BSDE and is a highly effective and accurate. Similar to Example 1, from Table 7, 8 we can also observe that the results can be improved by increasing the number of steps. And the estimated convergences rate are mostly consistent with the theoretical convergence order. Moreover, we observe that all estimated convergence rates are around $2.5$ for $K\geq 5.$ The reason for this is that the approximations (when $K\geq 5$ ) are too precise with $N_{T}=8.$ For this case we need to consider a greater value for $N_{T}$ in order to obtain an estimated rate close to $3.$ The plots of $\log_{2}\left(|Y_{0}-y_{0}^{0}|\right)$ and $\log_{2}\left(|Z_{0}-z_{0}^{0}|\right)$ with respect to $\log_{2}(N_{T})$ are displayed in Figure 2.

(a)

Y

-component

(b)

Z

-component

Figure 2: Plots of

\log_{2}\left(|Y_{0}-y_{0}^{0}|\right)

and

\log_{2}\left(|Z_{0}-z_{0}^{0}|\right)

with respect to

\log_{2}(N_{T})

for

K=1,\cdots 6

for Example 2.

The Black-Scholes model

In this example we compute the price of a European call option $V(t,S_{t})$ by a BSDE where the underlying asset follows a geometric Brownian motion

dS_{t}=\mu S_{t}\,dt+\sigma S_{t}dW_{t}.

(82)

We assume that the asset pays dividends with the rate $d.$ The corresponding BSDE for the price of option can be derived by setting up a self-financing portfolio $Y_{t},$ which consists of $\pi_{t}$ assets and $Y_{t}-\pi_{t}$ bonds with risk-free return rate $r,$ which reads [Karoui et al., 1997b]

\left\{\begin{array}[]{l}dS_{t}=\mu S_{t}\,dt+\sigma S_{t}\,dW_{t},\\ -dY_{t}=\left(-rY_{t}-\frac{\mu-r+d}{\sigma}Z_{t}\right)\,dt-Z_{t}\,dW_{t},\\ \quad Y_{T}=\xi=\max(S_{T}-K,0).\end{array}\right.

(83)

$Y_{t}$ is the option value $V(t,S_{t}),$ $Z_{t}$ corresponds to the hedging strategy, $Z_{t}=\sigma S_{t}\pi_{t}.$ We see that $S_{t}$ in (83) is a forward process, this type of BSDEs is called (uncoupled) forward backward stochastic differential equation (FBSDE). The exact solution of (83) is given by the Black-Scholes model [Black and Scholes, 1973]. For $K=S=100,r=10\%,\mu=0.2,d=0,\sigma=0.25,T=0.1$ ²² 2 We take the parameter values which are used in [Ruijter and Oosterlee, 2015] for comparison purpose., one obtains the exact solution $(Y_{0},Z_{0})=\left(3.65997,14.14823\right).$ In our experiment, for each time step we generate the grid point for $S$ by using the analytic solution of the geometric Brownian motion

S_{i+1}=S_{i}\exp\left(\left(\mu-\frac{\sigma^{2}}{2}\right)h+\sigma\Delta x\right).

(84)

Generally, one can use, e.g., the Euler or the Milstein method to simulate the forward process when there is no analytic solution available.

Note that the error analysis for the proposed methods relies on the smoothness assumptions of the initial data. However, in European option pricing, the payoff function exhibits discontinuities at the strike price, this leads to a maximal error in the region of at-the-money. For this problem, the smooth technqiue proposed by Kreiss et al. in [Kreiss et al., 1970] has been widely used. To further reduce the error caused by the missing smoothness we can e.g., start the multi-step algorithm without the (smoothed) initial data. More precisely, we firstly smooth the initial data at $T.$ As mentioned before, for a $K$ -step scheme we need to start with $K=1$ and choose a extremely small time step $\Delta t$ to compute $(y_{\Gamma}^{N_{T}-j},z_{\Gamma}^{N_{T}-j})$ for $j=1,\cdots,K-1$ using the smoothed initial data. Then, for computing $(y_{\Gamma}^{N_{T}-K},z_{\Gamma}^{N_{T}-K})$ we use $y_{\Gamma}^{N_{T}-j}$ and $z_{\Gamma}^{N_{T}-j}$ only for $j=1,2,\cdots,K-1$ (without $j=0,$ namely without initial data), this computation is done by a $(K-1)$ -step scheme. Finally, we can run the $K$ -step scheme to compute $(y_{\Gamma}^{N_{T}-K-1},z_{\Gamma}^{N_{T}-K-1})$ based on $(y_{\Gamma}^{N_{T}-j},z_{\Gamma}^{N_{T}-j}),j=1,2,\cdots,K,$ and so on backwards until the initial time. We report our numerical results in Table 9 and 10.

$\|Y_{0}-y_{0}^{0}\|$
	$N=8$	$N=16$	$N=32$	$N=64$	$N=128$	CR
$K_{y}=1,K_{z}=1,q=1$	6.35e-04	2.88e-04	1.33e-04	6.78e-05	3.36e-05	1.06
$K_{y}=1,K_{z}=2,q=2$	8.63e-06	1.02e-06	3.83e-07	1.22e-07	2.46e-08	2.00
$K_{y}=2,K_{z}=1,q=1$	3.73e-04	1.70e-04	7.61e-05	3.92e-05	1.95e-05	1.06
$K_{y}=2,K_{z}=2,q=2$	4.83e-06	1.31e-06	3.13e-07	4.85e-08	2.13e-08	2.04
$K_{y}=2,K_{z}=3,q=3$	4.52e-09	3.83e-09	5.38e-10	7.70e-11	1.16e-11	2.29
$K_{y}=3,K_{z}=1,q=1$	3.11e-04	1.60e-04	7.89e-05	4.22e-05	2.15e-05	0.96
$K_{y}=3,K_{z}=2,q=2$	4.08e-06	8.78e-07	2.34e-07	8.79e-08	1.27e-08	2.00
$K_{y}=3,K_{z}=3,q=3$	2.43e-08	3.37e-09	4.23e-10	8.75e-11	7.13e-12	2.87
$K_{y}=3,K_{z}=4,q=3$	2.38e-08	3.33e-09	4.18e-10	8.69e-11	7.11e-12	2.87
$K_{y}=4,K_{z}=1,q=1$	2.30e-04	1.25e-04	5.25e-05	2.70e-05	1.32e-05	1.05
$K_{y}=4,K_{z}=2,q=2$	2.70e-06	6.17e-07	2.36e-07	5.80e-08	1.50e-08	1.84
$K_{y}=4,K_{z}=3,q=3$	1.04e-08	1.25e-09	3.00e-10	4.85e-11	4.80e-12	2.69
$K_{y}=4,K_{z}=4,q=3$	1.01e-08	1.22e-09	2.95e-10	4.79e-11	4.78e-12	2.68
$K_{y}=4,K_{z}=5,q=3$	1.01e-08	1.19e-09	2.92e-10	4.76e-11	4.77e-12	2.67
$K_{y}=5,K_{z}=5,q=3$	9.36e-09	1.68e-09	2.76e-10	2.97e-11	4.60e-12	2.78
$K_{y}=6,K_{z}=6,q=3$	2.85e-08	1.38e-09	3.14e-10	3.13e-11	2.12e-12	3.29

Table 9: Errors and convergence rates for the Black-Scholes model

$\|Z_{0}-z_{0}^{0}\|$
	$N=8$	$N=16$	$N=32$	$N=64$	$N=128$	CR
$K_{y}=1,K_{z}=1,q=1$	3.03e-03	1.45e-03	7.23e-04	3.70e-04	1.85e-04	1.00
$K_{y}=1,K_{z}=2,q=2$	9.36e-05	2.46e-05	6.67e-06	1.73e-06	4.36e-07	1.93
$K_{y}=2,K_{z}=1,q=1$	3.03e-03	1.46e-03	7.24e-04	3.71e-04	1.85e-04	1.00
$K_{y}=2,K_{z}=2,q=2$	9.36e-05	2.48e-05	6.66e-06	1.73e-06	4.35e-07	1.93
$K_{y}=2,K_{z}=3,q=3$	4.43e-08	5.05e-09	6.08e-10	7.92e-11	5.34e-12	3.20
$K_{y}=3,K_{z}=1,q=1$	3.04e-03	1.46e-03	7.24e-04	3.71e-04	1.85e-04	1.00
$K_{y}=3,K_{z}=2,q=2$	9.36e-05	2.48e-05	6.66e-06	1.73e-06	4.36e-07	1.93
$K_{y}=3,K_{z}=3,q=3$	4.47e-08	5.45e-09	6.17e-10	7.98e-11	5.30e-12	3.22
$K_{y}=3,K_{z}=4,q=3$	4.91e-08	9.42e-10	1.15e-10	1.08e-11	9.74e-12	3.10
$K_{y}=4,K_{z}=1,q=1$	3.04e-03	1.46e-03	7.24e-04	3.71e-04	1.85e-04	1.00
$K_{y}=4,K_{z}=2,q=2$	9.36e-05	2.48e-05	6.66e-06	1.73e-06	4.36e-07	1.93
$K_{y}=4,K_{z}=3,q=3$	4.45e-08	5.42e-09	6.15e-10	7.93e-11	5.27e-12	3.22
$K_{y}=4,K_{z}=4,q=3$	4.89e-08	1.07e-09	1.02e-10	1.12e-11	9.75e-12	3.12
$K_{y}=4,K_{z}=5,q=3$	2.77e-08	1.09e-09	2.18e-11	1.65e-11	6.88e-12	3.05
$K_{y}=5,K_{z}=5,q=3$	2.76e-08	1.49e-09	3.90e-11	1.61e-11	6.84e-12	3.05
$K_{y}=6,K_{z}=6,q=3$	2.89e-08	2.27e-09	2.32e-11	1.12e-11	7.50e-12	3.15

Table 10: Errors and convergence rates for the Black-Scholes model

From those tables, we clearly see that we have obtained surprisingly good accuracy. The estimated convergence rates are again consistent with the theoretical convergence order. Similar to the last two example, the approximation errors reduce mostly with the increasing number of steps $K.$ We draw the plots of $\log_{2}\left(|Y_{0}-y_{0}^{0}|\right)$ and $\log_{2}\left(|Z_{0}-z_{0}^{0}|\right)$ with respect to $\log_{2}(N_{T})$ in Figure 3.

(a)

Y

-component

(b)

Z

-component

Figure 3: Plots of

\log_{2}\left(|Y_{0}-y_{0}^{0}|\right)

and

\log_{2}\left(|Z_{0}-z_{0}^{0}|\right)

with respect to

\log_{2}(N_{T})

for

K=1,\cdots 6

for the example of the Black-Scholes model.

Two-dimensional example

For a two-dimensional example we consider the BSDE

\left\{\begin{array}[]{l}-dY_{t}=\left(Y_{t}-\frac{Z^{1}_{t}}{2}-\frac{Z^{2}_{t}}{2}\right)\,dt-Z^{1}_{t}\,dW^{1}_{t}-Z^{2}_{t}\,dW^{2}_{t},\\ \quad Y_{T}=\sin(W^{1}_{T}+W^{2}_{T}+T),\end{array}\right.

with the analytic solution

\left\{\begin{array}[]{l}Y_{t}=\sin(W^{1}_{t}+W^{2}_{t}+t),\\ Z_{t}=(\cos(W^{1}_{t}+W^{2}_{t}+t),\cos(W^{1}_{t}+W^{2}_{t}+t)),\end{array}\right.

The exact solution of $(Y_{0},Z^{1}_{0},Z^{2}_{0})$ is then $\left(0,1,1\right).$ The numerical approximations are reported in Table 11 and 12, which show that our multi-step scheme is still quite highly accurate for solving a two-dimensional BSDE.

$\|Y_{0}-y_{0}^{0}\|$
	$N=8$	$N=16$	$N=32$	$N=64$	$N=128$	CR
$K_{y}=1,K_{z}=1,q=1$	1.32e-02	6.46e-03	3.18e-03	1.57e-03	7.81e-04	1.02
$K_{y}=1,K_{z}=2,q=2$	4.72e-03	1.31e-03	3.45e-04	8.86e-05	2.24e-05	1.93
$K_{y}=2,K_{z}=1,q=1$	1.22e-02	6.31e-03	3.17e-03	1.58e-03	7.88e-04	0.99
$K_{y}=2,K_{z}=2,q=2$	1.83e-03	5.51e-04	1.48e-04	3.84e-05	9.82e-06	1.89
$K_{y}=2,K_{z}=3,q=3$	3.97e-04	6.77e-05	9.74e-06	1.30e-06	1.65e-07	2.82
$K_{y}=3,K_{z}=1,q=1$	8.59e-03	5.37e-03	2.94e-03	1.52e-03	7.76e-04	0.87
$K_{y}=3,K_{z}=2,q=2$	1.48e-03	5.01e-04	1.42e-04	3.76e-05	9.69e-06	1.82
$K_{y}=3,K_{z}=3,q=3$	3.94e-04	6.75e-05	9.72e-06	1.30e-06	1.64e-07	2.82
$K_{y}=3,K_{z}=4,q=3$	1.88e-04	3.76e-05	5.68e-06	7.73e-07	9.78e-08	2.74
$K_{y}=4,K_{z}=1,q=1$	5.44e-03	4.47e-03	2.70e-03	1.46e-03	7.61e-04	0.73
$K_{y}=4,K_{z}=2,q=2$	1.14e-03	4.54e-04	1.36e-04	3.68e-05	9.60e-06	1.74
$K_{y}=4,K_{z}=3,q=3$	2.91e-04	5.99e-05	9.21e-06	1.27e-06	1.63e-07	2.72
$K_{y}=4,K_{z}=4,q=3$	1.90e-04	3.78e-05	5.69e-06	7.73e-07	9.77e-08	2.75
$K_{y}=4,K_{z}=5,q=3$	1.42e-04	3.65e-05	5.99e-06	8.46e-07	1.09e-07	2.61
$K_{y}=5,K_{z}=5,q=3$	1.39e-04	3.65e-05	5.99e-06	8.46e-07	1.09e-07	2.61
$K_{y}=6,K_{z}=6,q=3$	8.12e-05	3.07e-05	5.49e-06	7.98e-07	1.05e-07	2.45

Table 11: Errors and convergence rates for the two-dimensional example

$\left(\|Z^{1}_{0}-z_{0}^{0,1}\|+\|Z^{2}_{0}-z_{0}^{0,2}\|\right)/2$
	$N=8$	$N=16$	$N=32$	$N=64$	$N=128$	CR
$K_{y}=1,K_{z}=1,q=1$	3.02e-02	4.77e-03	3.26e-03	1.87e-03	9.86e-04	1.12
$K_{y}=1,K_{z}=2,q=2$	8.40e-03	2.31e-03	6.05e-04	1.54e-04	3.92e-05	1.94
$K_{y}=2,K_{z}=1,q=1$	1.49e-02	3.92e-03	3.05e-03	1.82e-03	9.82e-04	0.90
$K_{y}=2,K_{z}=2,q=2$	9.07e-03	2.51e-03	6.60e-04	1.69e-04	4.27e-05	1.94
$K_{y}=2,K_{z}=3,q=3$	1.43e-03	2.08e-04	2.79e-05	3.59e-06	4.41e-07	2.92
$K_{y}=3,K_{z}=1,q=1$	6.47e-03	2.99e-03	2.78e-03	1.75e-03	9.67e-04	0.63
$K_{y}=3,K_{z}=2,q=2$	7.89e-03	2.37e-03	6.41e-04	1.66e-04	4.24e-05	1.89
$K_{y}=3,K_{z}=3,q=3$	1.43e-03	2.08e-04	2.79e-05	3.59e-06	4.40e-07	2.92
$K_{y}=3,K_{z}=4,q=3$	8.03e-04	1.23e-04	1.67e-05	2.16e-06	2.61e-07	2.90
$K_{y}=4,K_{z}=1,q=1$	6.76e-03	2.13e-03	2.50e-03	1.68e-03	9.46e-04	0.60
$K_{y}=4,K_{z}=2,q=2$	6.73e-03	2.21e-03	6.21e-04	1.64e-04	4.21e-05	1.84
$K_{y}=4,K_{z}=3,q=3$	1.26e-03	1.98e-04	2.73e-05	3.55e-06	4.40e-07	2.88
$K_{y}=4,K_{z}=4,q=3$	8.09e-04	1.23e-04	1.67e-05	2.16e-06	2.61e-07	2.90
$K_{y}=4,K_{z}=5,q=3$	7.48e-04	1.30e-04	1.83e-05	2.41e-06	2.97e-07	2.84
$K_{y}=5,K_{z}=5,q=3$	7.49e-04	1.30e-04	1.83e-05	2.41e-06	2.97e-07	2.84
$K_{y}=6,K_{z}=6,q=3$	5.98e-04	1.18e-04	1.73e-05	2.31e-06	2.87e-07	2.77

Table 12: Errors and convergence rates for the two-dimensional example

As we have concluded for the one-dimensional examples above, in this two-dimensional example we see that a smaller error value can be mostly achieved with a higher value of $K_{y},K_{z},$ i.e., more multi-steps. The convergence rates are roughly consistent with the theoretical results in Theorem 4.3. The slight deviation comes from the quadratures and the two-dimensional interpolations. The plots of $\log_{2}\left(|Y_{0}-y_{0}^{0}|\right)$ and $\log_{2}\left((|Z^{1}_{0}-z_{0}^{0,1}|+|Z_{0}^{2}-z_{0}^{0,2}|)/2\right)$ with respect to $\log_{2}(N_{T})$ are given in Figure 4.

(a)

Y

-component

(b)

Z

-component

Figure 4: Plots of

\log_{2}\left(|Y_{0}-y_{0}^{0}|\right)

and

\log_{2}\left((|Z^{1}_{0}-z_{0}^{0,1}|+|Z_{0}^{2}-z_{0}^{0,2}|)/2\right)

with respect to

\log_{2}(N_{T})

for

K=1,\cdots 6

for the two-dimensional example.

6 Conclusion

In this work, we adopt a multi-step scheme for solving BSDEs on time-space grids proposed in [Zhao et al., 2010] by using the cubic spline interpolating polynomials instead of the Lagrange interpolating polynomials in time. In [Zhao et al., 2010] the number of multi-steps are limited, because the stability condition cannot be satisfied for a high number of time levels. We find that our new proposed multi-step scheme allows for more multi-time-steps, which gives mostly a better approximation as our numerical results showed. However, the convergence order of our scheme equals the one of scheme in [Zhao et al., 2010]. The convergence order cannot be improved by using a higher value of $K.$ The reason for this is that a cubic spline is maximal fourth-order accurate. Several numerical examples are provided to demonstrate the highly effectiveness and accuracy of our multi-step scheme for solving BSDEs. In our proposed multi-step schemes, the computations among space grids at each time level are absolutly independent and should be thus parallelized. Therefore, a GPU-based parallel computing is desirable for higher dimensional problems. This will be the task of future work.

References

[Abramowitz and Stegun, 1972] Abramowitz, M. and Stegun, I. (1972). Handbook of Mathematical Functions. Dover Publications. Dover Books on Mathematics.
[Bally, 1997] Bally, V. (1997). Approximation scheme for solutions of bsde. In Karoui, N. E. and Mazliak, L., editors, Backward stochastic differential equations. Addison Wesley Longman, Harlow, UK.
[Bender and Steiner, 2012] Bender, C. and Steiner, J. (2012). Least-squares monte carlo for backward sdes. Numer. Methods Finance, 12:257–289.
[Bender and Zhang, 2008] Bender, C. and Zhang, J. (2008). Time discretization and markovian iteration for coupled fbsdes. Ann. Appl. Probab., 18:143–177.
[Black and Scholes, 1973] Black, F. and Scholes, M. (1973). The pricing of options and corporate liabilities. J. Political Economy, 81:637–654.
[Bouchard and Touzi, 2004] Bouchard, B. and Touzi, N. (2004). Discrete-time approximation and monte-carlo simulation of backward stochastic differential equations. Stoch. Proc. Appl., 111:175–206.
[Crisan and Manolarakis, 2010] Crisan, D. and Manolarakis, K. (2010). Solving backward stochastic differential equations using the cubature method: Application to nonlinear pricing. SIAM J. FINAN. MATH, 3(1):534–571.
[Douglas et al., 1996] Douglas, J., Ma, J., and Protter, P. (1996). Numerical methods for forward-backward stochastic differential equations. Ann. Appl. Probab., 6:940–968.
[Gobet et al., 2005] Gobet, E., Lemor, J. P., and Warin, X. (2005). A regression-based monte carlo method to solve backward stochastic differential equations. Ann. Appl. Probab., 15:2172–2202.
[Karoui et al., 1997a] Karoui, N. E., Kapoudjan, C., Pardoux, E., Peng, S., and Quenez, M. C. (1997a). Reflected solutions of backward stochastic differential equations and related obstacle problems for pdes. Ann. Probab., 25:702–737.
[Karoui et al., 1997b] Karoui, N. E., Peng, S., and Quenez, M. C. (1997b). Backward stochastic differential equations in finance. Math. Finance, 7(1):1–71.
[Kreiss et al., 1970] Kreiss, H. O., Thomée, V., and Widlund, O. (1970). Smoothing of initial data and rates of convergence for parabolic difference equations. Commun. Pure Appl. Math, 23(2):241–259.
[Lemor et al., 2006] Lemor, J., Gobet, E., and Warin, X. (2006). Rate of convergence of an empirical regression method for solving generalized backward stochastic differential equations. Bernoulli, 12:889–916.
[Lepeltier and Martin, 1997] Lepeltier, J. P. and Martin, J. S. (1997). Backward stochastic differential equations with continuous generator. Statist. Probab. Lett., 32(425–430).
[Ma et al., 2002] Ma, J., Protter, P., Martín, J. S., and Torres, S. (2002). Numerical method for backward stochastic differential equations. Ann. Appl. Probab., 12:302–316.
[Ma et al., 1994] Ma, J., Protter, P., and Yong, J. (1994). Solving forward-backward stochastic differential equations explicity-a four step scheme. Probab. Theory Related Fields, 98(3):339–359.
[Ma et al., 2009] Ma, J., Shen, J., and Zhao, Y. (2009). On numerical approximations of forward-backward stochastic differential equations. SIAM J. Numer. Anal., 46:2636–2661.
[Ma and Zhang, 2005] Ma, J. and Zhang, J. (2005). Representations and regularities for solutions to bsdes with reflections. Stoch. Proc. Appl., 115:539–569.
[Milsetin and Tretyakov, 2006] Milsetin, G. N. and Tretyakov, M. V. (2006). Numerical algorithms for forward-backward stochastic differential equations. SIAM J. SCI. COMPUT., 28:561–582.
[Pardoux and Peng, 1990] Pardoux, E. and Peng, S. (1990). Adapted solution of a backward stochastic differential equations. System and Control Letters, 14:55–61.
[Pardoux and Peng, 1992] Pardoux, E. and Peng, S. (1992). Backward stochastic differential equation and quasilinear parabolic partial differential equations. Lectures Notes in CSI., 176:200–217.
[Peng, 1991] Peng, S. (1991). Probabilistic interpretation for systems of quasilinear parabolic partial differential equations. Stochastics and Stochastic Reports, 37(1–2):61–74.
[Ruijter and Oosterlee, 2015] Ruijter, M. J. and Oosterlee, C. W. (2015). A fourier cosine method for an efficient computation of solutions to bsdes. SIAM J. SCI. COMPUT., 37(2):A859–A889.
[Teng, 2018] Teng, L. (2018). A review of tree-based approaches to solve forward-backward stochastic differential equations. Preprint 18/03, University of Wuppertal.
[Zhang, 2004] Zhang, J. (2004). A numerical scheme for bsdes. Ann. Appl. Probab., 14:459–488.
[Zhao et al., 2006] Zhao, W., Chen, L., and Peng, S. (2006). A new kind of accurate numerical method for backward stochastic differential equations. SIAM J. SCI. COMPUT., 28(4):1563–1581.
[Zhao et al., 2009] Zhao, W., Wang, J., and Peng, S. (2009). Error estimates of the theta-scheme for backward stochastic differential equations. Discrete Contin. Dyn. Syst. Ser. B, 12:905–924.
[Zhao et al., 2010] Zhao, W., Zhang, G., and Ju, L. (2010). A stable multistep scheme for solving backward stochastic differential equations. SIAM J. NUMER. ANAL., 48:1369–1394.