On multilevel Picard numerical approximations for high-dimensional nonlinear parabolic partial differential equations and high-dimensional nonlinear backward stochastic differential equations

Introduction and main results

Parabolic partial differential equations (PDEs) and backward stochastic differential equations (BSDEs) have a wide range of applications. To give specific examples we focus now on a number of applications in finance. There are several fundamental assumptions incorporated in the Black-Scholes model that are not met in the real-life trading of financial derivatives. A number of derivative pricing models have been developed in about the last four decades to relax these assumptions; see, e.g., for models taking into account the fact that the “risk-free” bank account has higher interest rates for borrowing than for lending, particularly, due to the default risk of the trader, see, e.g., for models incorporating the default risk of the issuer of the financial derivative, see, e.g., for models for the pricing of financial derivatives on underlyings which are not tradeable such as financial derivatives on the temperature or mortality-dependent financial derivatives, see, e.g., for models incorporating that the hedging strategy influences the price processes through demand and supply (so-called large investor effects), see, e.g., for models taking the transaction costs in the hedging portfolio into account, and see, e.g., for models incorporating uncertainties in the model parameters for the underlying. In each of the above references the value function $u$, describing the price of the financial derivative, solves a nonlinear parabolic PDE. Moreover, the PDEs for the value functions emerging from the above models are often high-dimensional as the financial derivative depends in several cases on a whole basket of underlyings and as a portfolio containing several financial derivatives must often be treated as a whole in the case where the above nonlinear effects are taken into account (cf., e.g., ). These high-dimensional nonlinear PDEs can typically not be solved explicitly and, in particular, there is a strong demand from the financial engineering industry to approximately compute the solutions of such high-dimensional nonlinear parabolic PDEs.

In the seminal papers , Pardoux & Peng developed the theory of nonlinear backward stochastic differential equations and, in particular, established a considerably generalized nonlinear Feynman-Kac formula to obtain an explicit representation of the solution of a nonlinear parabolic PDE by means of the solution of an appropriate BSDE; see also Cheridito et al. for second-order BSDEs. Discretizations of BSDEs, however, require suitable discretizations of nested conditional expectations (see, e.g., ). Discretization methods for these nested conditional expectations proposed in the literature include the ’straight forward’ Monte Carlo method, the quantization tree method (see ), the regression method based on Malliavin calculus or based on kernel estimation (see ), the projection on function spaces method (see ), the cubature on Wiener space method (see ), and the Wiener chaos decomposition method (see ). None of these discretization methods has the property that the computational effort of the method grows at most polynomially both in the dimension and in the reciprocal of the prescribed accuracy (see Subsections 4.1–4.6 below for a detailed discussion). We note that solving high-dimensional semilinear parabolic PDEs at single space-time points and solving high-dimensional nonlinear BSDEs at single time points is essentially equivalent due to the generalized nonlinear Feynman-Kac formula established by Pardoux & Peng. In recent years the concept of fractional smoothness in the sense of function spaces has been used for studying variational properties of BSDEs. This concept of fractional smoothness quantifies the propagation of singularities in time and shows that certain non-uniform time grids are more suitable in the presence of singularities; see, e.g., Geiss & Geiss , Gobet & Makhlouf or Geiss, Geiss, & Gobet for details. Also these temporal discretization methods require suitable discretizations of nested conditional expectations resulting in the same problems as in the case of uniform time grids.

In the recent article we proposed a family of approximation methods which we denote as multilevel Picard approximations (see (8) for its definition and Section 2 for its derivation). Corollary 3.18 in shows under suitable regularity assumptions (including smoothness and Lipschitz continuity) on the exact solution that the computational complexity of this algorithm 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 complement the theoretical complexity analysis of with a simulation study. Our simulations in Section 3 indicate that the computational complexity grows at most linearly in the dimension and quartically in the reciprocal of the prescribed accuracy also for several $100$-dimensional nonlinear PDEs from physics and finance with non-smooth and/or non-Lipschitz nonlinearities and terminal condition functions. The simulation results for these $100$-dimensional example PDEs are very satisfactory in terms of accuracy and speed.

Multilevel Picard approximations for high-dimensional semilinear PDEs

In Subsection 2.3 below we define multilevel Picard approximations (see (8) below) in the case of semilinear PDEs (cf. (5) in Subsection 2.2 below). These approximations have been introduced in . In Subsection 2.1 we explain the abstract idea behind multilevel Picard approximations. In Subsection 2.2 we derive a fixed-point equation for semilinear PDEs which is based on the Feynman-Kac and Bismut-Elworthy-Li formulas.

Roughly speaking, a key idea in our approach to solve high-dimensional nonlinear PDEs/BSDEs is to formulate the solution of the considered PDE/BSDE as the solution of a suitable fixed-point equation and then to approximate the fixed-point by suitable multilevel Picard approximations (see (4) below). We now first outline this idea in an abstract form (see (2)–(4) below) and, thereafter, we demonstrate (see Subsections 2.2–2.3) how this general idea is applied to high-dimensional semilinear PDEs and high-dimensional BSDEs.

2 A fixed-point equation for semilinear PDEs

To get a better understanding of the approximation scheme introduced in Subsection 2.3, we present in this subsection a rough derivation of a fixed-point equation on which the scheme (8) is based on. For this fixed-point equation we impose for simplicity of presentation appropriate additional hypotheses that are not needed for the definition of the scheme (8) (cf. (5)–(7) in this subsection with Subsection 2.3).

Combining (7) with the Feynman-Kac formula and the Bismut-Elworthy-Li formula ensures that ${\bf u}^{\infty}=\Phi({\bf u}^{\infty})$. Note that we have incorporated a zero expectation term in (7). The purpose of this term is to slightly reduce the variance when approximating the right-hand side of (7) by Monte Carlo approximations. Now we approximate the non-discrete quantities in (7) (expectation and time integral) by discrete quantities (Monte Carlo averages and quadrature formulas) with different degrees of discretization on different levels (cf. the remarks in Subsection 2.4 below). This yields a family of approximations of $\Phi$. With these approximations of $\Phi$ we finally define multilevel Picard approximations of ${\bf u}^{\infty}$ through (4) which results in the approximations (8).

3 The approximation scheme

In this subsection we introduce multilevel Picard approximations in the case of semilinear PDEs (see (8) below). To this end we consider the following setting.

4 Remarks on the approximation scheme

5 Special case: semilinear heat equations

In this subsection we specialize the numerical scheme (8) to the case of semilinear heat equations.

The proof of Proposition 2.1 is clear and therefore omitted.

6 Special case: geometric Brownian motion

In this subsection we specialize the numerical scheme (8) to the case of the forward diffusion being a geometric Brownian motion. This case often appears in the financial engineering literature.

Numerical simulations of high-dimensional nonlinear PDEs

In this section we apply the algorithm (8) to approximate the solutions of several nonlinear PDEs; see Subsections 3.1–3.5 below. The solutions of the PDEs in Subsections 3.1–3.4 are not known explicitly. The solution of the PDE in Subsection 3.5 is known explicitly. In Subsections 3.1–3.4 the algorithm is tested for a one-dimensional and a one hundred-dimensional version of a PDE. In the one-dimensional cases in Subsections 3.1–3.4 we present the error of our algorithm relative to a high-precision approximation of the exact solution of the PDE provided by a finite difference approximation scheme (see the left-hand sides of Figures 1, 2, 4, and 5 and Tables 1, 3, 5, and 7 below). In the one hundred-dimensional cases in Subsections 3.1–3.4 we present the approximation increments of our scheme to analyze the performance of our scheme in the case of high-dimensional PDEs (see the right-hand sides of Figures 1, 2, 4, and 5 and Tables 2, 4, 6, and 8 below). In Subsection 3.5 we employ the explicit formula for the solution of the considered one hundred-dimensional PDE (see (26) below) to present the error of our scheme relative to the explicitly known exact solution (see the left-hand side of Figure 7 and Table 9). Moreover, for each of the PDEs in Subsections 3.1–3.5 we illustrate the growth of the computational effort with respect to the dimension by running the algorithm for each PDE for every dimension $d\in\{5,6,\ldots,100\}$ and recording the associated runtimes (see Figures 3 and 6 and the right-hand side of Figure 7). All simulations are performed with Matlab on a 2.8 GHz Intel i7 processor with 16 GB RAM.

To obtain smoother results we average over 10 independent simulation runs. More precisely, for the numerical results in Subsections 3.1–3.3, for every $d\in\{1,100\}$ we run Matlab code 1 twice to produce one realization of

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. Moreover, for the numerical results in Subsections 3.4 and 3.5, 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.

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

Figures 1, 2, 4, 5, and the left-hand side of Figure 7 illustrate the empirical convergence of our scheme. In Figures 1, 2, and 4 (respectively 5) the left-hand side depicts for the settings of Subsections 3.1–3.3 (respectively 3.4) in the one-dimensional case the relative approximation errors

The right-hand side of the Figures 1, 2, and 4 (respectively 5) depicts for the settings of Subsections 3.1–3.3 (respectively 3.4) in the one hundred-dimensional case with $\rho_{\max}=7$ (respectively $\rho_{\max}=5$) the relative approximation increments

against the average runtime needed to compute the realizations $({\bf U}^{i,}_{\rho,\rho}(0,x_{0}))_{i\in\{1,2,\dots,10\}}$ for $\rho\in\{1,2,\dots,\rho_{\max}-1\}$. They are obtained by executing the command plotincrementvsruntime(value,time), where the Matlab function plotincrementvsruntime is presented in Matlab code 9.

Tables 1–9 present several statistics for the simulations. More precisely, Tables 1–6 (respectively 7–9) show for the settings of Subsections 3.1–3.3 (respectively 3.4) for all $d\in\{1,100\}$, $\rho\in\{1,2,\dots,7\}$ (respectively $d\in\{1,100\}$, $\rho\in\{1,2,\dots,5\}$) the average runtime needed to compute $({\bf U}^{i,}_{\rho,\rho}(0,x_{0}))_{i\in\{1,2,\dots,10\}}$, the empirical mean $\overline{{\bf U}}^{}_{\rho,\rho}(0,x_{0})=\frac{1}{10}\sum_{i=1}^{10}{\bf U}^{i,}_{\rho,\rho}(0,x_{0})$, and the empirical standard deviation $\sqrt{\frac{1}{9}\sum_{i=1}^{10}|{\bf U}^{i,}_{\rho,\rho}(0,x_{0})-\overline{{\bf U}}^{}_{\rho,\rho}(0,x_{0})|^{2}}$. Tables 1, 3, 5, 7, and 9 show additionally the relative approximation error (16). Furthermore, Tables 2, 4, 6, and 8 present the relative approximation increments (17).

Figures 3, 6, and the right-hand side of Figure 7 show the growth of the runtime of our algorithm with respect to the dimension for each of the example PDEs. More precisely, Figure 3 and the left-hand side of Figure 6 show for the settings in Subsections 3.1–3.3 the runtime needed to compute one realization of ${\bf U}^{1}_{6,6}(0,x_{0})$ against the dimension $d\in\{5,6,\ldots,100\}$. The left-hand side of Figure 3 is obtained by running Matlab code 10 in combination with Matlab codes 4, 11, and 12. The right-hand side of Figure 3 is obtained by running Matlab code 10 in combination with Matlab codes 4, 11, and 13. The left-hand side of Figure 6 is obtained by running Matlab code 10 in combination with Matlab codes 4, 11, and 14. The right-hand sides of Figures 6 and 7 show for the the settings in Subsections 3.4–3.5 the average runtime needed to compute 20 realizations of ${\bf U}^{1}_{4,4}(0,x_{0})$ against the dimension $d\in\{5,6,\ldots,100\}$. We average over 20 runs here to obtain smoother results. The right-hand side of Figure 6 is obtained by running Matlab code 10 (with line 4 in Matlab code 10 replaced by average=20; and line 5 in Matlab code 10 replaced by rhomax=4;) in combination with Matlab codes 4, 11, and 15 (with line 10 in Matlab code 4 replaced by Mf(rho,k)=rho^k;). The right-hand side of Figure 7 is obtained by running Matlab code 10 (with line 4 in Matlab code 10 replaced by average=20; and line 5 in Matlab code 10 replaced by rhomax=4;) in combination with Matlab codes 4, 11, and 16 (with line 10 in Matlab code 4 replaced by Mf(rho,k)=rho^k;).

In this subsection we discuss an example which is a special case of the recursive pricing model with default risk due to Duffie, Schroder, & Skiadas . The five-dimensional version of this example has also been used as a test example in the literature on numerical approximations of BSDEs (see, e.g., Bender, Schweizer, & Zhuo ).

2 Pricing with counterparty credit risk

In this subsection we present a numerical simulation of a semilinear PDE that arises in the valuation of derivative contracts with counterparty credit risk. The PDE is a special case of the PDEs that are, e.g., derived in Henry-Labordère and Burgard & Kjaer .

3 Pricing with different interest rates for borrowing and lending

We consider a pricing problem of an European option in a financial market with different interest rates for borrowing and lending. The model goes back to Bergman and serves as a standard example in the literature on numerical methods for BSDEs (see, e.g., ).

4 Allen-Cahn equation

In this subsection we consider the Allen-Cahn equation with a double well potential.

Matlab code 15 presents the parameter values in the case $d=100$. In the case $d=1$ line 3 of Matlab code 15 is replaced by dim=1;. The simulation results are shown in Figure 5, the right-hand side of Figure 6, and Tables 7 and 8. The left-hand side of Figure 5 suggests an empirical convergence rate close to $\nicefrac{{1}}{{4}}$ in the case $d=1$. Moreover, the right-hand side of Figure 5 suggests an empirical convergence rate close $\nicefrac{{1}}{{3}}$ in the case $d=100$.

5 An example with an explicit solution

In this subsection we discuss an example with an explicit solution whose three-dimensional version has been considered in Chassagneux .

Matlab code 16 presents the parameter values. The simulation results are shown in Figure 7 and Table 9. The left-hand side of Figure 7 suggests an empirical convergence rate close to $\nicefrac{{1}}{{4}}$.

Discussion of approximation methods from the literature

Deterministic methods for second-order parabolic PDEs are known to have exponentially growing computational effort. Since a program with $10^{80}$, say, floating point operations will never terminate (on a non-quantum computer), deterministic methods such as finite elements methods, finite difference methods, spectral Galerkin approximation methods, or sparse grid methods are not suitable for solving high-dimensional nonlinear second-order parabolic PDEs no matter what the convergence rate of the method is. For this reason we discuss only stochastic approximation methods for nonlinear second-order parabolic PDEs. In the literature we have found the following articles which propose (possibly non-implementable) stochastic approximation methods for nonlinear second-order parabolic PDEs. All of these methods except for exploit a stochastic representation with BSDEs due to Pardoux & Peng . Moreover, all of these methods except for can be described in two steps. In the first step, time in the corresponding BSDE is discretized backwards in time via an explicit or an implicit Euler-type method which was investigated in detail, e.g., in Bouchard & Touzi and Zhang . The resulting approximations involve nested conditional expectations and, therefore, are not implementable. In the second step, these conditional expectations are approximated by ’straight-forward’ Monte Carlo simulations, by the quantization tree method (proposed in ), by a regression method based on kernel-estimation or on Malliavin calculus (proposed in ), by projections on function spaces (proposed in ), or by the cubature method on Wiener space (developed in and proposed in ). The first step does not cause problems in high dimensions in the sense that the backward (explicit or implicit) Euler-type approximations converge under suitable assumptions with rate at least $0.5$ (see Theorem 5.3 in Zhang and Theorem 3.1 in Bouchard & Touzi for the backward implicit Euler-type method) and the computational effort (assuming the conditional expectations are known exactly) grows at most linearly in the dimension for fixed accuracy. For this reason, we discuss below in detail only the different methods for discretizing conditional expectations. In addition, we discuss the Wiener chaos decomposition method proposed in , the branching diffusion method proposed in , and methods based on density estimation proposed in .

A difficulty in our discussion below is that the discussed algorithms (except for the branching diffusion method) depend on different parameters and the optimal choice of these parameters is unknown since no lower estimates for the approximation errors are known. For this reason we will choose parameters which are optimal with respect to the best known upper error bound. For these parameter choices we will show below for the discussed algorithms (except for the branching diffusion method) that the computational effort fails to grow at most polynomially both in the dimension and in the reciprocal of the best known upper error bound.

see, e.g., Theorem 4.3 and Display (4.14) in Crisan & Manolarakis . Thus the computational effort $(Md)^{(\frac{1}{|\pi|^{-1/2}})^{2}}\geq(Md)^{(\frac{c}{c|\pi|^{-1/2}+c|\pi|M^{-1/2}})^{2}}$ grows at least exponentially in the reciprocal of the right-hand side of (27). This suggests an at most logarithmic convergence rate of the ’straight-forward’ Monte Carlo method. We are not aware of a statement in the literature claiming that the ’straight-forward’ Monte Carlo method has a polynomial convergence rate.

2 The quantization tree method

3 The Malliavin calculus based regression method

The Malliavin calculus based regression method has been introduced in Section 6 in Bouchard & Touzi and is based on the implicit backward Euler-type method. The algorithm involves iterated Skorohod integrals which by Display (3.2) in Crisan, Manolarakis, & Touzi can be numerically computed with $2^{d}$ many independent standard normally distributed random variables. In that case the computational effort grows exponentially fast in the dimension. We are not aware of an approximation method of the involved iterated Skorohod integrals whose computational effort does not grow exponentially fast in the dimension. Example 4.1 in Bouchard & Touzi also mentions a method for approximating all involved conditional expectations using kernel estimation. For this method we have not found an upper error estimate in the literature so that we do not known how to choose the bandwidth matrix of the kernel estimation given the number of time grid points.

4 The projection on function spaces method

5 The cubature on Wiener space method

In this form of the algorithm, the computational effort for calculating $Y^{\pi_{n}}$, which is at least the number $(N_{m,d})^{n}$ of paths to be used, grows exponentially in the reciprocal of the right-hand side of (29). To avoid this exponential growth of the computational effort in the number of cubature paths, Crisan & Manolarakis specify two methods (a tree based branching algorithm of Crisan & Lyons and the recombination method of Litterer & Lyons ) which reduce the number of nodes and which result in approximations which converge with polynomial rate; cf. Theorem 5.4 in. The constant in the upper error estimate in Theorem 5.4 in may depend on the dimension (cf. also (5.16) and the proof of Lemma 3.1 in ). Simulations in the literature on the cubature method were performed in dimension 1 (see Figures 1–4 in and Figure 1 in ) or dimension 5 (see Figures 5–6 in ). To the best of our knowledge, there exist no statement in the literature on the cubature method which asserts that the computational effort of the cubature method together with a suitable complexity reduction method grows at most polynomially both in the dimension of the PDE and in the reciprocal of the prescribed accuracy.

6 The Wiener chaos decomposition method

7 The branching diffusion method

The branching diffusion method has been proposed in Henry-Labordère ; see also the extensions to the non-Markovian case in and to nonlinearities depending on derivatives in . This method approximates the nonlinearity $f$ by polynomials and then exploits that the solution of a semilinear PDE with polynomial nonlinearity (KPP-type equations) can be represented as an expectation of a functional of a branching diffusion process due to Skorohod . This expectation can then be numerically approximated with the standard Monte Carlo method and pathwise approximations of the branching diffusion process. The branching diffusion method does not suffer from the ’curse of dimensionality by construction’ and works in all dimensions. It’s convergence rate is $0.5$ if the forward diffusion can be simulated exactly and, in general, its rate is $0.5-$ using a pathwise approximation of the forward diffusion and the multilevel Monte Carlo method proposed in Giles .

Table 10 shows that the branching diffusion approximations of $u^{\infty}(0,0)$ become poor as $g(0)$ increases from $0.1$ to $0.7$. Thus the branching diffusion method fails to produce good approximations for $u^{\infty}(0,0)$ in our example as soon as condition (30) is not satisfied.

8 Approximations based on density representations

This upper bound becomes only small if we choose the bandwidth ${\varepsilon}$ small and if $n$ and $N$ grow exponentially in the dimension. The upper bounds established in are less explicit in the dimension. However, following the estimates in the proofs in , it becomes apparent that the number of initial particles in branching particle system approximations defined on pages 30 and 18 in need to grow exponentially in the dimension.

This project has been partially supported through the research grants ONR N00014-13-1-0338 and DOE DE-SC0009248 and through the German Research Foundation via RTG 2131 High-dimensional Phenomena in Probability – Fluctuations and Discontinuity and via research grant HU 1889/6-1.

References