On Multilevel Picard Numerical Approximations for High-Dimensional Nonlinear Parabolic Partial Differential Equations and High-Dimensional Nonlinear Backward Stochastic Differential Equations

Abstract

Parabolic partial differential equations (PDEs) and backward stochastic differential equations (BSDEs) are key ingredients in a number of models in physics and financial engineering. In particular, parabolic PDEs and BSDEs are fundamental tools in pricing and hedging models for financial derivatives. The PDEs and BSDEs appearing in such applications are often high-dimensional and nonlinear. Since explicit solutions of such PDEs and BSDEs are typically not available, it is a very active topic of research to solve such PDEs and BSDEs approximately. In the recent article (E et al., Multilevel Picard iterations for solving smooth semilinear parabolic heat equations, arXiv:1607.03295) we proposed a family of approximation methods based on Picard approximations and multilevel Monte Carlo methods and showed under suitable regularity assumptions on the exact solution of a semilinear heat equation that the computational complexity is bounded by $O( d \, {\varepsilon }^{-(4+\delta )})$ for any $\delta \in (0,\infty )$ where d is the dimensionality of the problem and ${\varepsilon }\in (0,\infty )$ is the prescribed accuracy. In this paper, we test the applicability of this algorithm on a variety of 100-dimensional nonlinear PDEs that arise in physics and finance by means of numerical simulations presenting approximation accuracy against runtime. The simulation results for many of these 100-dimensional example PDEs are very satisfactory in terms of both accuracy and speed. Moreover, we also provide a review of other approximation methods for nonlinear PDEs and BSDEs from the scientific literature.

Appendix

Here we provide the Matlab codes needed to approximate the solutions of the example PDEs from Sect. 3. Throughout this section assume the setting in Sect. 2.3.

Matlab code 1 below produces one realization of

$$\begin{aligned}&\{ 1, 2, \dots , 7 \}\times \{ 1, 2, \dots ,5 \} \ni (\rho ,i) \mapsto \mathbf{U}^i_{ \rho , \rho }(0, x_0 )\nonumber \\&=(\mathbf{U}^{i,[1]}_{ \rho , \rho }( 0,x_0 ),\mathbf{U}^{i,[2]}_{ \rho , \rho }( 0,x_0 ),\ldots ,\mathbf{U}^{i,[d+1]}_{ \rho , \rho }(0, x_0 ))\in {\mathbb {R}}^{d+1}. \end{aligned}$$

(34)

For the numerical results in Sects. 3.1, 3.2, for every $d\in \{1,100\}$ we run Matlab code 1 twice where, in the second run, line 2 of Matlab code 1 is replaced by rng(2017) to initiate the pseudorandom number generator with a different seed. This way we obtain in total 10 independent simulation runs. Moreover, for the numerical results in Sects. 3.3, we run Matlab code 1 once, where lines 4, 5, and 14 are replaced by average=10;, rhomax=5;, and [a,b]=approximateUZabm(n(rho),rho,zeros(dim,1),0);, respectively.

Matlab code 1 calls the Matlab functions approximateUZgbm (respectively approximateUZabm), modelparameters, and approxparameters. The Matlab functions approximateUZgbm and approximateUZabm are presented in Matlab codes 2 and 3 and implement the schemes (12) and (11), respectively. More precisely, up to rounding errors and the fact that random numbers are replaced by pseudo random numbers, it holds for all $\theta \in \Theta$, $n\in {\mathbb {N}}_0$, $\rho \in {\mathbb {N}}$, $x\in {\mathbb {R}}^d$, $s\in [0,T)$ that $\texttt {approximateUZgbm}(n,\rho ,x,s)$ returns one realization of $\mathbf{U}^\theta _{ n, \rho }(s, x)$ satisfying (12). Moreover, up to rounding errors and the fact that random numbers are replaced by pseudo random numbers, it holds for all $\theta \in \Theta$, $n\in {\mathbb {N}}_0$, $\rho \in {\mathbb {N}}$, $x\in {\mathbb {R}}^d$, $s\in [0,T)$ that $\texttt {approximateUZabm}(n,\rho ,x,s)$ returns one realization of $\mathbf{U}^\theta _{ n, \rho }(s, x)$ satisfying (11).

The Matlab function modelparameters called in line 7 of Matlab code 1 returns the parameters $T\in (0,\infty )$, $d\in {\mathbb {N}}$, $f:[0,T]\times {\mathbb {R}}^d \times {\mathbb {R}}\times {\mathbb {R}}^{d} \rightarrow {\mathbb {R}}$, $g:{\mathbb {R}}^d \rightarrow {\mathbb {R}}$, $\eta :{\mathbb {R}}^d \rightarrow {\mathbb {R}}$, ${{\bar{\mu }}}\in {\mathbb {R}}$, and ${{\bar{\sigma }}} \in {\mathbb {R}}$ for each example considered in Sects. 3.1–3.3. Matlab code 4 presents the implementation for the setting in Sect. 3.1 in the case $d=100$ and $T=2$.

The Matlab function approxparameters called in line 8 of Matlab code 1 provides for every example considered in Sects. 3.1, 3.2 (respectively Sect. 3.3) and every $\rho \in \{1,2,\ldots , 7\}$ (respectively $\rho \in \{1,2,\ldots , 5\}$) the numbers of Monte-Carlo samples $(m^g_{k,l,\rho })_{k,l \in {\mathbb {N}}_0}$ and $(m^f_{k,l,\rho })_{k,l \in {\mathbb {N}}_0}$ and the quadrature formulas $(q^{k,l,\rho }_s)_{k,l \in {\mathbb {N}}_0, s\in [0,T)}$. More precisely, we assume for every $s \in [0,T], k,l \in {\mathbb {N}}_0, \rho \in {\mathbb {N}}$ with $k\ge l$ that $q^{ k, l, \rho }_s$ is the Gauss–Legendre quadrature formula on (s, T) with ${\text {round}}(\varphi (\rho ^{(k-l)/2}))$ nodes, where $\varphi :[1,\infty ) \rightarrow [2,\infty )$ is the approximation of the inverse gamma function provided by Matlab code 6. To compute the Gauss–Legendre nodes and weights we use the Matlab function lgwt that was written by Greg von Winckel and that can be downloaded from www.mathworks.com. In addition, for every $k,l \in {\mathbb {N}}_0, \rho \in {\mathbb {N}}$ we choose in Sects. 3.1, 3.2 that $m^f_{ k, l, \rho } = {\text {round}}(\rho ^{(k-l)/2})$ and $m^g_{ k, l, \rho } = \rho ^{k-l}$ and in Sect. 3.3 that $m^f_{ k,l, \rho } = \rho ^{k-l}$ and $m^g_{ k, l, \rho } = \rho ^{k-l}$. For the numerical results in Sects. 3.1, 3.2Matlab code 5 presents the implementation of approxparameters. For the numerical results in Sect. 3.3 line 10 in Matlab code 5 is replaced by $\texttt {Mf(rho,k)=rho}^{\wedge }{} \texttt {k;}$. The reason for choosing in Sects. 3.1, 3.2 fewer Monte-Carlo samples $(m^f_{k,l,\rho })_{k,l \in {\mathbb {N}}_0, \rho \in {\mathbb {N}}}$ than in Sect. 3.3 is that in the former cases for every $s \in [0,T)$ the variance ${{\text {Var}}}(f(s,X^{0,x_0}_s, {\mathbb {E}}[g(X^{s,x}_T)(1,\frac{W_T-W_s}{T-s})]\big |_{x=X^{0,x_0}_s}))$ of the nonlinearity is of smaller magnitude than the variance ${{\text {Var}}}(g(X^{0,x_0}_T))$ of the terminal condition. Therefore, the nonlinearity requires fewer Monte-Carlo samples to obtain a Monte-Carlo error of the same magnitude as the terminal condition. Averaging the nonlinearity less saves computational effort and allows to employ a higher maximal number of Picard iterations (7 in Sects. 3.1, 3.2 compared to 5 in Sect. 3.3).

Solutions of one-dimensional PDEs can be efficiently approximated by finite difference approximation schemes. Matlab code 7 implements such an approximation scheme in the setting of Proposition 2.2 and Matlab code 8 implements such an approximation scheme in the setting of Proposition 2.1.

The command ploterrorvsruntime(v,value,time) (the matrices value and time are produced in Matlab code 1 and the value v by Matlab code 7 or 8) plots the error (17) against the runtime (cf. the left-hand side of Figs. 1–3).

The command plotincrementsvsruntime(value,time) (the matrices value and time are produced in Matlab code 1) plots the increments (18) against the runtime (cf. the right-hand side of Figs. 1–3).

Fig. 5

Left: Runtime needed to compute one realization of $\mathbf{U}^1_{6,6}(0,x_0)$ against dimension $d\in \{5,6,\ldots ,100\}$ for the pricing with counterparty credit risk example in Sect. 3.1. Middle: Runtime needed to compute one realization of $\mathbf{U}^1_{6,6}(0,x_0)$ against dimension $d\in \{5,6,\ldots ,100\}$ for the pricing with different interest rates example in Sect. 3.2. Right: Average runtime needed to compute 20 realizations of $\mathbf{U}^1_{4,4}(0,x_0)$ against dimension $d\in \{5,6,\ldots ,100\}$ for the Allen–Cahn equation in Sect. 3.3

The three graphs of Fig. 5 are produced with the help of Matlab codes 11 and 12. More precisely, up to rounding errors and the fact that random numbers are replaced by pseudo random numbers, Matlab code 11 generates for every $d\in \{5,6,\ldots ,100\}$ one realization of $\mathbf{U}^0_{ 6, 6 }(0, x_0)$ with $x_0= (100,\ldots ,100)\in {\mathbb {R}}^d$ and records the associated runtimes. Matlab code 11 calls Matlab code 12 to plot the three graphs in Fig. 5 where, for the right-hand side of Fig. 5, lines 4 and 11 in Matlab code 11 are replaced by average=20; and rhomax=4;, respectively.