A gradient method for high-dimensional BSDEs

Kossi Gnameho; Mitja Stadje; Antoon Pelsser

doi:10.1515/mcma-2024-2002

Article Open Access

A gradient method for high-dimensional BSDEs

Kossi Gnameho , Mitja Stadje and Antoon Pelsser

Published/Copyright: February 14, 2024

Published by

Become an author with De Gruyter Brill

Submit Manuscript Author Information Explore this Subject

From the journal Monte Carlo Methods and Applications Volume 30 Issue 2

Abstract

We develop a Monte Carlo method to solve backward stochastic differential equations (BSDEs) in high dimensions. The proposed algorithm is based on the regression-later approach using multivariate Hermite polynomials and their gradients. We propose numerical experiments to illustrate its performance.

Keywords: Regression; BSDE; Monte Carlo; Hermite polynomials

MSC 2020: 65C05; 65C40; 60H10

1 Introduction

This paper deals with the numerical approximation of backward stochastic differential equations (BSDEs) in high dimensions on a certain time interval [ 0 , T ] . Backward stochastic differential equations were first introduced by Bismut [6] in the linear case and later developed by Pardoux and Peng [53]. In the past decade, BSDEs have attracted a lot of attention and have been intensively studied in mathematical finance, insurance and stochastic optimal control theory. For example, in a complete financial market, the price of a standard European option can be seen as the solution of a linear BSDE. Moreover, the price of an American option can be formulated as the solution of a reflected BSDE. BSDEs have also been widely applied for portfolio optimization, indifference pricing, modelling of convex risk measures and the modelling of ambiguity with respect to the stochastic drift and the volatility. See for instance El Karoui, Peng and Quenez [24], Cheridito et al. [13], Barles and Lesigne [3], Duffie and Epstein [22], Hamadène and Jeanblanc [36], Hu, Imkeller and Müller [39], Laeven and Stadje [47] and the references therein. In general, most of these equations do not have an explicit or closed form solution. Some efforts have been have been made to provide numerical solutions. A four-step scheme has been proposed by Ma, Protter and Yong in [48] to solve forward-backward stochastic differential equations (FBSDEs). In [2], Bally has proposed a random time discretization scheme. Discrete time approximation schemes have also been proposed by Bouchard and Touzi in [7] and Chevance in [15], among others. In Chevance’s work [15], strong regularity assumptions on the coefficients on the BSDE are required for the convergence results. Approximation methods for BSDEs based on cubature are designed in [16, 10, 17, 12], among others. Gobet, Lemor and Warin [33] presented a discrete algorithm based on the Monte Carlo method to solve BSDEs. Based on random walks, Cheridito and Stadje [14] studied the approximation of BSDEs for non-Lipschitz drivers. Fourier methods to solve FBSDEs were proposed by Huijskens, Ruijter and Oosterlee [40]. A convolution method was used by Hyndman and Ngou [44]. Briand and Labart [8] introduced an algorithm to solve BSDEs based on Wiener chaos and Picard’s iterations expansion. Gobet and Turkedjiev [34] designed a numerical scheme for solving BSDEs by using Malliavin weights and least-squares regression.

With Hermite martingales, the problem of solving an FBSDE is formulated as a system of ordinary differential equations problem by Pelsser and Gnameho [54]. A control variate technique to reduce the simulation error is proposed by Alanko and Avellaneda [1]. For recent approximations of the initial value of high-dimensional BSDEs through Monte Carlo simulation via Picard sequences, see for instance [42, 43] and the references within. Other references can be found in Chassagneux and Richou [11], Khedher and Vanmaele [45], Bender and Steiner [5], Gnameho [30], Zhao, Chen and Peng [63], Ventura and Korzeniowski [58], Gong and Rui [35], among others.

Monte Carlo simulation methods are nowadays at the forefront of learning algorithms. We propose in this paper a new Monte Carlo method which is based on a regression-later approach. Glasserman and Ye [29] show that the regression-later approach offers advantages in comparison to the classical regression technique. Stentohoff [56] discusses the convergence of the regression-now (cf. Glasserman and Yu [29]) technique in the case of the evaluation of an American option.

Numerical solutions for BSDEs have sparked a lot of interest especially in applied mathematics. Our work aims to contribute to the rich body of this literature. Under some regularity assumptions, the solution of an FBSDE can be represented by the solution of a regular quasi-linear parabolic partial differential equation (PDE). By exploiting the Markov property of the solution of the FBSDE, we have developed a probabilistic numerical regression called regression-later algorithm based on the least-squares Monte Carlo method and the connection between the quasi-linear parabolic PDE and the FBSDE in order to provide more accurate solutions. For most numerical algorithms, to solve BSDEs, we need to compute in general two conditional expectations at each step across the time interval. This computation can be very costly especially in high-dimensional problems. Moreover, the computational accuracy deriving the controlled process 𝑍 can be rather unstable even for deep learning methods. The proposed algorithm requires only one conditional expectation to be computed for each step across the time interval. Although like any algorithm built on basis functions, the curse of dimensionality cannot be fully overcome due to the number of basis functions growing with the dimensions, our algorithm yields good convergence and stability results in practice, and the computation of the controlled process 𝑍 is much simpler and easier in comparison to most existing algorithms.

Recently, deep BSDE methods have been investigated, among many others, for instance in [27, 37, 41, 25]. Deep BSDE methods only give a result for the initial single value. In this context, solutions are often susceptible to over-fitting. Furthermore, unfortunately hidden behavior renders the evaluation of the error complex. The methodology developed in our work is a direct approach and based on multivariate polynomial regression. It provides a dynamical solution and does not require any prior training procedure.

This paper is structured as follows. In the first part of our work, we introduce the basic theory of BSDEs, give some general results on the studies of FBSDEs and review the classical backward Euler–Maruyama scheme. In the next step, we describe in detail the regression-later algorithm and derive a convergence result of the scheme. Finally, we provide several numerical experiments to illustrate the performance of the regression-later algorithm in the context of linear and nonlinear FBSDEs.

Notation and assumptions

We will use in this paper the notation of El Karoui, Hamadène and Matoussi [23]. We consider a filtered probability space ( Ω , F , P , F ) with F = F T , F = ( F t ) 0 ≤ t ≤ T a complete natural filtration of a 𝑑-dimensional Brownian motion 𝑊 and 𝑇 a fixed finite horizon. For all m ∈ N * and x ∈ R m , | x | denotes the Euclidean norm of 𝑥. For the matrix A ∈ R m × d , we define its Frobenius norm by | A | := Trace ⁡ ( A ⁢ A * ) . The matrix 𝐴 can be considered as an element of the space R m × d .

Ł m 2 ⁢ ( F t ) := { ( X t ) t ∈ [ 0 , T ] ∈ R m , F t ⁢ -measurable and ⁢ ∥ X ∥ Ł 2 = E ⁢ [ | X t | 2 ] 1 / 2 < ∞ } , S 2 ⁢ ( R m ) := { ( Y t ) t ∈ [ 0 , T ] ∈ R m , continuous and adapted such that ⁢ ∥ Y ∥ S 2 2 = E ⁢ [ sup t ∈ [ 0 , T ] | Y t | 2 ] < ∞ } , H 2 ⁢ ( R m ) := { ( Z t ) t ∈ [ 0 , T ] ∈ R m , continuous and adapted such that ⁢ ∥ Z ∥ H 2 2 = E ⁢ [ ( ∫ 0 T | Z s | 2 ⁢ d s ) ] < ∞ } .

All the equalities and the inequalities between random variables are understood in almost sure sense unless explicitly stated otherwise.
C b l , k ⁢ ( [ 0 , T ] × R m ) is the set of real-valued functions which are 𝑙 times continuously differentiable in their first coordinate and 𝑘 times in their second coordinate with bounded partial derivatives up to order 𝑘.
C k ⁢ ( Ω ) is the set of 𝑘 times continuously differentiable functions on Ω.
For x ∈ R m , ∇ x := ( ∂ ∂ x 1 , … , ∂ ∂ x m ) . The operator ∇ x is called the gradient with respect to 𝑥. In the one-dimensional case, we will use the same notation when there is no confusion.
For x , y ∈ R m , x . y denotes the usual inner product on the space R m .

Equipped with their norm, each of the three spaces Ł m 2 ⁢ ( F t ) , S 2 ⁢ ( R m ) and H 2 ⁢ ( R m ) define a Banach space.

2 Definitions and estimates

In this section, we introduce the general concept of backward stochastic differential equations and forward-backward stochastic differential equations with respect to the standard Brownian motion. In the last part of the section, we recall some classical estimates from the theory of BSDEs.

2.1 Backward stochastic differential equations

In the filtered probability space ( Ω , F , P , F ) , backward stochastic differential equations are a special class of stochastic differential equations. The main difference is that these equations are specified with a prescribed terminal value as shown in the following equation:

(2.1) { − d ⁢ Y t = f ⁢ ( t , Y t , Z t ) ⁢ d ⁢ t − Z t ⁢ d ⁢ W t , 0 ≤ t < T , Y T = ξ .

The preceding system can be written equivalently as

(2.2) Y t = ξ + ∫ t T f ⁢ ( s , Y s , Z s ) ⁢ d s − ∫ t T Z s ⁢ d W s ,

where

𝜉 is the terminal condition of equation (2.1) and is assumed to be an F T -measurable and a square integrable random variable,
the measurable mapping ( t , y , z ) ↦ f ⁢ ( t , y , z ) is generally called the generator or driver of (2.1).

A solution of the backward stochastic differential equation (2.1) is a couple of progressively measurable processes ( Y , Z ) such that

∫ 0 T | Z s | 2 ⁢ d s < ∞ and ∫ 0 T | f ⁢ ( s , Y s , Z s ) | ⁢ d s < ∞ ,
( Y t , Z t ) satisfies equation (2.1).

In general, equation (2.1) does not admit a unique solution. The existence and uniqueness of a solution can be shown under the conditions given by Pardoux and Peng [52], which involves the Lipschitz continuity of the driver function. In that case, we have

( Y t , Z t ) 0 ≤ t ≤ T ∈ S 2 ⁢ ( R m ) × H 2 ⁢ ( R m × d ) .

Remark 2.1

If the generator function is identically equal to zero, BSDE (2.2) is reduced to the following classical stochastic equation:

Y t = ξ − ∫ t T Z s ⁢ d W s .

This previous simplification can be associated with the martingale representation theorem in the filtration generated by the Brownian motion. In this case, the solution 𝑌 is a martingale and we have the explicit formula

Y t = E ⁢ ( ξ | F t ) .

BSDEs appear in numerous problems in finance, in insurance and especially in stochastic control. For example, in a complete financial market, the price of a standard European option can be seen as the solution of a BSDE. The interested reader can consult, among many others, the papers of El Karoui, Peng and Quenez [24], Duffie and Epstein [22], Hamadène and Jeanblanc [36] and the references therein for further details.

2.2 Forward-backward stochastic differential equations

We will consider decoupled forward-backward stochastic differential equations, which consist of a system of two equations given by

(2.3) { X t x = x + ∫ 0 t b ⁢ ( s , X s x ) ⁢ d s + ∫ 0 t σ ⁢ ( s , X s x ) ⁢ d W s , ( t , x ) ∈ [ 0 , T ] × R m , Y t x = ϕ ⁢ ( X T x ) + ∫ t T f ⁢ ( s , X s x , Y s x , Z s x ) ⁢ d s − ∫ t T Z s x ⁢ d W s .

The first component is a forward process and the second a backward process with respect to time. In general, system (2.3) does not admit a unique solution. The existence and uniqueness of a solution can be obtained under the conditions given by Pardoux and Peng [52], which involves the Lipschitz continuity property of the coefficient of system (2.3). We make the following regularity assumptions:

the functions ( t , x ) ↦ b ⁢ ( t , x ) , σ ⁢ ( t , x ) are uniformly Lipschitz in 𝑥 and satisfy
| b ⁢ ( t , x ) | + | σ ⁢ ( t , x ) | ≤ K ⁢ ( 1 + | x | ) ,
there exists a positive constant K > 0 such that
| f ⁢ ( t 1 , x 1 , y 1 , z 1 ) − f ⁢ ( t 2 , x 2 , y 2 , z 2 ) | ≤ K ⁢ ( | x 1 − x 2 | + | y 1 − y 2 | + | z 1 − z 2 | ) for any ⁢ ( t i , x i , y i , z i ) , i = 1 , 2 ,
there exist k σ , K σ > 0 such that, for all t ∈ [ 0 , T ] and x , ζ ∈ R m ,
k σ | ζ | 2 ≤ ∑ i , j [ σ σ * ] i , j ( t , x ) ζ i ζ j | ≤ K σ | ζ | 2 ,

and

there exists a positive constant K > 0 such that | f ⁢ ( t , 0 , 0 , 0 ) | ≤ K for every t ∈ [ 0 , T ] ,
the function x ↦ ϕ ⁢ ( x ) is Lipschitz and belongs to C 1 ⁢ ( R m , R ) , and we denote its Lipschitz constant by C ϕ ,
the driver f : [ 0 , T ] × R m × R × R d → R is continuously differentiable in ( x , y , z ) with uniformly bounded derivatives,
b ∈ C b 0 , 1 ⁢ ( [ 0 , T ] × R m , R m ) and σ ∈ C b 0 , 1 ⁢ ( [ 0 , T ] × R m , R m × d ) .

The last assumption (G4) means that the functions 𝑏 and 𝜎 are continuous in their first coordinate and continuously differentiable in the space variable with uniformly bounded derivatives. The condition ϕ ⁢ ( X T x ) ∈ Ł m 2 ⁢ ( F T ) and assumptions (H1)–(H3) and (G1) ensure the existence and uniqueness of the solution of the decoupled FBSDE (2.3). Based on these assumptions, we also have

( Y t x , Z t x ) 0 ≤ t ≤ T ∈ S 2 ⁢ ( R m ) × H 2 ⁢ ( R m × d ) .

As already mentioned in the introduction, for a large class of FBSDEs, we do not have an explicit solution. We therefore need approximation schemes to solve these equations numerically. Most of the existing numerical schemes are based on the Monte Carlo method. Our regression algorithm is based on the following theorem which establishes a link between the solution of the decoupled FBSDE (2.3) and the solution of the quasi-linear parabolic PDE (2.4). This theorem is one of the cornerstones for our numerical scheme.

Connection between quasi-linear PDE and forward-backward SDE

Theorem 2.2

Theorem 2.2 (Pardoux and Peng [53])

We assume that the following PDE has a classical solution 𝑢, where 𝑢 belongs to the space C 1 , 2 ⁢ ( [ 0 , T ] × R m ) and satisfies, for some constant C > 0 ,

| u ⁢ ( t , x ) | + | ( ∇ x u ) ⁢ ( t , x ) | ≤ C ⁢ ( 1 + | x | q ) ,

and the PDE is given by

(2.4) { ∂ u ∂ t ⁢ ( t , x ) + L ⁢ u ⁢ ( t , x ) + f ⁢ ( t , x , u ⁢ ( t , x ) , ∇ u ⁢ σ ⁢ ( t , x ) ) = 0 , ( t , x ) ∈ [ 0 , T ) × R m , u ⁢ ( T , x ) = ϕ ⁢ ( x ) , x ∈ R m ,

with ℒ being defined by

L ψ = : b ∇ ψ + 1 2 Trace ( A ∇ 2 ψ ) for any ψ ∈ C 1 , 2 ( [ 0 , T ] × R m ) ,

and A = σ ⁢ σ * (the matrix σ * denotes the transpose matrix of 𝜎).

Then the solution of system (2.3) can be represented as

Y t x = u ⁢ ( t , X t x ) and Z t x = σ ⁢ ( t , X t x ) * ⁢ ∇ x u ⁢ ( t , X t x ) for all ⁢ t ∈ [ 0 , T ] .

FBSDEs and their properties are well documented in the literature. Due to their importance, we need robust approximation schemes to solve these equations. Monte Carlo methods are useful tools in this context. Our work will focus on the numerical solution of FBSDE (2.3). We end up this section by providing a key lemma from Zhang [62] which establishes a path regularity result of the martingale integrand Z x . This result is known as the L 2 -time regularity property of Z x . For the reader’s convenience, we recall this result which will often be used in Section 4.2.

Lemma 2.3

Lemma 2.3 (Zhang [62])

Let 𝜋 be a partition of the interval [ 0 , T ] defined as π : 0 = t 0 < t 1 < ⋯ < t N = T , with the mesh Δ i := t i + 1 − t i and | π | := max ⁡ { Δ i ; 0 ≤ i ≤ N − 1 } . Under the assumptions of [62, Theorem 3.1], there exists a positive constant C > 0 independent of 𝜋 such that

∑ i = 1 N E ⁢ ∫ t i − 1 t i | Z s x − Z t i − 1 x | 2 + | Z s x − Z t i x | 2 ⁢ d ⁢ s ≤ C ⁢ ( 1 + | x | 2 ) ⁢ | π | .

3 Implicit backward Euler–Maruyama scheme

In this section, we review the Euler–Maruyama scheme of the forward-backward stochastic differential equation (2.3). As already mentioned in the introduction, there are several algorithms to solve BSDEs numerically. One of the difficulties is to solve a dynamic programming problem which involves the computation of conditional expectations at each step across the time interval. This computation can be very costly especially in high-dimensional problems. For most numerical algorithms, it is important to note that the process Z x is more challenging to compute accurately than the process Y x . Following the work of Gobet, Lemor and Warin [33], let us consider the one-dimensional discrete-time approximation of equation (2.3). We build a partition 𝜋 of the interval [ 0 , T ] defined as π : 0 = t 0 < t 1 < ⋯ < t N = T , with the mesh Δ i := t i + 1 − t i and | π | := max ⁡ { Δ i ; 0 ≤ i ≤ N − 1 } . Let ( X π , Y π , Z π ) be an approximation of the triplet ( X x , Y x , Z x ) defined as follows. The forward component X x of FBSDE (2.3) is approximated by the classical Euler–Maruyama scheme which is given by

{ X 0 π = x , X t i + 1 π = X t i π + Δ i ⁢ b ⁢ ( t i , X t i π ) + σ ⁢ ( t i , X t i π ) ⁢ ( W t i + 1 − W t i ) , 0 < i < N .

By integrating the second equation of system (2.3) from the discretization time t i to t i + 1 ,

Y t i x = Y t i + 1 x + ∫ t i t i + 1 f ⁢ ( s , X s x , Y s x , Z s x ) ⁢ d s − ∫ t i t i + 1 Z s x ⁢ d W s .

An Euler–Maruyama approximation of the previous stochastic integral is defined as

Y t i π = Y t i + 1 π + f ⁢ ( t i , X t i π , Y t i π , Z t i π ) ⁢ Δ i − Z t i π ⁢ Δ ⁢ W t i , Δ ⁢ W t i := W t i + 1 − W t i .

From the previous equalities, Bouchard and Touzi [7] derive the following backward scheme:

(S.I) { Y N π = ϕ ⁢ ( X T π ) , Z t i π = 1 Δ i ⁢ E ⁢ [ Y t i + 1 π ⁢ ( W t i + 1 − W t i ) | F t i ] , 0 ≤ i ≤ N − 1 , Y t i π = E ⁢ [ Y t i + 1 π | F t i ] + Δ i ⁢ f ⁢ ( t i , X t i π , Y t i π , Z t i π ) , 0 ≤ i ≤ N − 1 .

Scheme (S.I) is the standard backward Euler–Maruyama scheme for the backward component of system (2.3). Bouchard and Touzi [7] simulate the conditional expectations using Malliavin calculus techniques. In the spirit of the Longstaff–Schwartz algorithm for American option pricing, Gobet, Lemor and Warin [33] have used regression techniques to approach the solution of scheme (S.I). Their approach is based on the regression-now technique. The numerical scheme (S.I) is widely documented in the literature. The control of the simulation error has been analyzed in several papers. As with many other existing algorithms, the implementation of the previous scheme is not explicit. The work of Gobet, Lemor and Warin [33] provides an implementation of the numerical scheme (S.I) and derives an analytic convergence rate.

Theorem 3.1

Theorem 3.1 (Gobet, Lemor and Warin [33])

Under the assumptions of [33, Theorem 1], there exists a constant C > 0 such that, for | π | small enough,

max 0 ≤ i < N ⁡ E ⁢ | Y t i x − Y t i π | 2 + E ⁢ ∑ i = 0 N − 1 ∫ t i t i + 1 | Z t x − Z t i π | 2 ⁢ d t ≤ C ⁢ ( 1 + | x | 2 ) ⁢ | π | + C ⁢ E ⁢ | ϕ ⁢ ( X T x ) − ϕ ⁢ ( X t N π ) | 2 .

4 Regression-later algorithm

In order to describe our regression-later algorithm, we introduce the equivalent explicit scheme (S.II) below which governs our regression-later algorithm. The regression-later algorithm is devoted to solve numerically the forward-backward stochastic differential equation (2.3). Glasserman and Yu [29] have shown that the regression-later approach offers advantages in comparison to the regression-now technique. For the sake of clarity, we consider the one-dimensional discrete time approximation of system (2.3), where the partition is given by 𝜋. In the new scheme (S.II) below, we denote conventionally by ( X π , Y π , Z π ) an approximation of the triplet ( X x , Y x , Z x ) via our scheme. It is important to note that the family { ( Y π , Z π ) } defined below is different from the one defined in Section 3. Only the forward component X x of system (2.3) is approximated by the same Euler–Maruyama discretization scheme described in the previous Section 3. The other components are obtained as follows: due to the Markov property of our Euler–Maruyama scheme, there exist two measurable deterministic functions u t i π and v t i π such that, for every t i ∈ π , one has Y t i π = u t i π ⁢ ( X t i π ) and Z t i π = v t i π ⁢ ( X t i π ) almost surely. We use then the following scheme:

(S.II) { Y N π = ϕ ⁢ ( X T π ) , Z N π = σ ⁢ ( T , X T π ) ⁢ ( ∇ x ϕ ) ⁢ ( X T π ) , Y t i π = E ⁢ [ Y t i + 1 π | F t i ] + Δ i ⁢ E ⁢ [ f ⁢ ( t i + 1 , X t i + 1 π , Y t i + 1 π , Z t i + 1 π ) | F t i ] , 0 ≤ i ≤ N − 1 , Z t i π = σ ⁢ ( t i , X t i π ) * ⁢ ∇ x Y t i π , 0 ≤ i ≤ N − 1 .

The couple of discrete processes ( Y π , Z π ) is obviously adapted to our filtration by definition. In the sequel, it will be important to control the error of the numerical estimation of the couple ( Y x , Z x ) . On the other hand, the error analysis of the Euler approximation for the forward process X x is well documented and understood.

4.1 Description of the algorithm

We notice that the solution of system (2.3) has its value in an infinite-dimensional space. In order to compute the conditional expectations in our algorithm, for each iteration i ∈ { 0 , … , N } , we define a family ( e j i ) 1 ≤ j ≤ k of truncated orthogonal basis functions of the space Ł m 2 ⁢ ( F t i ) , where ( j , k ) ∈ N * × N * and m = 1 . The integer 𝑘 denotes the number of basis functions. Our algorithm admits five major steps of calculations. We define an orthogonal projection onto the linear subspace generated by the family ( e j i ) 1 ≤ j ≤ k . Each basis function is assumed to be at least differentiable and continuous in the space variable. In our numerical implementation, we will consider a sequence of Hermite polynomials. We start with the same partition 𝜋 of the time interval [ 0 , T ] as in the previous section. We denote by ( Y π , k , Z π , k ) the numerical approximation of the solution on the discretization grid of the partition 𝜋. We also assume that we have the Euler–Maruyama approximation of the forward process 𝑋 at our disposal on the same discretization grid. Moreover, the family of functions ( e j i ) 1 ≤ j ≤ k is selected such that the conditional expectation can be computed exactly. In other words, during the regression-later algorithm below, the conditional expectation term E ⁢ [ e i ⁢ ( X t i + 1 π ) | F t i ] is assumed to be known explicitly via the selected basis functions.

Algorithm 1

Algorithm 1 (Regression-later algorithm)

Initialization: Approximate the terminal condition Y T π , k = Y T π = ϕ ⁢ ( X T π ) .

For i = ( N − 1 ) to 0,

compute the vector α k i + 1 ∈ R k by projection of Y t i + 1 π in (S.II),
{ find α k i + 1 ∈ R k such that J ( α k i + 1 ) = inf α ∈ R k E [ | α . e i ( X t i + 1 π ) − Y t i + 1 π , k | 2 ] , with e i = ( e 1 i ⋮ e k i ) ;
compute Z t i + 1 π , k by the following formal derivation:
Z t i + 1 π , k = α k i + 1 ⁢ ∇ x e i ⁢ ( X t i + 1 π ) ⁢ σ ⁢ ( t i + 1 , X t i + 1 π ) ;
compute the vector β k i + 1 ∈ R k by the following optimization problem:
{ find β k i + 1 ∈ R k such that J ( β k i + 1 ) = inf β ∈ R k E [ | β . e i ( X t i + 1 π ) − f ( t i + 1 , X t i + 1 π , k , Y t i + 1 π , k , Z t i + 1 π , k ) | 2 ] ;
evaluate
Y t i π , k = ( α k i + 1 + β k i + 1 ⁢ Δ i ) . E ⁢ [ e i ⁢ ( X t i + 1 π ) | F t i ] .

End.

The regression-later scheme admits several advantages. The primary advantage is that, at each time step of the algorithm, the scheme requires only the computation of one conditional expectation. A second advantage is that the basis functions ( e j i ) 1 ≤ j ≤ k in the algorithm are selected such that the conditional expectation can be computed exactly. Therefore, the term E ⁢ [ e i ⁢ ( X t i + 1 π ) | F t i ] is known explicitly. These facts decrease significantly the time of computation and accelerate the convergence of the algorithm especially in high-dimensional frameworks where the curse of dimensionality problem occurs. Our numerical implementations yield good convergence results in practice.

4.2 Convergence results

By definition, the couple of discrete processes ( Y π , Z π ) is well defined and adapted to our filtration. Due to the Markov property of scheme (S.II), there exist two measurable deterministic functions u t i π and v t i π such that, for every t i ∈ π , one has Y t i π = u t i π ⁢ ( X t i π ) and Z t i π = v t i π ⁢ ( X t i π ) almost surely. It is important to control the error due to the estimation of the couple ( Y π , Z π ) . This control provides a convergence rate of the regression-later algorithm. Obtaining general convergence results for BSDEs is challenging and can only be done in certain situations for scheme (S.II).

The following theorem provides a convergence rate.

Theorem 4.1

Under assumptions (H1)–(H3) and (G1)–(G4) and inequality (4.15), there exists a positive constant 𝐶 independent of the partition 𝜋 such that, for | π | small enough,

max 0 ≤ i < N ⁡ E ⁢ | Y t i π − Y t i x | 2 + E ⁢ ∑ i = 0 N − 1 ∫ t i t i + 1 | Z s x − Z t i π | 2 ⁢ d s ≤ C ⁢ ( 1 + | x | 2 ) ⁢ | π | + C ⁢ E ⁢ | ϕ ⁢ ( X T x ) − ϕ ⁢ ( X t N π ) | 2 .

Proof

The proof consist of two parts. In the first part, we will prove that

max 0 ≤ i < N ⁡ E ⁢ | Y t i π − Y t i x | 2 ≤ C ⁢ ( 1 + | x | 2 ) ⁢ | π | + C ⁢ E ⁢ | ϕ ⁢ ( X T x ) − ϕ ⁢ ( X t N π ) | 2

and in the second step derive the existence of the constant C > 0 such that

E ⁢ ∑ i = 0 N − 1 ∫ t i t i + 1 | Z s x − Z t i π | 2 ⁢ d s ≤ C ⁢ ( 1 + | x | 2 ) ⁢ | π | + C ⁢ E ⁢ | ϕ ⁢ ( X T x ) − ϕ ⁢ ( X t N π ) | 2 .

During the proof, the constant 𝐶 may take different values from line to line, but will be independent from the partition 𝜋. Let us first remark that, along the time period [ t i , t i + 1 ] ,

(4.1) Y t i x = Y t i + 1 x + ∫ t i t i + 1 f ⁢ ( X s x , Y s x , Z s x ) ⁢ d s − ∫ t i t i + 1 Z s x ⁢ d W s .

Taking the conditional expectation with respect to F t i ,

Y t i x = E ⁢ ( Y t i + 1 x + ∫ t i t i + 1 f ⁢ ( X s x , Y s x , Z s x ) ⁢ d s | F t i ) .

As defined in scheme (S.II), one can compute an approximation of the process Y x at the given time t i as the following:

Y t i π = E ⁢ ( Y t i + 1 π + Δ i ⁢ f ⁢ ( t i + 1 , X t i + 1 π , Y t i + 1 π , Z t i + 1 π ) | F t i ) .

By a backward induction, one can derive from the preceding equality that Y t i π belongs to the space Ł 1 2 ⁢ ( F t i ) . Let us consider

U i = E ⁢ ( | Y t i x − Y t i π | 2 + ∫ t i t i + 1 | Z s x − Z t i π | 2 ⁢ d s ) , 0 ≤ i ≤ N − 1 .

We also define δ ⁢ f i , s π = f ⁢ ( s , X s , Y s x , Z s x ) − f ⁢ ( t i + 1 , X t i + 1 π , Y t i + 1 π , Z t i + 1 π ) , s ∈ [ t i , t i + 1 ] .

Remark 4.2

Note the following.

( Y t i x − Y t i π ) and ( ∫ t i t i + 1 ( Z s x − Z t i π ) ⁢ d W s ) are uncorrelated.
By the martingale representation theorem, there exists an ( F s ) t i ≤ s ≤ t i + 1 -adapted and square integrable process ( Z ̄ s π ) t i ≤ s < t i + 1 such that, for t ∈ [ t i , t i + 1 ] ,
(4.2) Y t π = Y t i + 1 π + ∫ t t i + 1 f ⁢ ( t i + 1 , X t i + 1 π , Y t i + 1 π , Z t i + 1 π ) ⁢ d s − ∫ t t i + 1 Z ̄ s π ⁢ d W s .
The process Z ̄ π is càdlàg and is equal to Z π only on the time instances of the partition 𝜋.

From equations (4.1) and (4.2),

(4.3) Y t i x − Y t i π + ∫ t i t i + 1 ( Z s x − Z t i π ) ⁢ d W s = Y t i + 1 x − Y t i + 1 π + ∫ t i t i + 1 δ ⁢ f i , s π ⁢ d s + ∫ t i t i + 1 ( Z ̄ s π − Z t i π ) ⁢ d W s .

From Lemma A.1, we have, for all a , b , c ∈ R and α > 0 ,

(4.4) ( a + b + c ) 2 ≤ ( 1 + α ) ⁢ a 2 + ( 1 + 2 α ) ⁢ b 2 + ( 1 + α ) ⁢ c 2 + 2 ⁢ a ⁢ c .

Using equation (4.3),

U i = E ⁢ [ Y t i + 1 x − Y t i + 1 π + ∫ t i t i + 1 δ ⁢ f i , s π ⁢ d s + ∫ t i t i + 1 ( Z ̄ s π − Z t i π ) ⁢ d W s ] 2 .

By the Itô isometry formula and the quadratic inequality (4.4), we derive from the above remark that, for every ϵ > 0 ,

U i ≤ E { ( 1 + Δ i ϵ ) | Y t i + 1 π − Y t i + 1 x | 2 + ( 1 + Δ i ϵ ) ∫ t i t i + 1 | Z ̄ s π − Z t i π | 2 d s + ( 1 + 2 ⁢ ϵ Δ i ) ( ∫ t i t i + 1 f ( s , X s x , Y s , Z s ) − f ( t i + 1 , X t i + 1 π , Y t i + 1 π , Z t i + 1 π ) d s ) 2 } + 2 ⁢ E ⁢ { ( Y t i + 1 π − Y t i + 1 x ) ⁢ ( ∫ t i t i + 1 ( Z ̄ s π − Z t i π ) ⁢ d W s ) } .

We know that M t = ∫ t i t ( Z ̄ s π − Z t i π ) ⁢ d W s , t ∈ [ t i , t i + 1 ] , defines a martingale in the Brownian filtration. Plugging equation (4.2) into the last term of the previous inequality and taking the conditional expectation according to F t i and noticing that ( ∫ t t i + 1 δ ⁢ f i , s π ⁢ d s ) t i ≤ t ≤ t i + 1 is of finite variation, we have, from the Itô isometry,

(4.5) U i ≤ E { ( 1 + Δ i ϵ ) | Y t i + 1 π − Y t i + 1 x | 2 + ( 1 + Δ i ϵ ) ∫ t i t i + 1 | Z ̄ s π − Z t i π | 2 d s + ( 1 + 2 ⁢ ϵ Δ i ) ⁢ ( ∫ t i t i + 1 f ⁢ ( s , X s x , Y s x , Z s x ) − f ⁢ ( t i + 1 , X t i + 1 π , Y t i + 1 π , Z t i + 1 π ) ⁢ d ⁢ s ) 2 − 2 E ∫ t i t i + 1 ( Z s x − Z ̄ s π ) ( Z ̄ s π − Z t i π ) d s } .

By noticing that ( Z s x − Z ̄ s π ) ⁢ ( Z ̄ s π − Z t i π ) = ( Z s x − Z t i π ) ⁢ ( Z ̄ s π − Z t i π ) − ( Z ̄ s π − Z t i π ) 2 and from the inequality 2 ⁢ a ⁢ b ≤ 1 θ ⁢ a 2 + θ ⁢ b 2 (with a , b ∈ R for any θ > 0 ), we have, from inequality (4.5),

U i ≤ E { ( 1 + Δ i ϵ ) | Y t i + 1 π − Y t i + 1 x | 2 + ( 3 + Δ i ϵ + θ ) ∫ t i t i + 1 | Z ̄ s π − Z t i π | 2 d s + ( 1 + 2 ⁢ ϵ Δ i ) ( ∫ t i t i + 1 f ( s , X s x , Y s x , Z s x ) − f ( t i + 1 , X t i + 1 π , Y t i + 1 π , Z t i + 1 π ) d s ) 2 + 1 θ E ∫ t i t i + 1 | Z s x − Z t i π | 2 d s } .

By Hölder’s inequality,

U i ≤ E { ( 1 + Δ i ϵ ) | Y t i + 1 π − Y t i + 1 x | 2 + ( 3 + Δ i ϵ + θ ) ∫ t i t i + 1 | Z ̄ s π − Z t i π | 2 d s + ( Δ i + 2 ϵ ) ∫ t i t i + 1 | f ( s , X s x , Y s x , Z s x ) − f ( t i + 1 , X t i + 1 π , Y t i + 1 π , Z t i + 1 π ) | 2 d s + 1 θ E ∫ t i t i + 1 | Z s x − Z t i π | 2 d s } .

By the Lipschitz condition of the driver function 𝑓 and inequality (A.1),

(4.6) U i ≤ E { ( 1 + Δ i ϵ ) | Y t i + 1 π − Y t i + 1 x | 2 + ( 3 + Δ i ϵ + θ ) ∫ t i t i + 1 | Z ̄ s π − Z t i π | 2 d s + 2 ⁢ K 2 ⁢ ( Δ i + 2 ⁢ ϵ ) ⁢ ( ∫ t i t i + 1 | Y s x − Y t i + 1 π | 2 ⁢ d s + ∫ t i t i + 1 | Z s x − Z t i + 1 π | 2 ⁢ d s ) + 1 θ E ∫ t i t i + 1 | Z s x − Z t i π | 2 d s + 2 K 2 ( Δ i + 2 ϵ ) ∫ t i t i + 1 | X s x − X t i + 1 π | 2 d s } .

It is known for instance from [62, Lemma 3.2] or [32, Proposition 5] that there exists a constant C > 0 such that

E ⁢ ∫ t i t i + 1 | X s x − X t i + 1 π | 2 ⁢ d s ≤ 2 ⁢ ∫ t i t i + 1 E ⁢ | X s x − X t i + 1 x | 2 ⁢ d s + 2 ⁢ Δ i ⁢ E ⁢ | X t i + 1 x − X t i + 1 π | 2 ≤ C ⁢ ( 1 + | x | 2 ) ⁢ | π | 2 .

Moreover, we have

By the preceding decomposition, we obtain, from inequality (4.6),

(4.7) U i ≤ E { C ϵ , i K | Y t i + 1 π − Y t i + 1 x | 2 + 4 K 2 ( Δ i + 2 ϵ ) ∫ t i t i + 1 | Z s + Δ i x − Z t i + 1 π | 2 d s + 4 ⁢ K 2 ⁢ ( Δ i + 2 ⁢ ϵ ) ⁢ ∫ t i t i + 1 | Y s x − Y t i + 1 x | 2 ⁢ d s + ( 3 + Δ i ϵ + θ ) ⁢ ∫ t i t i + 1 | Z ̄ s π − Z t i π | 2 ⁢ d s + 4 K 2 ( Δ i + 2 ϵ ) ∫ t i t i + 1 | Z s x − Z s + Δ i x | 2 d s + 1 θ E ∫ t i t i + 1 | Z s x − Z t i π | 2 d s } + C ϵ ( 1 + | x | 2 ) | π | 2 ,

where C ϵ , i K = ( 1 + Δ i / ϵ + 4 ⁢ Δ i ⁢ K 2 ⁢ ( Δ i + 2 ⁢ ϵ ) ) and C ϵ = 4 ⁢ C ⁢ K 2 ⁢ ( Δ i + 2 ⁢ ϵ ) . Clearly,

(4.8) E ⁢ ∫ t i t i + 1 | Z s + Δ i x − Z t i + 1 π | 2 ⁢ d s = E ⁢ ∫ t i + 1 t i + 2 | Z s x − Z t i + 1 π | 2 ⁢ d s .

From equality (4.8), inequality (4.7) becomes

(4.9) U i ≤ E { C ϵ , i K | Y t i + 1 π − Y t i + 1 x | 2 + 4 K 2 ( Δ i + 2 ϵ ) ∫ t i + 1 t i + 2 | Z s x − Z t i + 1 π | 2 d s + C ϵ ⁢ ( 1 + | x | 2 ) ⁢ | π | 2 + 4 ⁢ K 2 ⁢ ( Δ i + 2 ⁢ ϵ ) ⁢ ∫ t i t i + 1 | Y s x − Y t i + 1 x | 2 ⁢ d s + 4 K 2 ( Δ i + ϵ ) ∫ t i t i + 1 | Z s x − Z s + Δ i x | 2 d s + ( 3 + Δ i ϵ + θ ) ∫ t i t i + 1 | Z ̄ s π − Z t i π | 2 d s + 1 θ ∫ t i t i + 1 | Z s x − Z t i π | 2 d s } .

By [62, Lemma 3.2], there exists a positive constant 𝐶 such that

(4.10) E ⁢ ∫ t i t i + 1 | Y s x − Y t i + 1 x | 2 ⁢ d s ≤ C ⁢ ( 1 + | x | 2 ) ⁢ | π | 2 .

Inserting inequality (4.10) into (4.9) and setting ( ϵ ; 1 θ ) = ( 1 16 ⁢ K 2 ; 1 2 ) , we derive a constant C > 0 such that

U ̃ i ≤ ( 1 + C ⁢ Δ i ) ⁢ U ̃ i + 1 + C ⁢ ( 1 + Δ i ) ⁢ ( ( 1 + | x | 2 ) ⁢ | π | 2 + E ⁢ ∫ t i t i + 1 | Z s x − Z s + Δ i x | 2 ⁢ d s ) + C ⁢ ( 1 + Δ i ) ⁢ E ⁢ ∫ t i t i + 1 | Z ̄ s π − Z t i π | 2 ⁢ d s ,

where U ̃ i = U i − 1 2 ⁢ E ⁢ ∫ t i t i + 1 | Z s x − Z t i π | 2 ⁢ d s . From Lemma 2.3 and Lemma A.2, we infer that there exists a constant C > 0 such that, for | π | small enough,

(4.11) max 0 ≤ i ≤ N U ̃ i ≤ C E ( ϕ ( X T x ) − ϕ ( X T π ) 2 + C ∑ i = 0 N − 1 E ∫ t i t i + 1 | Z ̄ s π − Z t i π | 2 d s + C ( 1 + | x | 2 ) | π | .

The following argument concludes our proof. Interval-by-interval, given the fact that Z ̄ t i π = Z t i π , by Lemma 4.3 and Lemma 2.3, there exists a positive constant C > 0 independent of 𝜋 such that

(4.12) ∑ i = 0 N − 1 E ⁢ ∫ t i t i + 1 | Z ̄ s π − Z t i π | 2 ⁢ d s ≤ C ⁢ ( 1 + | x | 2 ) ⁢ | π | + C ⁢ E ⁢ ( ϕ ⁢ ( X T x ) − ϕ ⁢ ( X T π ) ) 2 .

Inserting inequality (4.12) into inequality (4.11), we obtain

(4.13) max 0 ≤ i ≤ N ⁡ U ̃ i ≤ C ⁢ E ⁢ ( ϕ ⁢ ( X T x ) − ϕ ⁢ ( X T π ) ) 2 + C ⁢ ( 1 + | x | 2 ) ⁢ | π | .

In particular, one can derive the following inequality which completes the first step of the proof of the theorem:

(4.14) max 0 ≤ i ≤ N ⁡ E ⁢ | Y t i π − Y t i x | 2 ≤ C ⁢ E ⁢ ( ϕ ⁢ ( X T x ) − ϕ ⁢ ( X T π ) ) 2 + C ⁢ ( 1 + | x | 2 ) ⁢ | π | .

From inequality (4.10) and Lemma 2.3, inequality (4.9) becomes, for | π | small enough, choosing ( ϵ , 1 θ ) = ( 1 32 ⁢ K 2 , 1 2 ) ,

U ̃ i − 1 + 1 4 ⁢ E ⁢ ∫ t i t i + 1 | Z s x − Z t i π | 2 ⁢ d s ≤ ( 1 + C ⁢ Δ i ) ⁢ U ̃ i + C ⁢ ( 1 + | x | 2 ) ⁢ | π | 2 + C ⁢ ( 1 + Δ i ) ⁢ E ⁢ ∫ t i − 1 t i | Z ̄ s π − Z t i − 1 π | 2 ⁢ d s ,

where C > 0 and we know that

E ⁢ | Y t i π − Y t i x | 2 = U ̃ i − 1 2 ⁢ E ⁢ ∫ t i t i + 1 | Z s x − Z t i π | 2 ⁢ d s .

Summing both sides of the previous inequality over the variable 𝑖 from 1 to N − 1 , and using inequality (4.12), there exists a positive constant C > 0 independent of 𝜋 such that

∑ i = 1 N − 1 U ̃ i − 1 + 1 4 ⁢ E ⁢ ∑ i = 1 N − 1 ∫ t i t i + 1 | Z s x − Z t i π | 2 ⁢ d s ≤ ∑ i = 1 N − 1 ( 1 + C ⁢ Δ i ) ⁢ U ̃ i + C ⁢ ( 1 + | x | 2 ) ⁢ | π | .

We deduce from the previous relation and inequality (4.13) that there exists a constant C > 0 independent of 𝜋 such that

∑ i = 0 N − 1 E ⁢ ∫ t i t i + 1 | Z s x − Z t i π | 2 ⁢ d s ≤ C ⁢ ( 1 + | x | 2 ) ⁢ | π | + C ⁢ E ⁢ | ϕ ⁢ ( X T x ) − ϕ ⁢ ( X t N π ) | 2 .

The previous inequality and (4.14) conclude. ∎

Lemma 4.3

We assume that there exists a positive constant γ > 0 such that

(4.15) ( y − y ′ ) . ( δ ⁢ f 1 , 2 ⁢ ( t ) ) ≤ γ ⁢ ( | y − y ′ | ⁢ | x 1 − x 2 | + | y − y ′ | ⁢ | y 1 − y 2 | ) ,

where δ ⁢ f 1 , 2 ⁢ ( t ) = f ⁢ ( t , x 1 , y 1 , z 1 ) − f ⁢ ( t , x 2 , y 2 , z 2 ) . Under assumptions (H1)–(H3) and (G1)–(G4), the function u t i π is uniformly Lipschitz for any t i ∈ π with the Lipschitz constant independent of the partition 𝜋.

Proof

During the proof, the constant 𝐶 may take different values from line to line, but is independent from the partition 𝜋. We consider the same partition 𝜋 of the interval [ 0 , T ] as described in algorithm (S.II). We recall that Y t i π defines the Euler approximation of Y t i x (the exact process at the time step t i ). As introduced previously,

Y t i π = E ⁢ ( Y t i + 1 π + Δ i ⁢ f ⁢ ( t i + 1 , X t i + 1 π , Y t i + 1 π , Z t i + 1 π ) | F t i ) .

By the martingale representation theorem, there exists an ( F s ) t i ≤ s ≤ t i + 1 -adapted and square integrable process ( Z ̄ s π ) t i ≤ s ≤ t i + 1 such that

(4.16) Y t π = Y t i + 1 π + ∫ t t i + 1 f ⁢ ( t i + 1 , X t i + 1 π , Y t i + 1 π , Z t i + 1 π ) ⁢ d s − ∫ t t i + 1 Z ̄ s π ⁢ d W s , t i ≤ t ≤ t i + 1 .

The preceding representation can be seen as a continuous version of a BSDE on the time interval [ t i , t i + 1 ] . Let us introduce the continuous Euler discretization of the process of X x in system (2.3) given by

(4.17) X s π = X t i π + ∫ t i s b ⁢ ( t i , X t i π ) ⁢ d u + ∫ t i s σ ⁢ ( t i , X t i π ) ⁢ d W u , s ∈ [ t i , t i + 1 ) .

We consider X π , i ( i = 1 , 2 ) two solutions of (4.17) associated with two initial conditions x i ( i = 1 , 2 ). To each X π , i , we associate the couple ( Y π , x i , Z ̄ π , x i ) solution of equation (4.16). Let us define

Δ ⁢ Y t 1 , 2 := Y t π , x 1 − Y t π , x 2 , Δ ⁢ Z t 1 , 2 := Z t π , x 1 − Z t π , x 2 , Δ ⁢ X t 1 , 2 := X t π , 1 − X t π , 2 and Δ ⁢ Z ̄ t 1 , 2 := Z ̄ t π , x 1 − Z ̄ t π , x 2 .

As highlighted above, due to the Markov property of our Euler scheme, there exist two measurable deterministic functions u t i π and v t i π such that, for every t i ∈ π , one has Y t i π = u t i π ⁢ ( X t i π ) and Z t i π = v t i π ⁢ ( X t i π ) almost surely. It is known that u t i is locally Lipschitz (see e.g. El Karoui, Hamadène and Matoussi [23]). We will then show that u t i is Lipschitz by a backward induction.

For i = N , the function x ↦ u T ⁢ ( x ) = ϕ ⁢ ( x ) is Lipschitz by assumption. We now suppose that the function u t i + 1 π is Lipschitz in the space variable with C i + 1 its Lipschitz constant. We will show that u t i π is Lipschitz.

Applying Itô’s formula to | Y π , x 1 − Y π , x 2 | 2 and taking the expectation, we obtain

where

δ ⁢ f i π = f ⁢ ( t i + 1 , X t i + 1 π , x 1 , Y t i + 1 π , x 1 , Z t i + 1 π , x 1 ) − f ⁢ ( t i + 1 , X t i + 1 π , x 2 , Y t i + 1 π , x 2 , Z t i + 1 π , x 2 ) .

From assumption (H2), the inequality a ⁢ b ≤ 1 2 ⁢ α ⁢ a 2 + 1 2 ⁢ α ⁢ b 2 , for every α > 0 and assumption (4.15),

Setting α = 1 γ and plugging the last inequality into (4.18), we have in particular

(4.19) E ⁢ | Δ ⁢ Y t 1 , 2 | 2 ≤ ( 1 + Δ i ) ⁢ E ⁢ | Δ ⁢ Y t i + 1 1 , 2 | 2 + Δ i ⁢ E ⁢ | Δ ⁢ X t i + 1 1 , 2 | 2 + 2 ⁢ γ 2 ⁢ ∫ t t i + 1 E ⁢ | Δ ⁢ Y s 1 , 2 | 2 ⁢ d s .

During our backward induction, we have assumed above that the function u t i + 1 π is Lipschitz. From equation (4.17), we have

E ⁢ | X t i + 1 π , x 1 − X t i + 1 π , x 2 | 2 ≤ ( 1 + C ⁢ Δ i ) ⁢ | x 1 − x 2 | 2 .

From (4.19), Gronwall’s inequality from Lemma A.3 applied to the function t ↦ E ⁢ | Δ ⁢ Y t 1 , 2 | 2 with t ∈ [ t i , t i + 1 ) leads to

E ⁢ | Δ ⁢ Y t 1 , 2 | 2 ≤ ( ( 1 + Δ i ) ⁢ ( 1 + C ⁢ Δ i ) ⁢ C i + 1 2 + Δ i ⁢ ( 1 + C ⁢ Δ i ) ) ⁢ exp ⁡ ( 2 ⁢ γ 2 ⁢ Δ i ) ⁢ | x 1 − x 2 | 2 .

We recall that our objective is to prove that the function u t i π is Lipschitz with a uniform Lipschitz constant in the space variable. We have

| u t i π ⁢ ( x 1 ) − u t i π ⁢ ( x 2 ) | 2 ≤ C i 2 ⁢ | x 1 − x 2 | 2 ,

where C i 2 = ( ( 1 + Δ i ) ⁢ ( 1 + C ⁢ Δ i ) ⁢ C i + 1 2 + Δ i ⁢ ( 1 + C ⁢ Δ i ) ) ⁢ exp ⁡ ( 2 ⁢ γ 2 ⁢ Δ i ) . It is then enough to show that the positive constant C i is uniformly bounded to conclude the backward induction result. Let us first remark that, in the neighborhood of zero, there exists a positive constant 𝐶 such that exp ⁡ ( x ) ≤ ( 1 + C ⁢ x ) . Hence, for Δ i small enough, there exists a positive constant 𝐶 such that C i 2 ≤ ( 1 + C ⁢ Δ i ) ⁢ C i + 1 2 + C ⁢ Δ i . By Lemma A.2, we have the following uniformly bounded inequality:

max 0 ≤ i ≤ N ⁡ C i 2 ≤ e C ⁢ T ⁢ ( C ϕ 2 + C ⁢ T ) ,

where C ϕ is the Lipschitz constant of the function 𝜙 in the forward-backward stochastic differential equation (2.3). Finally,

| u t i π ⁢ ( x 1 ) − u t i π ⁢ ( x 2 ) | 2 ≤ e C ⁢ T ⁢ ( C ϕ 2 + C ⁢ T ) ⁢ | x 1 − x 2 | 2 .

This completes the induction. From the previous inequality, the function x ↦ u t i π ⁢ ( x ) is Lipschitz with a uniform Lipschitz constant e C ⁢ T ⁢ ( C ϕ 2 + C ⁢ T ) . ∎

5 Applications

In this section, we provide three numerical experiments to illustrate the performance of the regression-later algorithm: the first in the context of option pricing, the second in the case where the terminal condition is a functional of a Brownian motion and the last on the backward Burgers equation. The first example is connected to the numerical approximation of a linear BSDE, while the second one analyzes nonlinear BSDEs.

As mentioned in the introduction, BSDEs appear in numerous problems in finance, in insurance and especially in stochastic control. A frequent problem in finance or in insurance is the valuation of a contract or the risk management of a portfolio (see e.g. [50, 26, 51, 31, 60, 20], among others). Linear and nonlinear BSDEs appear naturally in these situations. The interested reader can consult the papers of El Karoui, Peng and Quenez [24], Delong [21], Duffie and Epstein [22], Hamadène and Jeanblanc [36] and the references therein for further details. We will implement in the first example the linear case and investigate how fast our algorithm converges. In financial markets, the most popular contracts of derivative securities are European and American call and put options.

In our first example, we therefore evaluate standard European options. The algorithm can also be applied to compute the price of some non-path dependent insurance contracts. In our implementation, we consider the orthogonal Hermite polynomial family as basis in order to solve analytically the conditional expectations occurring at each iteration.

Application 1: Pricing

Our market model is composed of two financial assets: 𝑆 (risky asset) and S 0 (risk-less asset). Let

{ S t ⁢ is the price of ⁢ S ⁢ at time ⁢ t , S t 0 ⁢ is the price of ⁢ S 0 ⁢ at time ⁢ t .

Based on their assumptions, Black and Scholes have modelled the dynamic of the risky asset 𝑆 as a geometric Brownian motion. We denote by the constant 𝑟 the daily interest rate which is assumed to be constant. The process S 0 is governed by the following differential equation: d ⁢ S t 0 = r ⁢ S t 0 ⁢ d ⁢ t with the initial condition S 0 0 = 1 . The process S t follows the following linear SDE with constant coefficients:

{ d ⁢ S t S t = μ ⁢ d ⁢ t + σ ⁢ d ⁢ W t , S 0 = x ,

where μ ∈ R is a constant drift coefficient which represents the expected rate of return of 𝑆, S 0 is the initial value of the risky asset 𝑆 and 𝜎 is a constant positive volatility coefficient. Let us consider a European call option on the risky asset 𝑆 with characteristics ( K , T ) , where 𝑇 is the maturity date and 𝐾 is the strike value of the contract. The seller of the call option is committed to pay to the holder the sum ( S T − K ) + which represents the profit that allows to exercise the option. We build the following portfolio: at the time instance 𝑡, we invest a Δ t part of the risky asset and a β t part of the non-risky asset. Denoting 𝑌 the wealth process, we have, at time 𝑡,

Y t = Δ t ⁢ S t + β t ⁢ S t 0 .

A main assumption is that our strategy is self-financing, and in a context of continuously trading for the agent, a mathematical translation is given by

d ⁢ Y t = Δ t ⁢ d ⁢ S t + β t ⁢ d ⁢ S t 0 .

Denoting θ = μ − r σ and Z t = σ ⁢ Δ t ⁢ S t , the triplet ( S t , Y t , Z t ) solves

(E₁) { − d ⁢ Y t = f ⁢ ( t , S t , Y t , Z t ) ⁢ d ⁢ t − Z t ⁢ d ⁢ W t , Y T = ϕ ⁢ ( S T ) , d ⁢ S t = μ ⁢ S t ⁢ d ⁢ t + σ ⁢ S t ⁢ d ⁢ W t , S 0 = x 0 ,

where ϕ ⁢ ( x ) = ( x − K ) + and f ⁢ ( t , x , y , z ) = − ( r ⁢ y + θ ⁢ z ) . The value of the portfolio at 𝑡 is given by Y t and Z t is related to the hedging strategy. We can evaluate explicitly the value of the wealth process 𝑌 for a fixed time. In particular, at the time instance t = 0 , Y 0 = E ⁢ ( e − r ⁢ T ⁢ ϕ ⁢ ( S T ) ⁢ exp ⁡ ( − θ ⁢ W T + 1 2 ⁢ θ 2 ⁢ T ) ) . By evaluating the preceding expectation, we obtain the classical Black–Scholes formula

Y 0 = e − r ⁢ T ⁢ ( F T ⁢ N ⁢ ( d + ) − K ⁢ N ⁢ ( d − ) ) and Z 0 = σ ⁢ N ⁢ ( d + ) ⁢ S 0 ,

where

F T = S 0 ⁢ e r ⁢ T , d ± = log ⁡ ( F T / K ) ± 1 2 ⁢ σ 2 ⁢ T σ ⁢ T and N ⁢ ( x ) = 1 2 ⁢ π ⁢ ∫ − ∞ x e − t 2 2 ⁢ d t .

Moreover, Z 0 is the portfolio value at t = 0 . The function 𝒩 denotes the cumulative distribution function of the standard normal distribution. The objective is to provide a numerical solution of system (E₁ ). We are interested in the initial value of the couple ( Y , Z ) . We suppose that we have at our disposal the value of the forward process 𝑆 on the grid of the partition 𝜋.

In our numerical simulation, we consider a finite-dimensional system of normalized orthogonal Hermite polynomials. We choose a fixed number of basis functions and evaluate the couple ( Y , Z ) along the time period [ 0 , T ] .

We build the partition 𝜋 of the interval [ 0 , T ] defined as π : 0 = t 0 < ⋯ < t N = T , with Δ i := t i + 1 − t i and | π | := max ⁡ { Δ i ; 0 ≤ i ≤ N − 1 } . We set the following parameters:

𝑘 is the number of basis functions,
𝑀 is the number of simulated paths of the Brownian motion,
𝑁 is the number of the discretization points on 𝜋.

As input values, we use

T = 1 , r = 0.02 , S 0 = 100 , K = 100 , μ = 1 % , σ = 20 % .

Ignoring the domestic currency’s unity, the value of the European call option is Y 0 = 8.916 and Z 0 = 11.5852 . Figure 1 shows the logarithm of the relative error when estimating the couple ( Y 0 , Z 0 ) . Modulo the choice of | π | , 𝑘 and 𝑁, the error decreases significantly as we increase the number of simulations 𝑀. The error curves on the estimation of Z 0 seem to be more volatile in this particular case.

In the case of a European put option, the same argument as above leads to a similar conclusion and similar convergence results. The regression-later approach shows a stable convergence trend which is less volatile in comparison to scheme (S.I). In Monte Carlo frameworks, a few number of basis functions is sufficient to obtain convergence of the algorithm due to the exponential decay of Hermite coefficients to zero.

$Figure 1 Log absolute error on ( Y 0 , Z 0 ) (Y_{0},Z_{0}) , European call$

Figure 1

Log absolute error on ( Y 0 , Z 0 ) , European call

$Figure 2 Log absolute error curve to estimate ( Y 0 , Z 0 ) (Y_{0},Z_{0})$

Figure 2

Log absolute error curve to estimate ( Y 0 , Z 0 )

Application 2: Brownian functional case

In this example, we implement the regression-later algorithm in the context of nonlinear BSDEs. The underlying process is assumed to be a standard Brownian motion 𝑊 on the time interval [ 0 , 1 ] . The terminal function is a functional of the Brownian motion 𝑊. We consider the following BSDE:

(5.1) { − d ⁢ Y t = f ⁢ ( t , W t , Y t , Z t ) ⁢ d ⁢ t − Z t ⁢ d ⁢ W t , 0 ≤ t < 1 , Y 1 = ϕ ⁢ ( W 1 ) ,

where, for ( t , x , y , z ) ∈ [ 0 , 1 ] × R × R × R ,

ϕ ⁢ ( x ) = x ⁢ arctan ⁡ ( x ) − ln ⁡ ( 1 + x 2 ) and f ⁢ ( t , x , y , z ) = − 1 2 ⁢ ( 1 + tan 2 ⁡ ( z ) ) .

The exact solution of the above system is given by

( Y t , Z t ) = ( − 1 2 ⁢ ln ⁡ ( 1 + W t 2 ) + W t ⁢ arctan ⁡ ( W t ) , arctan ⁡ ( W t ) ) .

By noting that the function x ↦ ln ⁡ ( x ) satisfies the linear growth condition and the function x ↦ arctan ⁡ ( x ) is bounded, the unique solution of (5.1) satisfies

( Y t , Z t ) 0 ≤ t ≤ T ∈ S 2 ⁢ ( R ) × H 2 ⁢ ( R ) .

The exact value of the couple ( Y , Z ) at the time point t 0 is ( Y 0 , Z 0 ) = ( 0 , 0 ) . Figure 2 shows the empirical logarithm of the absolute error when computing ( Y 0 , Z 0 ) .

Modulo the choice of | π | and the number of the basis function 𝑘, the graphics show a stable convergence result. This leads us to the same conclusion as above. Nevertheless, the estimation of the initial value Y 0 is more stable and quicker than the estimation of the initial value Z 0 in the first example.

Application 3: A high-dimensional case

For numerical computation in high dimensions, we consider an example based on the backward Burgers equation (e.g. [9, 28, 18], among others). In practice, one of the difficulties is to construct a suitable system of basis functions in high dimensions especially in anisotropic spaces. Solutions in high dimensions are more costly numerically in general. Multivariate Hermite polynomials appear naturally when approximating multidimensional functions or distributions and will therefore serve as a system of basis functions. Hermite polynomials are used in the Schrödinger theory of quantum physics. These polynomials appear also naturally in the study of the propagation of the heat equation, in the study of quantum harmonic oscillator, etc. In high dimensions, the definition of these polynomials can simply be done through multivariate normal distributions. Based on the inputs of our underlying example, we will use the well-known tensorized construction. In this case, in a 𝑑-dimensional space, generalized Hermite polynomials for any multi-index n = ( n 1 , … , n d ) ∈ N d are defined by

H n ⁢ ( x ) = ∏ i = 1 d H n i ⁢ ( x i ) , with ⁢ x ∈ R d ,

where

H n i ⁢ ( x ) = ( − 1 ) n i ⁢ e x 2 / 2 ⁢ d n i d ⁢ x n i ⁢ e − x 2 / 2 , x ∈ R , n i ∈ N , n i ≥ 1 and H 0 ⁢ ( x ) = 1 .

The above equality defines the one-dimensional probabilist’s Hermite polynomials using Rodrigues’s formula. The components of the sequence ( H n ⁢ ( x ) ) n ∈ N are orthogonal polynomials with respect to the Gaussian weight function

ϕ ⁢ ( x ) = exp ⁡ { − 1 2 ⁢ x 2 } / 2 ⁢ π , x ∈ R .

Hence, for ( n , m ) ∈ N 2 , any pair ( H n ⁢ ( x ) , H m ⁢ ( x ) ) satisfies the orthogonality relationship

∫ − ∞ ∞ H n ⁢ ( x ) ⁢ H m ⁢ ( x ) ⁢ ϕ ⁢ ( x ) ⁢ d x = n ! ⁢ δ n ⁢ m .

A simple expression for multivariate Hermite polynomials is given by Withers [59]. Other methods for computing multivariate Hermite polynomials exist (see e.g. [4, 19, 38], among others). Differentiation of these multivariate polynomials is based on the chain rule argument and easy due to the independence structure of H n and the chain rule argument. We will consider an example based on the backward Burgers equation in high dimensions. Burgers’ equation is a fundamental PDE appearing in various areas of mathematics (see e.g. [9, 28, 18], among others). In this example, the underlying forward process is driven by a standard three-dimensional Brownian motion 𝑊 on the time interval [ 0 , T ] . As in the one-dimensional case, the approximation is performed sequentially by solving a dynamical programming problem. The numerical computation in high dimensions is the same as the one-dimensional case but computationally more expensive. Our backward Burgers’ equation is governed by the following equations:

{ − d ⁢ Y t = f ⁢ ( t , S t , Y t , Z t ) ⁢ d ⁢ t − Z t ⁢ d ⁢ W t , Y T = ϕ ⁢ ( S T ) and 0 ≤ t < T , d ⁢ S t = σ ⁢ d ⁢ W t , S 0 = 0 and σ > 0 ,

where the function f : [ 0 , T ] × R d × R × R d ↦ R is the driver defined by

f ⁢ ( t , x , y , z ) = ( y − 2 + d 2 ⁢ d ) ⁢ ∑ i = 1 d z i ,

and x ↦ ϕ ⁢ ( x ) = exp ⁡ ( T + 1 d ⁢ ∑ i = 1 d x i ) / ( 1 + exp ⁡ ( T + 1 d ⁢ ∑ i = 1 d x i ) ) , x ∈ R d .

The analytical solution of the above system of equations is given by

( Y t , Z t ) = ( exp ⁡ ( t + 1 d ⁢ ∑ i = 1 d S t i ) ( 1 + exp ⁡ ( t + 1 d ⁢ ∑ i = 1 d S t i ) ) , σ ⁢ exp ⁡ ( t + 1 d ⁢ ∑ i = 1 d S t i ) d ⁢ ( 1 + exp ⁡ ( t + 1 d ⁢ ∑ i = 1 d S t i ) ) 2 ⁢ 1 d ) ,

where 1 d ∈ R d takes the value one at each coordinate.

The algorithm could be summarized by solving sequentially a system of linear equations in high dimensions. The implementation is performed with the normalized system of the generalized Hermite polynomials by truncation. As the dimension of the problem could be very high, this rises naturally the cost of computation.

Nevertheless, in our framework, only a few number of basis functions have sizable influence on the solution and are sufficient to obtain convergence of the algorithm. This is a consequence of the exponential decay of Hermite coefficients when solving the projection problems by least-squares technique at each iteration.

The computation of the conditional exception in high dimensions could be unstable and challenging. The computation of 𝑍 involves also the computation of a system of gradient terms. In our algorithm, the conditional exception and the gradient term at each iteration are known explicitly which decreases significantly the time of computation and thus accelerates the speed of convergence of the algorithm especially in high-dimensional frameworks where the curse of dimensionality problem occurs. As pointed out by Xiu and Karniadakis [61], Hermite expansion has been quite effective in solving or computing the functional of Markov processes with Gaussian inputs. This holds true as well as for certain types of non-Gaussian inputs.

For low-dimensional implementation of generalized Hermite polynomials, one could use the package Hermite.py in Python software. For higher dimensions, it is straightforward to implement these polynomials directly (e.g. [55, 57, 46], among others). As input values for our numerical experiments, we define the following parameters:

T = 1 , d = 3 , σ = 1 , k = 56 .

The integer 𝑘 is the number of selected basis functions. Usually, in practice, we retain all polynomials with the total degree less than or equal to certain degree which in our example is equal to 5. We then have 56 basis functions using the Cartesian product. In Monte Carlo frameworks, a few number of basis functions is sufficient to obtain convergence of the algorithm due to the exponential decay of the Hermite coefficients to zero. The exact value of 𝑌 at t 0 = 0 is Y t 0 = 0.5 . The coordinates of the analytical solution of Z t 0 are given by

Z t 0 = ( Z t 0 1 Z t 0 2 Z t 0 3 ) = 1 12 ⁢ ( 1 1 1 ) .

$Figure 3 Log-error curve to estimate Y t 0 Y_{t_{0}}$

Figure 3

Log-error curve to estimate Y t 0

Figure 3 shows the empirical logarithm of the absolute error of the truncated solution induced by the numerical approximation of Y t 0 . The graphic shows a clear convergence of the estimated Y t 0 .

$Figure 4 Log-error curve to estimate Z ̂ t 0 = 1 3 ⁢ ( Z t 0 1 + Z t 0 2 + Z t 0 3 ) \hat{Z}_{t_{0}}=\frac{1}{3}(Z^{1}_{t_{0}}+Z^{2}_{t_{0}}+Z^{3}_{t_{0}})$

Figure 4

Log-error curve to estimate Z ̂ t 0 = 1 3 ⁢ ( Z t 0 1 + Z t 0 2 + Z t 0 3 )

Figure 4 shows the empirical logarithm of the absolute error of the truncated solution induced by the numerical approximation of Z ̂ t 0 = 1 3 ⁢ ( Z t 0 1 + Z t 0 2 + Z t 0 3 ) . As in previous cases, our numerical results show a stable convergence regarding the estimation of the solution ( Y , Z ) . Nevertheless, the estimation of Y t 0 is slower than the estimation of Z ̂ t 0 . The error curve on the estimation of Z ̂ t 0 seems to be more volatile but more precise. The computation time is low, making the algorithm potentially suitable for the class of high-dimensional Markovian problems.

6 Conclusion

In this paper, we have discussed a new numerical scheme for backward stochastic differential equations (BSDEs). The scheme is based on the regression-later approach. In the first part of our work, we introduced the theory of BSDEs, gave some general background on their studies and reviewed the classical backward Euler–Maruyama scheme. In the next step, we described our regression-later algorithm in detail and derived a convergence result. Finally, we provided numerical experiments to illustrate the performance of the regression-later algorithm: a first experiment in the context of option pricing, a second in the case where the terminal condition is a functional of a Brownian motion and a third based on the backward Burgers equation. Solutions in high dimensions are numerically more costly. Multivariate Hermite polynomials appear naturally when approximating multidimensional functions or distributions. Modulo a suitable choice of the number of discretization points and the number of basis functions, our numerical results show a stable convergence. In many numerical algorithms for solving BSDEs, one of the difficulties is to solve a dynamic programming problem which involves often the computation of conditional expectations at each step of the time interval. We remark that, in many alternative algorithms, the numerical computation of 𝑍 is more challenging than the computation of the process 𝑌, leading to potential numerical instabilities especially in higher dimensions. It is interesting to note that our algorithm circumvents this difficulty, by obtaining the numerical approximation of 𝑍 directly from the approximation of 𝑌.

A Appendix

Lemma A.1

For any constant α > 0 and for any a , b ∈ R ,

(A.1) ( a + b ) 2 ≤ ( 1 + α ) ⁢ a 2 + ( 1 + 1 α ) ⁢ b 2 .

Proof

The result is a direct consequence of Young’s inequality. ∎

Let us now recall the classical discrete Gronwall lemma (see e.g. [62, 49]).

Lemma A.2

Lemma A.2 (Gronwall Inequality A)

Let us consider the partition π : 0 = t 0 < ⋯ < t N = T of the interval [ 0 , T ] and let Δ i be its mesh. We also consider the families ( a k ) 0 ≤ k ≤ N and ( b k ) 0 ≤ k ≤ N , assumed to be non-negative such that, for some positive constant γ > 0 , we have a k − 1 ≤ ( 1 + γ ⁢ Δ i ) ⁢ a k + b k , k = 1 , … , N . Then

max 0 ≤ i ≤ N ⁡ a i ≤ e γ ⁢ T ⁢ ( a N + ∑ i = 1 N b i ) .

Lemma A.3

Lemma A.3 (Gronwall Inequality B)

Let y , b , a : [ 0 , T ] ↦ R be three continuous functions such that 𝑏 is non-negative and

y ⁢ ( t ) ≤ a ⁢ ( t ) + ∫ 0 t b ⁢ ( s ) ⁢ y ⁢ ( s ) ⁢ d s , 0 ≤ t ≤ T .

Then

y ⁢ ( t ) ≤ a ⁢ ( t ) + ∫ 0 t a ⁢ ( s ) ⁢ b ⁢ ( s ) ⁢ exp ⁡ ( ∫ s t b ⁢ ( u ) ⁢ d u ) ⁢ d s , 0 ≤ t ≤ T .

In addition, if the function 𝑎 is non-decreasing, then

y ⁢ ( t ) ≤ a ⁢ ( t ) ⁢ exp ⁡ ( ∫ 0 t b ⁢ ( s ) ⁢ d s ) , 0 ≤ t ≤ T .

References

[1] S. Alanko and M. Avellaneda, Reducing variance in the numerical solution of BSDEs, C. R. Math. Acad. Sci. Paris 351 (2013), no. 3–4, 135–138. 10.1016/j.crma.2013.02.010Search in Google Scholar

[2] V. Bally, Approximation scheme for solutions of BSDE, Backward Stochastic Differential Equations, Pitman Res. Notes Math. Ser. 364, Longman, Harlow (1997), 177–191. Search in Google Scholar

[3] G. Barles and E. Lesigne, SDE, BSDE and PDE, Backward Stochastic Differential Equations, Pitman Res. Notes Math. Ser. 364, Longman, Harlow (1997), 47–80. Search in Google Scholar

[4] O. E. Barndorff-Nielsen, Asymptotic techniques; for use in statistics, Technical report, 1989. 10.1007/978-1-4899-3424-6Search in Google Scholar

[5] C. Bender and J. Steiner, A posteriori estimates for backward SDEs, SIAM/ASA J. Uncertain. Quantif. 1 (2013), no. 1, 139–163. 10.1137/120878689Search in Google Scholar

[6] J.-M. Bismut, Conjugate convex functions in optimal stochastic control, J. Math. Anal. Appl. 44 (1973), 384–404. 10.1016/0022-247X(73)90066-8Search in Google Scholar

[7] B. Bouchard and N. Touzi, Discrete-time approximation and Monte-Carlo simulation of backward stochastic differential equations, Stochastic Process. Appl. 111 (2004), no. 2, 175–206. 10.1016/j.spa.2004.01.001Search in Google Scholar

[8] P. Briand and C. Labart, Simulation of BSDEs by Wiener chaos expansion, Ann. Appl. Probab. 24 (2014), no. 3, 1129–1171. 10.1214/13-AAP943Search in Google Scholar

[9] J. M. Burgers, A mathematical model illustrating the theory of turbulence, Advances in Applied Mechanics, Academic Press, New York (1948), 171–199. 10.1016/S0065-2156(08)70100-5Search in Google Scholar

[10] J.-F. Chassagneux and C. A. Garcia Trillos, Cubature method to solve BSDEs: Error expansion and complexity control, Math. Comp. 89 (2020), no. 324, 1895–1932. 10.1090/mcom/3522Search in Google Scholar

[11] J.-F. Chassagneux and A. Richou, Numerical simulation of quadratic BSDEs, Ann. Appl. Probab. 26 (2016), no. 1, 262–304. 10.1214/14-AAP1090Search in Google Scholar

[12] P. E. Chaudru de Raynal and C. A. Garcia Trillos, A cubature based algorithm to solve decoupled McKean–Vlasov forward-backward stochastic differential equations, Stochastic Process. Appl. 125 (2015), no. 6, 2206–2255. 10.1016/j.spa.2014.11.018Search in Google Scholar

[13] P. Cheridito, H. M. Soner, N. Touzi and N. Victoir, Second-order backward stochastic differential equations and fully nonlinear parabolic PDEs, Comm. Pure Appl. Math. 60 (2007), no. 7, 1081–1110. 10.1002/cpa.20168Search in Google Scholar

[14] P. Cheridito and M. Stadje, BSΔEs and BSDEs with non-Lipschitz drivers: Comparison, convergence and robustness, Bernoulli 19 (2013), no. 3, 1047–1085. 10.3150/12-BEJ445Search in Google Scholar

[15] D. Chevance, Résolution numérique des équations différentielles stochastiques rétrogrades, PhD thesis, Universtité de Provence, 1997. Search in Google Scholar

[16] D. Crisan and K. Manolarakis, Solving backward stochastic differential equations using the cubature method: Application to nonlinear pricing, SIAM J. Financial Math. 3 (2012), no. 1, 534–571. 10.1137/090765766Search in Google Scholar

[17] D. Crisan and K. Manolarakis, Second order discretization of backward SDEs and simulation with the cubature method, Ann. Appl. Probab. 24 (2014), no. 2, 652–678. 10.1214/13-AAP932Search in Google Scholar

[18] A. B. Cruzeiro and E. Shamarova, On a forward-backward stochastic system associated to the Burgers equation, Stochastic Analysis with Financial Applications, Progr. Probab. 65, Birkhäuser/Springer, Basel (2011), 43–59. 10.1007/978-3-0348-0097-6_4Search in Google Scholar

[19] G. Dattoli, S. Lorenzutta, G. Maino and A. Torre, Ann. Numer. Math. 2 (1995), no. 1–4, 211–232. Search in Google Scholar

[20] G. Deelstra, P. Devolder, K. Gnameho and P. Hieber, Valuation of hybrid financial and actuarial products in life insurance by a novel three-step method, Astin Bull. 50 (2020), no. 3, 709–742. 10.1017/asb.2020.25Search in Google Scholar

[21] L. Delong, Backward Stochastic Differential Equations with Jumps and Their Actuarial and Financial Applications, Springer, London, 2013. 10.1007/978-1-4471-5331-3Search in Google Scholar

[22] D. Duffie and L. G. Epstein, Stochastic differential utility, Econometrica 60 (1992), no. 2, 353–394. 10.2307/2951600Search in Google Scholar

[23] N. El Karoui, S. Hamadène and A. Matoussi, Backward stochastic differential equations and applications, Indifference Pricing Theory Appl. 8 (2008), 267–320. 10.1515/9781400833115.267Search in Google Scholar

[24] N. El Karoui, S. Peng and M. C. Quenez, Backward stochastic differential equations in finance, Math. Finance 7 (1997), no. 1, 1–71. 10.1111/1467-9965.00022Search in Google Scholar

[25] M. Fujii, A. Takahashi and M. Takahashi, Asymptotic expansion as prior knowledge in deep learning method for high dimensional BSDES, Asia-Pacific Financial Markets 26 (2019), no. 3, 391–408. 10.1007/s10690-019-09271-7Search in Google Scholar

[26] H. Geman, N. El Karoui and J.-C. Rochet, Changes of numéraire, changes of probability measure and option pricing, J. Appl. Probab. 32 (1995), no. 2, 443–458. 10.2307/3215299Search in Google Scholar

[27] M. Germain, H. Pham and X. Warin, Approximation error analysis of some deep backward schemes for nonlinear PDEs, SIAM J. Sci. Comput. 44 (2022), no. 1, 28–56. 10.1137/20M1355355Search in Google Scholar

[28] O. Glass and S. Guerrero, On the uniform controllability of the Burgers equation, SIAM J. Control Optim. 46 (2007), no. 4, 1211–1238. 10.1137/060664677Search in Google Scholar

[29] P. Glasserman and B. Yu, Simulation for American options: Regression now or regression later?, Monte Carlo and Quasi-Monte Carlo methods 2002, Springer, Berlin (2004), 213–226. 10.1007/978-3-642-18743-8_12Search in Google Scholar

[30] K. Gnameho, Numerical solution of backward stochastic differential equations and volatility risk premium, PhD thesis, Maastricht University, 2016. Search in Google Scholar

[31] K. Gnameho, J. Kanniainen and Y. Yue, Modeling variance risk premium, Mathematical and Statistical Methods for Actuarial Sciences and Finance, Springer, Cham (2017), 129–141. 10.1007/978-3-319-50234-2_11Search in Google Scholar

[32] E. Gobet and C. Labart, Error expansion for the discretization of backward stochastic differential equations, Stochastic Process. Appl. 117 (2007), no. 7, 803–829. 10.1016/j.spa.2006.10.007Search in Google Scholar

[33] E. Gobet, J.-P. Lemor and X. Warin, A regression-based Monte Carlo method to solve backward stochastic differential equations, Ann. Appl. Probab. 15 (2005), no. 3, 2172–2202. 10.1214/105051605000000412Search in Google Scholar

[34] E. Gobet and P. Turkedjiev, Approximation of backward stochastic differential equations using Malliavin weights and least-squares regression, Bernoulli 22 (2016), no. 1, 530–562. 10.3150/14-BEJ667Search in Google Scholar

[35] B. Gong and H. Rui, One order numerical scheme for forward-backward stochastic differential equations, Appl. Math. Comput. 271 (2015), 220–231. 10.1016/j.amc.2015.08.127Search in Google Scholar

[36] S. Hamadène and M. Jeanblanc, On the starting and stopping problem: Application in reversible investments, Math. Oper. Res. 32 (2007), no. 1, 182–192. 10.1287/moor.1060.0228Search in Google Scholar

[37] J. Han and J. Long, Convergence of the deep BSDE method for coupled FBSDEs, Probab. Uncertain. Quant. Risk 5 (2020), Paper No. 5. 10.1186/s41546-020-00047-wSearch in Google Scholar

[38] B. Holmquist, The 𝑑-variate vector Hermite polynomial of order 𝑘, Linear Algebra Appl. 237 (1996), 155–190. 10.1016/0024-3795(95)00595-1Search in Google Scholar

[39] Y. Hu, P. Imkeller and M. Müller, Utility maximization in incomplete markets, Ann. Appl. Probab. 15 (2005), no. 3, 1691–1712. 10.1214/105051605000000188Search in Google Scholar

[40] T. P. Huijskens, M. J. Ruijter and C. W. Oosterlee, Efficient numerical Fourier methods for coupled forward-backward SDEs, J. Comput. Appl. Math. 296 (2016), 593–612. 10.1016/j.cam.2015.10.019Search in Google Scholar

[41] C. Huré, H. Pham and X. Warin, Deep backward schemes for high-dimensional nonlinear PDEs, Math. Comp. 89 (2020), no. 324, 1547–1579. 10.1090/mcom/3514Search in Google Scholar

[42] M. Hutzenthaler, A. Jentzen, T. Kruse and T. A. Nguyen, Overcoming the curse of dimensionality in the numerical approximation of backward stochastic differential equations, J. Numer. Math. 31 (2023), no. 1, 1–28. Search in Google Scholar

[43] M. Hutzenthaler, T. Kruse and T. A. Nguyen, On the speed of convergence of Picard iterations of backward stochastic differential equations, Probab. Uncertain. Quant. Risk 7 (2022), no. 2, 133–150. 10.3934/puqr.2022009Search in Google Scholar

[44] C. B. Hyndman and P. O. Ngou, Global convergence and stability of a convolution method for numerical solution of bsdes, preprint (2014), https://arxiv.org/abs/1410.8595. Search in Google Scholar

[45] A. Khedher and M. Vanmaele, Discretisation of FBSDEs driven by càdlàg martingales, J. Math. Anal. Appl. 435 (2016), no. 1, 508–531. 10.1016/j.jmaa.2015.10.022Search in Google Scholar

[46] T. Kollo and D. von Rosen, Advanced Multivariate Statistics with Matrices, Math. Appl. (New York) 579, Springer, Dordrecht, 2005. 10.1007/1-4020-3419-9Search in Google Scholar

[47] R. J. A. Laeven and M. Stadje, Robust portfolio choice and indifference valuation, Math. Oper. Res. 39 (2014), no. 4, 1109–1141. 10.1287/moor.2014.0646Search in Google Scholar

[48] J. Ma, P. Protter and J. M. Yong, Solving forward-backward stochastic differential equations explicitly—a four step scheme, Probab. Theory Related Fields 98 (1994), no. 3, 339–359. 10.1007/BF01192258Search in Google Scholar

[49] J. Mémin, S. Peng and M. Xu, Convergence of solutions of discrete reflected backward SDE’s and simulations, Acta Math. Appl. Sin. Engl. Ser. 24 (2008), no. 1, 1–18. 10.1007/s10255-006-6005-6Search in Google Scholar

[50] R. C. Merton, Applications of option-pricing theory: Twenty-five years later, Amer. Econ. Rev. 88 (1998), no. 3, 323–349. Search in Google Scholar

[51] T. Møller, On valuation and risk management at the interface of insurance and finance, British Actuarial J. 8 (2002), no. 4, 787–827. 10.1017/S1357321700003913Search in Google Scholar

[52] E. Pardoux and S. Peng, Adapted solution of a backward stochastic differential equation, Systems Control Lett. 14 (1990), no. 1, 55–61. 10.1016/0167-6911(90)90082-6Search in Google Scholar

[53] E. Pardoux and S. Peng, Backward stochastic differential equations and quasilinear parabolic partial differential equations, Stochastic Partial Differential Equations and Their Applications, Lect. Notes Control Inf. Sci. 176, Springer, Berlin (1992), 200–217. 10.1007/BFb0007334Search in Google Scholar

[54] A. Pelsser and K. Gnameho, A Monte Carlo method for backward stochastic differential equations with Hermite martingales, Monte Carlo Methods Appl. 25 (2019), no. 1, 37–60. 10.1515/mcma-2019-2028Search in Google Scholar

[55] B. Ribeiro, R. F. Albrecht, A. Dobnikar, D. W. Pearson and N. C. Steele, Adaptive and Natural Computing Algorithms, Springer, Vienna, 2005. 10.1007/b138998Search in Google Scholar

[56] L. Stentoft, Convergence of the least squares Monte Carlo approach to American option valuation, Manag. Sci. 50 (2004), no. 9, 1193–1203. 10.1287/mnsc.1030.0155Search in Google Scholar

[57] B. Sudret, Global sensitivity analysis using polynomial chaos expansions, Reliab. Eng. Syst. Safety 93 (2008), no. 7, 964–979. 10.1016/j.ress.2007.04.002Search in Google Scholar

[58] W. A. Ventura and A. Korzeniowski, On discretely reflected backward stochastic differential equations, Stoch. Anal. Appl. 34 (2016), no. 1, 1–23. 10.1080/07362994.2015.1094670Search in Google Scholar

[59] C. S. Withers, A simple expression for the multivariate Hermite polynomials, Statist. Probab. Lett. 47 (2000), no. 2, 165–169. 10.1016/S0167-7152(99)00153-4Search in Google Scholar

[60] M. V. Wüthrich and M. Merz, Financial Modeling, Actuarial Valuation and Solvency in Insurance, Springer Finance, Springer, Heidelberg, 2013. 10.1007/978-3-642-31392-9Search in Google Scholar

[61] D. Xiu and G. E. Karniadakis, Modeling uncertainty in flow simulations via generalized polynomial chaos, J. Comput. Phys. 187 (2003), no. 1, 137–167. 10.1016/S0021-9991(03)00092-5Search in Google Scholar

[62] J. Zhang, A numerical scheme for BSDEs, Ann. Appl. Probab. 14 (2004), no. 1, 459–488. 10.1214/aoap/1075828058Search in Google Scholar

[63] W. Zhao, L. Chen and S. Peng, A new kind of accurate numerical method for backward stochastic differential equations, SIAM J. Sci. Comput. 28 (2006), no. 4, 1563–1581. 10.1137/05063341XSearch in Google Scholar

Received: 2023-08-29

Revised: 2024-01-20

Accepted: 2024-01-29

Published Online: 2024-02-14

Published in Print: 2024-06-01

This work is licensed under the Creative Commons Attribution 4.0 International License.

Articles in the same Issue

https://doi.org/10.1515/mcma-2024-2002

Keywords for this article

Regression; BSDE; Monte Carlo; Hermite polynomials

Creative Commons

BY 4.0