A Numerical Scheme For High-dimensional Backward Stochastic Differential Equation Based On Modified Multi-level Picard Iteration

TL;DR

This paper introduces a novel numerical scheme for high-dimensional backward stochastic differential equations using a modified multi-level Picard iteration, improving computational complexity and providing explicit error estimates.

Contribution

It presents a new numerical scheme based on a modified multi-level Picard iteration that enhances efficiency for high-dimensional BSDEs and includes explicit error bounds.

Findings

Improved complexity over traditional methods

Explicit error estimates for generator-independent cases

Effective handling of high-dimensional problems

Abstract

In this paper, we propose a new kind of numerical scheme for high-dimensional backward stochastic differential equations based on modified multi-level Picard iteration. The proposed scheme is very similar to the original multi-level Picard iteration but it differs on underlying Monte-Carlo sample generation and enables an improvement in the sense of complexity. We prove the explicit error estimates for the case where the generator does not depend on control variate.

Equations

$$y_t=\xi +\underset{t}{\overset{T}{\int }}f(s,y_s,z_s)ds-\underset{t}{\overset{T}{\int }}{z}_{s}dWs,t\in [0,T]$$

$$y_t=u(t,W_t),z_t={\nabla }_{x}u(t,W_t),\forall t\in [0,T)$$

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

$$u(t,x)$$

$$\frac{\partial u}{\partial x_k}(t,x)$$

$$y_0(t,x)=0,z_0(t,x)=0$$

$$y_n(t,x)$$

$$+$$

$$z_n(t,x)$$

$$+$$

$$y_0(t,x)=0,z_0(t,x)=0$$

$$y_1(t,x)=\frac{1}{M}\sum _{i=1}^M\phi \left(x+W_{T-t}^{(1,i)}\right)+\sum _{i=1}^Q{w}_{1,j}f\left(t_{1,j},0,0\right)$$

$$z_1(t,x)=\frac{1}{M}\sum _{i=1}^M\left(\phi \left(x+W_{T-t}^{(1,i)}\right)-\phi (x)\right)\frac{W_{T-t}^{(1,i)}}{T-t}+$$

$$+\sum _{i=1}^Q\frac{w_{1,j}}{M}\sum _{i=0}^{M}f\left(t_{1,j},0,0\right)\frac{W_{t_{1,j}-t}^{(1,i)}}{t_{1,j}-t}$$

$$y_n(t,x)=\frac{1}{M}\sum _{i=1}^M{y}_{n-1}^i(t,x)$$

$$+\sum _{j=1}^Q\frac{w_{n,j}}{M}\sum _{i=1}^M\left[\left(f\left(y_{n-1},z_{n-1}\right)-f(y_{n-2},z_{n-2})\right)\left(t_{n,j},x+W_{t_{n,j}-t}^{(n,i)}\right)\right]$$

$$z_n(t,x)=\frac{1}{M}\sum _{i=1}^M{z}_{n-1}^i(t,x)\frac{W_{T-t}^{(n,i)}}{T-t}$$

$$+\sum _{j=1}^Q\frac{w_{n,j}}{M}\sum _{i=1}^M\left[\left(f\left(y_{n-1},z_{n-1}\right)-f(y_{n-2},z_{n-2})\right)\left(t_{n,j},x+W_{t_{n,j}-t}^{(n,i)}\right)\frac{W_{t_{n,j}-t}^{(n,i)}}{t_{n,j}-t}\right]$$

$$y_n(t,x)$$

$$L^q(x)=\frac{1}{2^{q}q!}\frac{d^q}{d{x}^q}[(x^2-1{)}^q]$$

$$\underset{a}{\overset{b}{\int }}g(t)dt-\sum _{j=1}^q{w}_{j}g(t_j)=\frac{{\left[q!\right]}^4(b-a{)}^{2q+1}}{(2q+1){\left[(2q)!\right]}^3}{g}^{(2q)}(\xi ),\xi \in [a,b]$$

$$t_j=\frac{c_j^q(b-a)+(a+b)}{2},w_j=\underset{a}{\overset{b}{\int }}\left[\prod _{i\in \{1,\dots ,q\},t_i\ne t_j}\frac{2x-(b-a)c_i^q-(a+b)}{2t_j-(b-a)c_i^q-(a+b)}\right]dx$$

$$\underset{[a,b],q}{\overset{\sim }{\int }}g(t)dt=\sum _{j=1}^q{w}_{j}g(t_j)$$

$$\forall n\in \mathbb{N},\forall t\in [0,T],g^{(2n)}(t)\ge 0$$

$$\forall q\in \mathbb{N},\underset{[a,b],q}{\overset{\sim }{\int }}g(t)dt=\sum _{j=1}^q{w}_{j}g(t_j)\le \underset{a}{\overset{b}{\int }}g(t)dt,$$

$$y_1(t,x)=\frac{1}{M}\sum _{i=1}^M\phi \left(x+W_{T-t}^{(1,i)}\right)+\sum _{j=1}^Q{w}_{1,j}f(t_{1,j},0)$$

$$\left({\Omega }_1,{\mathcal{F}}_1,P_1\right)=\left(\Omega ,\mathcal{F},P\right).$$

$$y_2(t,x)=\frac{1}{M}\sum _{i=1}^M{y}_1^i(t,x)+\sum _{j=1}^Q\frac{w_{2,j}}{M}\sum _{i=1}^M\left[f\left(t_{2,j},y_1\left(W_{t_{2,j}-t}^{(2,i)}\right)\right)-f(t_{2,j},0)\right]$$

$$\forall n\in \mathbb{N},\sqrt{2\pi n}{\left(\frac{n}{e}\right)}^n\le n!\le \sqrt{2\pi n}{\left(\frac{n}{e}\right)}^n{e}^{\frac{1}{12}}$$

$$\forall n\in\mathbb{N},k<n,\left(\begin{array}[]{c}n\\ k\end{array}\right)<2^{n}$$

$$\exists C_f>0,\forall y_1,y_2\in \mathbb{R},\forall s\in [0,T],|f(s,y_1)-f(s,y_2)|\le C_f|y_1-y_2|$$

$$\exists C_0>0,\forall s\in [0,T],|f(s,0)|\le C_0$$

$$\exists C_{\phi }>0,\forall x\in {\mathbb{R}}^d,|\phi (x)|\le C_{\phi }$$

$$\exists C_y>0,\forall t\in [0,T],\forall x\in {\mathbb{R}}^d,|y(t,x)|\le C_y$$

$$F_t(s):=E\left[f\left(s,y\left(s,x+W_{s-t}\right)\right)\right],G_t(s):=E\left[f\left(s,y\left(s,x+W_{s-t}\right)\right)\frac{W_{s-t}}{s-t}\right]$$

Full text

We prove the explicit error estimates for the case where the generator does not depend on control variate.

Keywords backward stochastic differential equations; numerical scheme; curse-of-dimensionality; error estimates

1 Introduction

Let $\left(\Omega,\mathcal{F},P\right)$ be a probability space, $T>0$ a finite time and $\left\{\mathcal{F}_{t}\right\}_{0\leq t\leq T}$ a filtration satisfying the usual conditions. Let $\left(\Omega,\mathcal{F},P,\left\{\mathcal{F}_{t}\right\}_{0\leq t\leq T}\right)$ be a complete, filtered probability space on which a standard $d$-dimensional Brownian motion $W_{t}=\left(W_{t}^{1},W_{t}^{2},\ldots,W_{t}^{d}\right)^{T}$ is defined and $\mathcal{F}_{0}$ contains all the $P-$null sets of $\mathcal{F}$. Let $L^{2}=L_{\mathcal{F}}^{2}(0,T)$ be the set of all $\{\mathcal{F}_{t}\}-$adapted mean$-$square$-$integrable processes.

We consider the backward stochastic differential equation (BSDE)

[TABLE]

where the generator $f=f(t,y,z)$ is a vector function valued in $\mathbb{R}^{m}$ and is $\mathcal{F}_{t}-$adapted for each $(y,z)$ and the terminal variable $\xi\in L^{2}$ is $\mathcal{F}_{T}-$ measurable.

A process $(y_{t},z_{t}):[0,T]\times\Omega\rightarrow\mathbb{R}^{m}\times\mathbb{R}^{m\times d}$ is called an $L^{2}-$solution of the BSDE (1.1) if it is $\{\mathcal{F}_{t}\}-$adapted, square integrable, and satisfies the equation. In 1990, Pardoux and Peng first proved in [6] the existence and uniqueness of the solution of general nonlinear BSDEs and afterwards there has been very active research in this field with many applications.([5, 7, 9, 10] )

In this paper we assume that the terminal condition is a function of $W_{T}$, i.e. $\xi=\varphi(W_{T})$ and the BSDE (1.1) has a unique solution $(y_{t},z_{t})$.

It was shown in [7] that the solution $(y_{t},z_{t})$ of (1.1) can be represented as

[TABLE]

where $u(t,x)$ is the solution of the following parabolic partial differential equation

[TABLE]

with the terminal condition $u(T,x)=\varphi(x)$, and $\nabla_{x}u$ is the gradient of $u$ with respect to the spacial variable $x$. The smoothness of $u$ depends on $f$ and $\varphi$.

Furthermore, the celebrated Feynman-Kac formula holds as follows.

[TABLE]

Moreover, regarding to the gradient of solution, Bismut-Elworthy-Li formula ([4]) holds as follows.

[TABLE]

Although BSDEs have very important applications in many fields such as mathematical finance and stochastic control, it is well known that it is difficult to obtain analytic solutions except some special cases and there have been many works on numerical methods to get approximate solution.

In [5] a numerical method by using binomial approach was proposed and in [9, 10] a $\theta$-scheme which is based on time-space grid is proposed. Besides, there are several kinds of numerical methods for BSDEs (see references in [2]).

All the algorithms mentioned above suffer from curse-of-dimensionality, i.e., the computational complexity grows exponentially in dimensionality. Most of the problems from real applications are high-dimensional and it has been very challenging to build a scheme of which complexity grows polynomially in dimension.

Recently Martin et al. [1, 2] proposed such a scheme for the first time based on Monte-Carlo simulation and multi-level Picard iteration. Their multi-level approach is very effective and it shows applications in many important numeric fields. They proved the error estimates of the multi-level Picard approximation rigorously in [3].

In this paper, we propose a scheme which is based on modified multi-level Picard iteration. The proposed scheme is very similar to the original multi-level scheme but makes use of different underlying Monte-Carlo samples and it enables a slight improvement in complexity. We note that the improvement just lies in constant scale (approximately $1.5\sim 2$ times). But we still believe the idea should be effective and applicable to other applications.

The rest of this paper is organized as follows. In Section 2 we introduce the scheme based on original multi-level Picard iteration and present a new scheme. In Section 3, we give error estimates of the proposed scheme theoretically under some assumptions for a limited case. In Section 4, some conclusions are given.

2 Multi-level Picard iteration and modification

First we introduce the multi-level Picard iteration method proposed in [1].

Initial step ($n=0$)

[TABLE]

Iteration step ($n\geq 1$)

[TABLE]

In the above scheme, $W_{s-t}^{(n,i)}$ denotes the $i$-th $d$-dimensional Gauss random variable that is used for Monte-Carlo approximation of the expectation at $n$-th iteration stage and $\left\{W_{s-t}^{(n,i)}\right\}$ are mutually independent.

Gauss-Legendre quadrature with $Q-$points is used for the approximation of the time integral. $(t_{n,1},\ldots,t_{n,Q})\in[t,T]$ are the quadrature points and $w_{n,j}$ denotes the weight for the point $t_{n,j}$.

We propose a modified multi-level Picard approximation scheme as follows.

Initial step ($n<2$)

[TABLE]

Iteration step ($n\geq 2$)

[TABLE]

[TABLE]

In the above scheme, $\{y_{n}^{i}(t,x)\}$ are the i.i.d. random variables that is identically distributed with $y_{n}(t,x)$.

The proposed scheme (2.4)-(2.7) is very similar to the original multi-level Picard iteration (2.1)-(2.3). Actually if we expand $y_{n}(t,x)$ one step further, we have the following.

[TABLE]

Likewise if we repeatedly expand we have a very similar one to (2.1)-(2.3). But note that the proposed scheme makes use of Monte-Carlo samples of different level to calculate $y_{n}(t,x)$ while (2.2)-(2.3) uses the ones of the same level. More importantly, we need to calculate $y_{n-1}(t,x)$ to get $y_{n}(t,x)$. So for the evaluation of $f(y_{n-2},z_{n-2})$ at $\left(t_{n,j},x+W_{t_{n,j}-t}^{(n,i)}\right)$ on the right hand of (2.6)-(2.7) we can use the one that we used for evaluation of $f(y_{n-1},z_{n-1})$ at the same time-space point. This reduces the computational complexity slightly, to half roughly.

3 Error estimates

In this section, we prove the convergence of the proposed scheme for a limited case where the driver $f(s,y,z)$ does not depend on control variate. We also note that we just present the proof for $y$ component.

First let us point out some results on Gauss-Legendre quadrature.

Let $q\in\mathbb{N}$ , let $(C_{j}^{q})^{j-1,q}\subset[-1,1]$ be the roots of $q$-th order Legendre polynomial

[TABLE]

. For the real function $g(s):[a,b]\rightarrow\mathbb{R}$, its integral $\int\limits_{a}^{b}g(t)dt$ can be approximated using $q$ points Gauss-Legendre quadrature as follows.

[TABLE]

where $(t_{j})_{j=1\cdots q}$ are the quadrature points and $(w_{j})_{j=1\cdots q}$ are their corresponding weights defined as follows.

[TABLE]

For the simplicity, in the rest of the paper we denote the $q$-points Gauss-Legendre quadrature approximation of $\int\limits_{a}^{b}g(t)dt$ by $\int\limits_{[a,b],q}^{\sim}g(t)dt$ . Moreover as long as the number of quadrature points does not change, we simply write as $\int\limits_{[a,b]}^{\sim}g(t)dt$.

[TABLE]

The following proposition is a direct result from (3.1) and it will be used in many places later.

Proposition 3.1**.**

For a sufficiently smooth real function $g:[a,b]\rightarrow\mathbb{R}$ which satisfies

[TABLE]

the following inequality holds

[TABLE]

where $\int\limits_{[a,b],q}^{\sim}g(t)dt=\sum_{j=1}^{q}w_{j}g(t_{j})$ denotes the $q$-points Gauss-Legendre quadrature approximation of $\int\limits_{a}^{b}g(t)dt$.

Now we make it clear the meaning of $y_{n}(t,x),z_{n}(t,x)$.

For $n=1$ , we have

[TABLE]

and $y_{1}(t,x)$ is the random variable in the probability space

[TABLE]

For $n=2$ , we have

[TABLE]

and the right hand includes $y_{1}(t_{2,j},x+W_{t_{2,j}-t}^{(2,i)})$ . So $y_{2}(t,x)$ is the random variable in the probability space$\left(\Omega_{2},\mathcal{F}_{2},P_{2}\right)=\left(\Omega\times\Omega,\mathcal{F}\otimes\mathcal{F},P\times P\right)$ which depends on $W^{(2,i)}$ and $W^{(1,i)}$.

Likewise for the $n$-th iteration stage, $y_{n}(t,x)$ is the random variable in the probability space $\left(\Omega_{n},\mathcal{F}_{n},P_{n}\right)=\left(\underbrace{\Omega\times\Omega\times\cdots\times\Omega}_{n},\underbrace{\mathcal{F}\otimes\mathcal{F}\otimes\cdots\otimes\mathcal{F}}_{n},\underbrace{P\times P\times\cdot\times P}_{n}\right)$ . We denote the expectation, variance and $L_{p}$-norm in $\left(\Omega_{n},\mathcal{F}_{n},P_{n}\right)$ by $E^{n}[\cdot]$, $var^{n}[\cdot]$ and $\|\|_{n,p}$ respectively.

The following propositions will be used in the error estimates

Proposition 3.2**.**

[8]**

[TABLE]

Proposition 3.3**.**

[TABLE]

Now we address some necessary assumptions for the error estimates.

Assumption 1**.**

The generator $f(s,y)$ is Lipschitz continuous in $y$ and $f(s,0)$ is globally bounded.

[TABLE]

Assumption 2**.**

The terminal condition $\varphi(x)$ is bounded.

[TABLE]

Assumption 3**.**

The true solution $y(t,x)$ is bounded

[TABLE]

Assumption 4**.**

For any $t\in[0,T],x\in\mathbb{R}^{d}$ , if we define two real functions as follows

[TABLE]

, then $F_{t}(s),G_{t}(s)$ are bounded, smooth enough and all of their derivatives are also bounded.

[TABLE]

It is not difficult to check that the Assumption (4) holds if the generator and the true solution are smooth enough, bounded and all of their derivatives are also bounded. (See [10])

Now we state the error estimates of the proposed scheme (2.4)-(2.7).

Theorem 3.4**.**

Let $(y(t,W_{t}),z(t,W_{t}))$ be the solution of the following backward stochastic differential equation.

[TABLE]

Let $(y_{n},z_{n})$ be the approximation series by the scheme (2.4)-(2.7). Under the Assumptions (1)-(4), the following estimate holds.

[TABLE]

[TABLE]

*where $C_{1},C_{2}$ are some constants that depend only on $T,C_{f},C_{y},C_{0},C_{d},C_{\varphi}$. *

Proof 3.5**.**

For $n\geq 2$ from (3.2) it holds that

[TABLE]

At the last step of the above equality the following reasoning was applied

[TABLE]

where $p_{t_{n,j}-t}(u)$ denotes the probability density function of Gaussian random variable $W_{t_{n,j}-t}^{(n)}$.

[TABLE]

Repeating the similar procedure, we deduce the following result for $n\geq 2$.

[TABLE]

*Note that (3.4) holds for $n=1$ from (2.5).

For all $n\in\mathbb{N}$ , from (3.4) using triangle inequality and Assumption (1) it holds that*

[TABLE]

where

[TABLE]

From (3.1) and Assumption(4) we have

[TABLE]

Now let us evaluate the variance in each iteration step. From the independence of $(y_{n-1}^{i})_{i=1,\ldots,M}$ and the independence of $(W_{t_{n,j}-t}^{n,i})_{i=1,\ldots,M}$, using the Liptschtz property of $f$ we deduce that

[TABLE]

Taking squared roots of the both side, it holds that

[TABLE]

On the other hand, from (2.6) it holds that

[TABLE]

Let us define

[TABLE]

Then for any $n\geq 2$ , from (3.8) and (3.9) it holds that

[TABLE]

[TABLE]

For the initial step, from the Assumption (1), Assumption (2) and (2.1)-(2.2) we deduce that

[TABLE]

[TABLE]

Now let us define $a:=1+\dfrac{1}{\sqrt{M}}$ and substituting (3.11) into (3.10) repeatedly it holds that

[TABLE]

Because it holds that $\dfrac{d^{2q}}{ds^{2q}}\left(\dfrac{(T-t)^{k}}{k!}\right)\geq 0$ for any $k,q\in\mathbb{N}$ , from (3.13) and Proposition 3.1 we deduce that

[TABLE]

For $n\geq 2$, combining (3.14) and (3.15) it holds that

[TABLE]

*Now we prove the following inequality by induction on $2\leq n\in\mathbb{N}$.

[TABLE]

For the base case where $n=2$, from (3.10)-(3.13) it holds that

[TABLE]

*Now let us assume that (LABEL:3.17) holds for $v_{1}(t),\ldots,v_{n-1}(t)$.

*Then, for any $1\leq l\leq n-1$ the following inequality holds.

[TABLE]

If we apply this to each term on the right hand of (LABEL:3.17) then the coefficient of $\dfrac{C_{\varphi}}{\sqrt{M}}(T-t)^{k}$ from the first part is as follows.

[TABLE]

From Pascal’s formula we deduce that

[TABLE]

*Likewise the coefficient of $\dfrac{C_{\varphi}+C_{0}T}{\sqrt{M}}(T-t)^{k}$ from the second part becomes $\binom{n-2}{k-1}$ and (LABEL:3.17) holds for $n$.

Now from (3.13) we have $C_{\varphi}\leq C_{\varphi}+C_{0}T$ and it holds that*

[TABLE]

Now setting $C_{1}=\max\{2a,C_{\varphi}+C_{0}T\}$ , from the Proposition 3.3 the following inequality holds.

[TABLE]

Let us define $\mu_{n}(t):=\sup_{x\in\mathbb{R}^{d}}\left[|(y-E^{n}[y_{n}])(t,x)|\right]$ then from the Assumption 3 we have $\mu_{0}(t)\leq C_{y}$ and from (3.5) we deduce that

[TABLE]

From (3.6) and the Proposition 3.1-3.2, the first sum on the right hand of (3.20) satisfies the following inequality

[TABLE]

where

[TABLE]

From the Taylors expansion of $e^{x}$ at 0, it holds that

[TABLE]

From (3.19) it holds that

[TABLE]

[TABLE]

In a similar way, one can easily check the following inequality by induction on $k\in\mathbb{N}$.

[TABLE]

Summing up (LABEL:3.22) for $k=1,\ldots,n$, it holds that

[TABLE]

Likewise the last term of (3.20) satisfies the following inequality.

[TABLE]

Substituting (LABEL:3.21),(LABEL:3.23),(LABEL:3.24) into (3.20), it holds that

[TABLE]

*(3.19) and (LABEL:3.25) proves the result of the theorem. *

Note that the main contribution of the error estimates given by the Theorem 3.4 is from $\left(\dfrac{C_{1}}{\sqrt{M}}\right)^{n}$.

Because the complexity of the scheme (2.4)-(2.7) grows at about $d(MQ)^{n}$, we can see that the complexity grows polynomially in both dimension and error by choosing $M$ and $Q$ properly. (See [1, 3] for details)

4 Conclusion

In this paper we proposed a modified multi-level Picard iteration scheme for backward stochastic differential equation and presented an explicit error estimates of the scheme. The proposed scheme is very similar to the original multi-level Picard iteration scheme but it differs slightly and enables improvement in computational complexity. Our further interests will lie on the error estimates of the general case where the generator depends on control variate. Application of multi-level Picard approximation in other fields would also be our interests.