Radon Series Comp. Appl. Math 8, 1–34

c de Gruyter 2009

Discrete-time approximation of BSDEs andprobabilistic schemes for fully nonlinear PDEs
Bruno Bouchard, Romuald Elie, and Nizar Touzi

Abstract. The aim of this paper is to provide a survey on recent advances on probabilistic numericalmethods for nonlinear PDEs, which serve as an alternative to classical deterministic schemes andallow to handle a large class of multidimensional nonlinear problems. These probabilistic schemes
are based on the stochastic representation of semilinear PDEs by means of backward SDEs, which
can be viewed as an extension of the well-known Feynman–Kac formula to the semilinear case. In
this context, we first explain how smoothness properties can be obtained for non-reflected BSDEs set
in the whole domain by using purely probabilistic techniques introduced in Ma and Zhang [59], andhow they can be exploited to provide convergence rates of discrete time Euler-type approximation
schemes. We then present some recent extensions to BSDEs with jumps, reflected BSDEs, BSDEs
with horizon given by a finite stopping time, and second order BSDEs. The extension to fully nonlinear PDEs requires to enlarge the definition of backward SDEs. However, a natural probabilistic
numerical scheme can be introduced formally by evaluating the nonlinear PDE along the trajectoryof some chosen underlying diffusion. This point of view shows an intimate connection between theprobabilistic numerical schemes reviewed in this paper and the standard finite differences methods.
In particular, convergence and error estimates are provided by using the monotone schemes methods
from the theory of viscosity solutions.
Key words. Backward stochastic differential equations, numerical resolution, conditional expectation, monte carlo methods.
AMS classification. 65C05, 60H35, 60G20

1Introduction

The application field offered by the resolution of BSDEs is wide and diverse. As noticed by Bismut [12], Backward Stochastic Differential Equations (hereafter BSDEs)
cover in particular the field of stochastic optimal control problems and the corresponding Hamilton–Jacobi–Bellman equations. See also [67], [72] for optimal control of diffusions with jumps, and [15] and [51] for optimal switching problems. Consequently,
these techniques are of great interest for applications in economics. This is particularly
true in finance, where their use was popularised by El Karoui, Peng and Quenez [48].
In incomplete markets, it can be used to solve optimal management problems by working along the lines presented in, e.g., El Kaouri, Peng and Quenez [49] for recursive
utilities, El Karoui and Rouge [50] for indifference pricing, see also the more recent
Third author: This research is part of the Chair Financial Risks of the Risk Foundation sponsored by SociétéGénérale, the Chair Derivatives of the Future sponsored by the Fédération Bancaire Française, the Chair Financeand Sustainable Development sponsored by EDF and Calyon.

2

B. Bouchard and R. Elie and N. Touzi

paper [46] and the references therein, and Peng [68] for risk measures. Adding jumps
to financial assets dynamics can be done as in [10] or [36]. BSDEs with jumps arealso used for the valuation of financial derivatives with default risk. The valuation ofAmerican or game options follows from the resolution of reflected BSDEs, see [47].
The list of potential applications is long.
More generally, it is well known that the so-called Feynman Kac representation extends to second order semilinear PDEs which can be interpreted in terms of decoupled
Forward-Backward SDEs, see Pardoux and Peng [63, 64] and Pardoux [62]. The term
“decoupled” refers to the fact that the dynamic of the forward process does not dependon the solution of the Backward SDE. Hence, solving a BSDE or the related semilinear
PDE is essentially the same. It was actually first proposed by Ma, Protter and Yong[58], see also [34], to use an estimation of the solution of the semilinear PDE to provide
an approximation of the solution of Forward-Backward SDEs. However, solving PDEs
in high dimension by finite-difference or finite-element schemes is a challenging task
and any deterministic procedure suffers from the so-called “curse of dimensionality”
which, very quickly, prevents from a good approximation of the solution of the PDE.
As a powerful alternative to classical purely deterministic techniques, more probabilistic approaches, directly based on the resolution of BSDEs, have been introduced
these last years. We can essentially classify them in two categories. The first one concerns the approximation of the driving Brownian motion by a random walk on a finite
tree. Such techniques have been widely studied, see e.g. [1], [28], [29] or [57], however they can be efficient only in low dimension since the size of the tree explodes
exponentially with the dimension. The second one is based on Monte-Carlo type techniques. The main idea consists in simulating the Forward process by using standard
Euler schemes and then computing backward the solution of the BSDE associated tothe simulated forward paths by: 1. Discretising in time the BSDE and writing an Eulertype backward scheme which involves conditional expectations, 2. Approximating theconditional expectations by using non-parametric regression methods, see [22], [56],[27], [43] or integration by parts based estimators, see [55], [17] and [20]. The main
advantage of such methods is that, like other Monte-Carlo type methods, the convergence rate does not depend a-priori on the dimension of the problem and they should
therefore less suffer from the curse of dimensionality. In practice, the variance and/or
the convergence rate of the estimators of the conditional expectations usually depends
on the complexity/dimension of the problem, and the above comment has to be tempered. There is actually a third class of methods which lies in between the two previous
ones and which were developed in a series of papers by Bally, Pagès and Printems, seee.g. [3], [4] and the references therein. In this approach, the original forward process
is replaced by a quantised version in discrete time, i.e. a projection on a finite grid,
whose associated Backward SDE can be computed by backward induction. The difference with the random walk technique is that the grid is computed in some optimal way,
which allows to consider high dimension problems with a limited number of points,
see the above quoted papers.
We shall not detail the pure numerical parts in this paper but rather insists on thefirst step of all these algorithms: the discrete time approximation of the BSDE whichleads to a Backward Euler-type scheme.

3

Discrete-time approximation of BSDEs

Not surprisingly, the discrete time approximation error depends heavily on the regularity of the solution itself. The good news is that it essentially only depends on it.
Therefore, all the analysis breaks down to that of the regularity of the solution of the
BSDE. This has been first observed by Zhang [73, 74] and Bouchard and Touzi [20],who exploited the regularity results of Ma and Zhang [59] to provide a |π|1/2 converge
rate for decoupled Forward-Backward SDEs with Lipschitz coefficients, for a timegrid modulus |π|. Various extensions have been then developed in the literature and arereviewed in this paper: reflected BSDEs, the presence of a jump component, and thecorresponding Dirichlet problem.
More recently, the backward schemes suggested by the discretisation of BSDEs havealso been extended to the case of fully nonlinear PDEs by Fahim, Touzi and Warin[37]. The connection with an extended class of BSDEs, namely second order BSDEs,
was outlined in [25]. The last section of this paper provides a presentation of theprobabilistic numerical schemes which avoids BSDEs and only requires to evaluate thenonlinear PDE along the trajectory of some chosen underlying diffusion. In particular,
this point of view highlights the relationship of the probabilistic schemes studied inthis paper with the corresponding standard deterministic schemes. The convergence
results and the corresponding error estimates in the fully nonlinear case are reviewed
in the final section of this paper, and are based on the monotone schemes method
initiated within the theory of viscosity solutions by Barles and Souganidis [9], Krylov[52, 53, 54], and Barles and Jakobsen [6, 7, 8].
To conclude, we should note that the case of coupled Forward-Backward SDEs has
been considered in e.g. [31], [32], see also the references therein. The nature of theanalysis being very different from the ones carried out here, we shall not discuss thiscase here (although it is important). See also [11] for Picard type iterations techniques.
Notations: Any element of Rd is viewed as a column vector and we denote by ′ thetransposition. Md (resp. Sd ) denotes the set of d × d (resp. symmetric) matrices. Fora smooth function (t, x) 7→ ϕ(t, x), we denote by ∂t ϕ its first derivative with respect
to t, and by Dϕ and D2 ϕ its Jacobian and Hessian matrix with respect to x. It thefunction depends on more than one space variable, we use Dx , Dy , etc . . . to denotethe partial Jacobian with respect to a specific argument.

2

Numerical approximation of decoupled forward
backward SDE

Let us first consider the simplest case of decoupled Forward-Backward SDE on thewhole domain:
Xt

= x+

Z t
0

Yt

b(Xs )ds −= g(X1 ) +

Z 1
t

Z t

σ(Xs )dWs

(2.1)

0

h(Xs , Ys , Zs )ds −Z 1
t

Zs′ dWs ,

t≤1

(2.2)

4

B. Bouchard and R. Elie and N. Touzi

where W is d-dimensional Brownian Motion on some probability space (Ω, F , P),
whose associated complete and right-continuous filtration is denoted by F = (Ft )t≤1 .
Here, b : Rd → Rd , σ : Rd → Md , g : Rd → R and h : Rd × R × Rd → R are some
Lipschitz continuous functions.

2.1 Discrete time approximation
The discrete-time approximation is constructed as follows.
Euler scheme for the Forward process. Given a time-grid π := {0 = t0 < t1 < · · · <tn = 1} of [0, 1] with modulus |π|, we first approximate the Forward process X by its
Euler scheme X π defined as
X0π := x , and

Xtπi +1 := b(Xtπi )∆ti + σ(Xtπi )∆Wti for i < n ,(2.3)

where ∆ti := ti+1 − ti and ∆Wti := Wti+1 − Wti .
Euler scheme for the Backward process. As for the solution (Y, Z) of the BSDE, we
first approximate the terminal condition Y1 = g(X1 ) by simply replacing X by its Euler
scheme : Y1 ≃ g(X1π ). Then, motivated by the formal Euler discretisation:
Z ti+1
Z ti+1Zs′ dWs
h(Xs , Ys , Zs )ds −=Yti+1 +
Yti
ti

ti

′′

≃′′

Yti+1 + h(Xtπi , Yti , Zti )∆ti − Zt′i ∆Wti ,we define the discrete-time approximation of the BSDE as follows:
1. Taking expectation conditionally to Fti on both sides leads to
Yti ≃ E Yti+1 | Fti + h(Xtπi , Yti , Zti )∆ti .
2. Pre-multiplying by ∆Wti and then taking conditional provides
0 ≃ E Yti+1 ∆Wti | Fti − Zti ∆ti .

This formal approximation argument leads to a Backward Euler scheme (Ytπi , Z̄tπi ) of
the formih1π Z̄tπ =
|
F
∆W
E
Y
ti
t
ti+1
i
∆thi
i i
(2.4)

 Y π = E Y π | Fti + h X π , Y π , Z̄ π ∆ti
t
t
t
t
t
i

i+1

i

i

i

with terminal condition given by Y1π = g(X1π ). The reason why we denote the Z component by Z̄ π , and not simply Z π , will become clear in Remark 2.1.2 below.
Remark 2.1.1 (Explicit or implicit schemes). Note that the above scheme is implicit
as Ytπi appears in both sides of the equation. Since h is assumed to be Lipschitz andsince it is multiplied by ∆ti , which is intended to be small, the equation can however
be solved numerically very quickly by standard fixed point methods. Nevertheless, weDiscrete-time approximation of BSDEs

5

could also consider an explicit version of this scheme by replacing the second equationin (2.4) byi
hYtπi = E Ytπi+1 + h Xtπi , Ytπi+1 , Z̄tπi ∆ti | Fti .

This will not change the convergence rate, see [18]. Both can thus be used indifferently
depending on the exact nature of the coefficients.

Remark 2.1.2 (Continuous time version of the Euler scheme). Note that the martingale representation theorem implies that we can find a square integrable process Z π
such that
Z ti+1πππππ(Zsπ )′ dWs .
Yti = Yti+1 + h Xti , Yti , Z̄ti ∆ti −tiThis allows to consider a continuous version of the Backward Euler scheme by setting
Ytπ := Ytπi+1 +Z ti+1th

Xtπi , Ytπi , Z̄tπi−Z ti+1t(Zsπ )′ dWs for t ∈ [ti , ti+1 ) .

(2.5)

It then follows from the Itô isometry thatZ̄tπi =

1E∆ti

Z ti+1ti

Zsπ ds | Fti



.

(2.6)

In other words, the step constant process Z̄ π , defined as Z̄tπ := Z̄tπi for t ∈ [ti , ti+1 ),can be seen as the best L2 (Ω × [0, 1], dP ⊗ dt) approximation of Z π in the class ofadapted processes which are constant on each interval [ti , ti+1 ).

2.2 Error decomposition
In order to compare (Y, Z) and (Y π , Z π ), let us introduceZ̄ti :=

1E∆ti

Z ti+1ti

Zs ds | Fti

and observe that, since h is Lipschitz continuous,−dYt=h(Xt , Yt , Zt )dt − Zt′ dWth(Xtπ , Yti , Z̄ti )dt − Zt′ dWt

for t ∈ [ti , ti+1 )(2.7)
up to an error term of order C |Xt − Xtπi | + |Yt − Yti | + |Zt − Z̄ti | dt for some constant C > 0. Next observe from the definition of Z̄ti and Z̄tπi and Jensen’s inequalitythat
Z ti+1
1E |Z̄ti − Z̄tπi |2 ≤ E|Zt − Ztπ |2 dt .
(2.8)
∆ti ti≃6

B. Bouchard and R. Elie and N. Touzi

It then follows from (2.5)–(2.7), recalling the involved error terms, and standard techniques for BSDEs with Lipschitz coefficients that
Z ti+1π 2
π 2|Zt − Zt | dt
max E |Yti − Yti | +
0≤i<n

≤ CE

ti

"

sup

ti ≤t≤ti+1|Xt − Xtπi |2 +2supti ≤t≤ti+1|Yt − Yti | +Z ti+1ti

2|Zt − Z̄ti | dt#

for some C > 0. Using (2.8) again and standard estimates for Forward Euler schemes,
it follows that##"n−1 Z"
X ti+1π 2
π 22|Zt − Z̄ti | dt
Err(π) := max E
sup |Yt − Yti | + E
0≤i<n

≤

ti ≤t≤ti+1

2 max E
0≤i<n

"

sup
ti ≤t≤ti+1|Zt − Z̄ti |

:=

"

ti

i=0

≤

|Yti − Ytπi |2 + |Yt − Yti |2"n−1 Z
X ti+1

+2 E

ti

i=0

2

#


+ |Z̄ti − Z̄tπi |2 dt#!


C |π| + R2Y (π) + R2Z (π)

for some C > 0, where

R2Y (π)R2Z (π)

:=

max E

0≤i<n

E

sup

ti ≤t≤ti+1"n−1 Z
X ti+1i=0ti

2|Yt − Yti |2|Zt − Z̄ti | dt

#

and

#

can be seen as (squared) modulus of regularity of Y and Z .
Proposition 2.2.1 ([20], [74]). There is a constant C > 0, which is independent on πsuch that:1
Err(π) ≤ C |π| + R2Y (π) + R2Z (π) 2 .

Remark 2.2.2. 1. Note that, by continuity of the process Y , the quantity supti ≤t≤ti+1
P − a.s. to 0 as |π| → 0. Since the coefficients are Lipschitz, we|Yt − Yti |2 convergesalso have E supt≤1 |Yt |2 < ∞. It thus follows from dominated convergence that
R2Y (π) −→ 0 as |π| → 0 .

2. Also note that the step constant process Z̄ , defined as Z̄t := Z̄ti for t ∈ [ti , ti+1 ), canbe seen as the best L2 (Ω × [0, 1], dP ⊗ dt) approximation of Z in the class of adapted
processes which are constant on each interval [ti , ti+1 ). Since such processes can beused to approximate Z in L2 (Ω × [0, 1], dP ⊗ dt), it follows that
R2Z (π) −→ 0 as |π| → 0 .

7Discrete-time approximation of BSDEs

Combining Proposition 2.2.1 and the last remark, implies the first converge result.
Corollary 2.2.3. We have: Err(π)2 −→ 0 as |π| → 0 .
Remark 2.2.4. It follows from the interpretation of Z̄ in terms of the best approximation of Z by step constant processes that"n−1 Z#
X ti+1π 2E|Zt − Z̄ti | dt ≥ R2Z (π) .
i=0ti

In the case where h ≡ 0, then Y is a martingale and Yti is therefore the best L2 (Ω, dP)
approximation of Yt for t ∈ [ti , ti+1 ] by a Fti -measurable random variable. The
Burckolder–Davis–Gundy’s inequality then implies that
"
#
"
#
E

sup

ti ≤t≤ti+1|Yt − Ytπi |2 ≥ c E

supti ≤t≤ti+1|Yt − Yti |2,

for some c > 0. In view of Proposition 2.2.1, this shows thatc R2Y (π) + R2Z (π) ≤ Err(π)2 ≤ C |π| + R2Y (π) + R2Z (π) .

Up to a |π| term this can actually be shown even if h 6≡ 0. This result is not surprising. As usual, the convergence rate should depend on the regularity of the functions/processes to be approximated.

2.3 Path regularity and convergence rate
In view of Proposition 2.2.1, the whole analysis of the convergence rate of the discretetime approximation breaks down to the analysis of the (squared) modulus of continuity
R2Y (π) and R2Z (π). The regularity of Y is standard in our context where the coefficientsare Lipschitz continuous.
Proposition 2.3.1. There exists C > 0 independent on π such that R2Y (π) ≤ C |π|.
Note that this implies some regularity on the (unique viscosity) solution v of thesemilinear PDE−Lv − h(·, v, σ ′ Dv)= 0 on [0, 1) × Rdv(T −, ·) = g on Rd ,where L is the Dynkin operator associated to X ,1 
Lv := ∂t v + hb, Dvi + Tr σσ ′ D2 v ,2since it is related to (Y, Z) by the mapping Y = v(·, X) and (at least formally) Z =σ(X)′ Dv(·, X), see e.g. [64]. As a matter of fact, it follows from our Lipschitz continuity assumption that this solution is Lipschitz continuous in space and 1/2-Hölder8

B. Bouchard and R. Elie and N. Touzi

in time. Since X is itself 1/2-Hölder in time in L2 , this actually implies the regularity
statement of Proposition 2.3.1.
It is however less standard to study the regularity of the Z -component. If Dv , whenever it is well-defined, were Lipschitz continuous in space and 1/2-Hölder in time, wewould easily deduce that R2Z (π) is of order |π|. But such a result is not true in general,and indeed much stronger than what we need.
A pure probabilistic analysis of the modulus of regularity of Z has been carried out
by Ma and Zhang [59] and Zhang [74], actually in a more general setting. It is basedon the following analysis (d = 1 to avoid notational complexities):
Step 1. The first step consists in assuming that g and σ are Cb1 and using the fact that
(Y, Z) are differentiable in the Malliavin sense and that there exists a version of their
Malliavin derivative (DY, DZ) which solves
Z 1
Z 1
Ds Yt = Dg(X1 )Ds X1 +Dh(Θr )Ds (Xr , Yr , Zr )dr −Ds Zr dWr , s ≤ 1 ,ttwhere Θ := (X, Y, Z), see e.g. [61]. Moreover, Itô’s formula and uniqueness of solutions of BSDEs with Lipschitz coefficients imply that it is related to the solution(∇Y, ∇Z) of
Z 1
Z 1∇Yt = Dg(X1 )∇X1 +Dh(Θr )∇(Xr , Yr , Zr )dr −∇Zr dWr , s ≤ 1 ,ttwhere ∇X denotes the first variation process of X , by the relation Ds (X, Y, Z) =1[s,1] ∇(X, Y, Z)(∇Xs )−1 σ(Xs ).
Step 2.

On the other hand, the identity
Z t
Z t
Yt = Y0 −h(Θr )dr +Zr dWr00

implies that Dt Yt = Zt , up to passing to a suitable version of these processes. It thus
follows that


Z 1
Z 1
∇Zr dWr (∇Xt )−1 σ(Xt ) .
Zt = Dg(X1 )∇X1 +Dh(Θr )∇Θr dr −tt

Settingη := ∇X exp

Z ·0

Dz f (Θr )dWr −

Z ·
0

 
1
|Dz f (Θr )|2 + Dy f (Θr ) dr ,2

we then deduce from Itô’s formula that
Z 1
Zt = E Dg(X1 )η1 +Dx h(Θr )ηr dr | Ft ηt−1 σ(Xt ) .
t

(2.9)

(2.10)

9

Discrete-time approximation of BSDEs

Step 3. The rest of the proof is based on a martingale argument. For sake of simplicity, let us consider the case Dh ≡ 0, ∇X ≡ 1 and σ is constant. Then
Zt = E [Dg(X1 ) | Ft ] σ ,and, for t ∈ [ti , ti+1 ],E |Zt − Zti |2 ≤ E |Zti+1 |2 − |Zti |2

so that, by 2. of Remark 2.2.2,
R2Z (π) ≤n−1
X Z ti+1
i=0

ti







E |Zt − Zti |2 dt ≤ |π|E |Z1 |2 = |π|E |Dg(X1 )|2 .

Also more technical, the general case is essentially treated similarly, by using the
Lipschitz continuity of the coefficients. The above proof is performed for Cb1 coefficients, but the final bound depends only on the Lipschitz constants. It can then be
extended by standard approximation techniques.
Theorem 2.3.2 ([74]). If the coefficients are Lipschitz, then R2Z (π) = O(|π|).
1

Combining the last theorem with Propositions 2.2.1 and 2.3.1 provides a |π| 2 convergence rate for the Backward Euler scheme.
1

Corollary 2.3.3. If the coefficients are Lipschitz, then Err(π) ≤ O(|π| 2 ).
Remark 2.3.4. Importantly, all these results are obtained under the single assumptionthat the coefficients are Lipschitz continuous. In particular, no ellipticity condition is
imposed on σ .
1Remark 2.3.5. The |π| 2 convergence rate is clearly optimal and is similar to the one
obtained for Forward SDEs.
Remark 2.3.6. Since no ellipticity assumption is made on σ , one component of X
can be chosen so as to coincide with time. It is therefore not a restriction to consider
a-priori time homogeneous dynamics. Note however, that a similar analysis could be
performed under the assumption that the b, σ and f are only 1/2-Hölder in the timecomponent, as opposed to Lipschitz continuous.
Remark 2.3.7. Similar results are obtained in [74] in the case where the terminal condition is of the form f (X) where f is a function of the whole path of X on [0, 1], undersuitable Lipschitz continuity assumptions on the space of continuous paths. The caseof non-regular terminal conditions has been discussed in [44]. The grid has then to
be refined near the terminal time t = 1 in order to compensate for the explosion of
the gradient Dv . However, their proof relies on PDE arguments and requires strongsmoothness assumptions on the coefficients as well as a uniform ellipticity conditionon σ , a condition that was not required so far.

10

B. Bouchard and R. Elie and N. Touzi

2.4 Weak expansion in the smooth case
So far, we have only discussed the strong convergence rate. But, as for linear problems,
which correspond to the case where h does not depend on Y and Z , we can also study
the weak approximation error. Such an analysis as been performed recently by Gobet
and Labart [42], thus extending previous results of Chevance [26] who considered thecase where h does not depend on Z .
Their result is of the form:
Yt − Ytπ

Zt − Z̄tπ= Dv(t, Xt )(Xt − Xtπ ) + O(|π|) + O |Xt − Xtπ |2′= D [Dv(t, Xt )σ] (Xt − Xtπ ) + O(|π|) + O |Xt − Xtπ |2 ,where (Xtπ )t≤1 is the continuous version of the Euler scheme of X and (Y π , Z π ) aredefined consistently with the explicit Backward Euler scheme byihYtπ = E Ytπi+1 + h Xtπi , Ytπi+1 , Z̄tπi (ti+1 − t) | Ftih1
E Ytπi+1 (Wti+1 − Wt ) | Ft , t ∈ [ti , ti+1 ) .
Z̄tπ =
ti+1 − t
In particular, this implies that v(0, X0 ) − Y0π = Y0 − Y0π = O(|π|). Moreover, when
the grid is uniform, this allows to deduce from the weak convergence of the process11|π|− 2 (X − X π ) that the process |π|− 2 (X − X π , Y − Y π , Z − Z̄ π ) weakly converges
too.

As usual, this requires much stronger assumptions than only our previous Lipschitz
continuity conditions, in particular they need to assume that σ is uniformly elliptic.

3

Extensions to other semilinear cases

The analysis of the previous Sections 2.2 and 2.3 consists in two parts:
1. Control the convergence rate of the discrete time approximation in terms of themodulus of regularity of Y and Z .
2. Study these modulus of regularity. Here, the difficult task consists in controlling
R2Z (π). In the above classical case, it is obtained by applying martingale-type techniques to a suitable representation of Z .
In the following, we show how these techniques can be adapted to more general situations corresponding to BSDEs with jumps, reflected BSDEs or BSDEs with random
time horizon.

3.1 BSDEs with jumps
The analysis of Chapter 2 can be extended in a very natural way to decoupled Forward Backward SDEs with jumps. Such equations appear naturally in optimal control

11

Discrete-time approximation of BSDEs

for jump diffusion processes, see [72], hedging problems, see [36], or portfolio management, see [10], in finance. From the PDE point of view, they have two important
applications.
Example 3.1.1 (Semilinear PDEs with integral term). Let us consider the solution(X, Y, Z, U ) of
Z t
Z t
Z tZ
Xt = X0 +b(Xr )dr +σ(Xr )dWr +β(Xr− , e)µ̄(de, dr) ,
0

Yt

=

g(X1 ) +0

Z 1
th (Θr ) dr −

0

Z 1
t

Zr′ dWr −

E

Z 1Z
t

Ur (e)µ̄(de, dr)

(3.1)

E

where µ̄(de, dt) = µ(de, dt) − λ(de)dt is a compensated compound point measure,
independent on W , with mark-space E ⊂ Rκ , and Θ := (X, Y, Z, Γ) withZ
Γ :=ρ(e)U (e)λ(de)E

for some bounded map ρ : E 7→ Rκ . Then, in the case where all the coefficients are
Lipschitz continuous, uniformly in the e-component, (Y, Z, U ) can be associated to the(unique viscosity) solution v of
−Lv − h (·, v, σ ′ Dv, I[v]) = 0, v(1−, ·) = g ,(3.2)

where L is the Dynkin operator associated to X andZ{v(t, x + β(x, e)) − v(t, x)} ρ(e) λ(de) ,I[v](t, x) :=Ethough the mapping, Yt = v(t, Xt ), Ut (e) = v(t, Xt− + β(Xt− , e)) − v(t, Xt− ) and
Zt = σ(Xt− )′ Dv(t, Xt− ), whenever it is well defined, see [5] and [18].
Example 3.1.2 (Systems of semilinear PDEs). More interestingly, it was shown by [65]and [71] that BSDEs of the form
Z t
Z t
Xt = X0 +b(Mr , Xr )dr +σ(Mr , Xr )dWr ,
Mt

= 1+

0

0

 Z t κ−1Xjµ({j}, dt) [modulo κ]

0 j=1

Yt

= g(M1 , X1 ) +

Z 1
t

−

Z 1Z
t

h (Mr , Xr , Yr , Zr , Ur ) dr −Ur (e)µ̄(de, dr)

Z 1

Zr dWr

t

(3.3)

E

allow to provide a probabilistic representation for systems ofPsemilinear PDEs. Herethe mark space is finite E = {1, . . . , κ − 1} and λ(de) ≡ λ κ−1
k=1 δk (e), where λ > 012

B. Bouchard and R. Elie and N. Touzi

is a constant. Namely, given a family (b̃m , σ̃m , h̃m , g̃m )m≤κ of Lipschitz continuousfunctions, if we set(b(m, x), σ(m, x), g(m, x)) = (b̃m (x), σ̃m (x), g̃m (x))

and
κ
 
X
 
γj ,
h(m, x, y, z, γ) = h̃m x, y + γ κ−m+1 , . . . , y + γ κ , y , y + γ 1 , . . . , z − λ
|{z}
j=1

m

then one can show that Yt = v Mt (t, Xt ) where v = (v 1 , . . . , v κ ) is the unique viscosity
solution of the system−Lm v m − hm (·, v, Dv m σm ) = 0, v m (1, ·) = gm (·) ,with Lm denoting the Dynkin operator associated to the Forward SDE with drift b̃mand volatility σ̃m .
Obviously, the formulations (3.1) and (3.3) can be combined so as to provide systems of semilinear PDEs with integral term.
The Backward Euler scheme for BSDEs with jumps of the form (3.1) was introduced
by Bouchard and Elie [18] and follows the ideas of Section 2.1. The formulation (3.3)
can be treated similarly, see Elie [35] for more details.
First, we define the Euler scheme X π of X as usual. The Euler scheme (Y π , Z π , Γπ )for (Y, Z, Γ) is then constructed as above by the Backward induction

i
h
1
π
π

=
Z̄
)
|
F
−
W
(W
E
Y

ti
ti
t
∆ti

 ti
i
h ti+1 R i+11ππΓ̄ti = ∆ti E Yti+1 E ρ(e)µ̄(de, (ti , ti+1 ]) | Fti

i
h



 Ytπ = E Ytπ | Fti + h Xtπ , Ytπ , Z̄tπ , Γ̄πt ∆ti
i

i+1

i

i

i

i

with terminal condition Y1π := g(X1π ).

Remark 3.1.3. In the case where E is a finite set {e1 , . . . , eκ } , one can always choose
ρ of the form ρi (ej ) = 1i=j /λ({ei }) so that Γ = (U (ei ))i≤κ . The above procedure thenallows to approximate U directly. It can obviously be extended to the approximationof any integrated version of e 7→ U (e) with respect to λ. We only concentrate here on
Γ as defined above because this is the term which appears in h and its approximationis therefore necessary to approximate Y .
As in the no-jump case, rather standard BSDE techniques allow to show that theapproximation error
#"
X Z ti+1|Zt − Z̄tπi |2 dt
Err(π)2 := max E
sup |Yt − Ytπi |2 + E 
i≤n−1



+E 

t∈[ti ,ti+1 ]X Z ti+1

i≤n−1

ti

i≤n−1



|Γt − Γ̄πti |2 dt ,tiDiscrete-time approximation of BSDEs

13

can be controlled in terms of |π|, R2Y (π), R2Z (π) and an additional term that takes
jumps into accountZ ti+1
X Z ti+11
Γs ds | Fti .
E
R2Γ (π) := E |Γt − Γ̄ti |2 dt with Γ̄ti :=∆titii≤n−1 ti

Under our Lipschitz continuity assumptions, the 1/2-Hölder continuity in time andLipschitz continuity in space of the deterministic map v is rather standard. In view of
the identities Yt = v(t, Xt ) and Ut (e) = v(t, Xt− + β(Xt− , e)) − v(t, Xt− ), see [18]for details, this readily implies that R2Y (π) + R2Γ (π) ≤ C|π| for some C > 0.
It thus remains to study the squared modulus of continuity R2Z (π) of Z . Following
the approach of Zhang [74], it was shown by Bouchard and Elie [18] that R2Z (π) can
be controlled by |π| whenever one of the following assumptions holds:
Ha) For all e ∈ E , the map x ∈ Rd 7→ β(x, e) admits a Jacobian matrix Dβ(x, e) suchthat(x, ξ) ∈ Rd × Rd 7→ a(x, ξ; e) := ξ ′ (Dβ(x, e) + Id )ξ

satisfies, uniformly in (x, ξ) ∈ Rd × Rda(x, ξ; e) ≥ |ξ|2 K −1ora(x, ξ; e) ≤ −|ξ|2 K −1 ,for some K > 0.
Hb) b, σ , β(·, e), h and g are Cb1 with Lipschitz first derivatives, uniformly in e ∈ E .

Remark 3.1.4. The condition Ha) ensures that the tangent process ∇X is invertible,
with inverse satisfying suitable integrability conditions, which is necessary to reproduce the argument presented in Section 2.3, recall the definition of the term η . Under
Hb), it can be shown, see [35], that Dv is 1/2-Hölder-continuous in time and Lipschitzcontinuous in space. This allows to show that R2Z (π) can be controlled by |π| by simply
using the 1/2-Hölder regularity of X in L2 .
1

This allows to provide a |π| 2 convergence rate.
Theorem 3.1.5 ([18]). Let Ha) or Hb) hold, assume that λ is bounded and that the1
coefficients are Lipschitz, uniformly in the e-variable. Then, Err(π) ≤ O(|π| 2 ).
Remark 3.1.6. So far, it seems that only the case of jump measures with finite activity has been treated. In [18], the compensatory measure is assumed to be bounded,
but it can be generalised without much difficulties under suitable integrability conditions. Moreover, pure jump Lévy processes can be approximated by compound point
measures by truncating the small jumps. It would therefore be natural to approximate a BSDE driven by a Lévy process with infinite activity by a sequence of BSDEs
driven by compound point measures to which the above analysis would apply. This
would introduce an additional error related to the approximation of the original random measure. We refer to [70] and the references therein for the approximation ofLévy processes.

14

B. Bouchard and R. Elie and N. Touzi

3.2 Reflected BSDEs
Let us now turn to the approximation of the solution (Y, Z, A) of a reflected BSDE ofthe form
Z 1
Z 1
Yt = g(X1 ) +h (Xs , Ys , Zs ) ds −Zs′ dWs + A1 − At
(3.4)
t

Yt

t

≥ f (Xt ) , 0 ≤ t ≤ 1 ,where A is a cadlag adapted non-decreasing process satisfying
Z 1(Yt − f (Xt )) dAt = 0 .
0In the case where h = 0 and f = g , Y is the Snell envelope of g(X), which, infinance, corresponds to the super-hedging price of the American option with payoff g .
More generally, the solution (Y, Z, A) of (3.4) is related to semilinear PDEs with freeboundaries of the formmin {−Lv − h(·, v, σ ′ Dv) , v − f } = 0 , v(1−, ·) = g(3.5)

see [47].
For such BSDEs, it is natural to modify the Backward Euler scheme (2.4) as follows
so as to take into account the reflection:ihZ̄tπi = (ti+1 − ti )−1 E Ytπi+1 ∆Wti | FtiihỸtπi = E Ytπi+1 | Fti + h(Xtπi , Ỹtπi , Z̄tπi )∆ti


(3.6)
Ytπi = R ti , Xtπi , Ỹtπi , i ≤ n − 1 ,with the terminal condition Y1π = g(X1π ) and whereR(t, x, y) :=y + [f (x) − y]+ 1{t∈ℜ\{0,T }} , (t, x, y) ∈ [0, 1] × Rd+1

for a partition ℜ := {0 = r0 < r1 < . . . < rκ = 1} ⊃ π .
Numerical procedures for such Backward Euler schemes have been first studied in
Bally and Pagès [3] in the case where π = ℜ, h does not depend on Z and g = f . The
fact that h is independent on Z allows them to use the formulation of Y in terms of
optimal stopping problem to study its discrete time approximation error. They retrieve1
a |π|− 2 control for the convergence rate when the coefficients are Lipschitz continuous
and improve it to |π| when g is semiconvex and X can be perfectly simulated on thetime-grid. The error on Z is not discussed.
A more general analysis has then been performed by Ma and Zhang [60], who ob1
tained a bound for the convergence rate of order |π|− 4 when f depends on Z whenever
it is Cb2 and σ is uniformly elliptic. Note that the arguments leading to (2.10) cannot

Discrete-time approximation of BSDEs

15

be used here, since this would involve taking the Malliavin derivative of A, which has
no reason to be well-defined. Instead, they rely on a representation of Z of the form
(d = 1 to avoid notational complexities):


Z 1
Z 1
t
t
t
(3.7)
Nr dAr | Ft σ(Xt )h (Xr , Yr , Zr ) Nr dr +Zt = E g(X1 )N1 +
t

t

where
Nst :=

1
s−t

Z s
t

σ(Xr )−1 ∇Xr dWr (∇Xt )−1 .

Their argument can be formally explained in three steps. First, they derive a representation of the Z -component associated to a discretely reflected BSDE of the form
Z 1
Z 1
Ytb = g(X1 ) +h(Θbr )dr −(Zrb )′ dWr + Ab1 − Abt
(3.8)
t

t

i+
h
P
b
1rj ≤t . The first represenwhere Θb := (X, Y b , Z b ) and Abt := κ−1
j=1 h(Xrj ) − Yrjtation is obtained by considering the Malliavin derivatives of (Y, Z) as in Section 2.3,
see Step 1, on each interval [rj , rj+1 )Z rj+1
Z rj+1Dt Zrb dWr , t ∈ [rj , rj+1 )Dt h(Θbr )dr −Ztb = Dt Ytb = Dt Yrbj+1 − +
t

t

i+
h
where Dt Yrbj+1 − = Dt Yrbj+1 + Dt h(Xrj+1 ) − Yrbj+1 since, by construction, Yrbj+1 =hi+
Yrbj+1 + h(Xrj+1 ) − Yrbj+1 . This implies that, for t < rj+1 ,Ztb = Dt g(X1 ) +Z 1t

Dt h(Θbr )dr −Z 1
t

Dt Zrb dWr +

κ
X

k=j+1


+
Dt h(Xrk ) − Yrbk.

The fact that the discretely reflected BSDE is differentiable in the Malliavin sense canbe easily proved by using a backward induction argument and the fact that z 7→ z + isa Lipschitz function. Second, they perform an integration by parts, in the Malliavin
sense, in order to get rid of the Malliavin derivatives, and obtain a formulation of the
form
#
"
Z 1
κ
X
b
t
b
t
t
b
Zt = E g(X1 )N1 +h(Θr )Nr dr +Nrk ∆Ark 1t<rk | Ft σ(Xt )
t

k=0


+
for t ∈ [rj , rj+1 ), where ∆Abrk := h(Xrk ) − Yrbk . Finally, they pass to the limit on
the grid of reflection times so as to recover (3.7) for the continuously reflected BSDE.
This formulation is in itself very interesting and can be viewed as a generalisationof the representation obtained in [38] for the linear non reflected case. However, it
requires a uniform ellipticity condition on σ . Moreover, the exploding behaviour of Nrt
as r → t is a major drawback in the regularity analysis of Z .

16

B. Bouchard and R. Elie and N. Touzi

Another regularity analysis was carried out in Bouchard and Chassagneux [16]. The
main difference is that they do not perform the integration by parts but rather provide
a representation of Z b in terms of the next reflection time:

b
(3.9)
Zt = E Dg(X1 )η1b 1τj =1 + Df (Xτj )ητbj 1τj <1
+

Z τj
t


Dx h(Θbr )ηrb dr | Ft (ηtb )−1 σ(Xt ) ,for t ∈ [rj , rj+1 ), where τj := inf{r ∈ ℜ | r ≥ rj+1 , f (Xr ) > Yrb } ∧ 1 and η b isdefined as in (2.9) above with (X, Y b , Z b ) in place of (X, Y, Z). This is the analogueof (2.10) up to the stopping time τj .
This formulation does not require any uniform ellipticity assumption on σ and turnout to be more flexible.
Let us expose the argument in the simpler case where d = 1, g ≡ f , h ≡ 0, ∇X ≡ Id
and σ is constant. Set
Vtj := E Df (Xτj ) | Ft σ so that Ztb = Vtj if t ∈ [rj , rj+1 ) ,and note that the martingale property of V j implies thatiihh
≤ E |Vtji+1 |2 − |Vtji |2 for t ∈ [ti , ti+1 ] .
E |Vtj − Vtji |2

Letting ij be defined by rj = tij , recall that π ⊂ ℜ, it follows that
−1

|π|

n−1
X
i=0

E

Z ti+1ti|Ztb − Ztbi |2 dt≤ Σ :=

κ−1
X−1
X ij+1
j=0 k=ijih
E |Vtjk+1 |2 − |Vtjk |2

where, by direct computations,
Σ

=

κ−1
X
j=0

κ−1
i
h
i
X h
2
j 2
|
|
−
|V
E |Vrj−1
E |Vrjj+1 |2 − |Vrjj |2 ≤ C +
rj
j
j=1

and C > 0 denotes a generic constant independent on π and ℜ. On the other hand, it
follows from the Cauchy–Schwarz inequality that, if Df is bounded,| ≤ C E β | Frj |Vrjj − Vrj−1| |Vrjj + Vrj−1|2 − |Vrjj |2 ≤ |Vrjj − Vrj−1|Vrj−1|jjjjwhere β := 1 + supt≤1 |Xt |. Moreover, if Df ∈ Cb2 , then
= E (Df (Xτj ) − Df (Xτj−1 ) | Frj σ ≤ C E (τj − τj−1 ) | Frj σ .
Vrjj − Vrj−1jCombining the above inequalities, this implies



κ−1
n−1
X  
X Z ti+1E E β | Frj E τj − τj−1 | Frj 
E|Ztb − Ztbi |2 dt ≤ C|π| 1 +
i=0

ti

j=1

17

Discrete-time approximation of BSDEs

which is of order of |π| since the expectation on the right-hand side simply equalsE [β(τj − τj−1 )].
The general case is treated along the same main argument: play with the differenceof the next reflection times τj − τj−1 . It allows to control the modulus of regularity
R2Z b (π) of Z b , defined as R2Z (π) for Z b is place of Z , as follows:
Theorem 3.2.1 ([16]). Assume that all the coefficients are Lipschitz continuous. Then,
R2Z b (π) ≤n−1
X Z ti+1
ti

i=0



E |Ztb − Ztbi |2 dt

≤

C |π|




1
1 + β|π|− 2 + α(ℜ) ,1

where α(ℜ) = |ℜ|− 2 and β = 1 under the assumption
H1 ) : f is Cb1 with Lipschitz first derivative
and, α(ℜ) = 1 and β = 0 under the assumption
H2 ) : f is Cb2 with Lipschitz second derivative and σ is Cb1 with Lipschitz first derivative.
This regularity property of Z b then allows to provide a convergence rate for the Euler
scheme (Y π , Z̄ π ) to (Y b , Z b ).
Theorem 3.2.2 ([16]). Assume that the coefficients are Lipschitz. Then,

max E
i<n

"

sup

t∈(ti ,ti+1 ]

n−1
X
i=0

E

|Ytπi+1 − Ytb |2

Z ti+1ti# 12|Z̄tπi − Ztb |2 dt



≤



1
1
O αY (ℜ) |π| 2 + β |π| 4

≤



1
O αZ (ℜ) |π| + β |π| 21with (αY (ℜ), αZ (ℜ)) = (|ℜ|− 4 , |ℜ|−1 ) and β = 1 under H1 ), and, (αY (ℜ), αZ (ℜ)) =(1, |ℜ|−1 ) and β = 0 under H2 ).
1

They also show that the above result holds true with αZ (κ) = |ℜ|− 2 under H1 ), andαZ (κ) = 1 under H2 ), when (Xtπi )i≤n is replaced by (Xti )i≤n in the Backward Euler
scheme (3.6).
To obtain a convergence rate of the approximation error of Y , it then essentially suf1
fices to consider the case ℜ = π and show that Y b approximates Y at least at a rate |ℜ| 2 .
The approximation of the Z -component is more involved and requires the introductionof another approximation scheme.
Theorem 3.2.3 ([16]). Assume that the coefficients are Lipschitz. Then, there exists

18

B. Bouchard and R. Elie and N. Touzi

C > 0, independent on π , such thatEmax |Ytπi − Yti |2
i≤n

 12

+ max Ei<nR2Z (π) ≤1"|Ytπi+1 − Yt |2

supt∈(ti ,ti+1 ]

n−1
X
i=0

E

Z ti+1ti# 21= O(α(π))

1|Z̄tπi − Zt |2 dt = O(|π| 2 ) ,1

where α(π) = |π| 4 under H1 ) and α(π) = |π| 2 under H2 ).
When H2 ) holds and (Xtπi )i≤n is replaced by (Xti )i≤n in the Backward Euler scheme(3.6), the (squared) error on Z and the (squared) modulus of regularity R2Z (π) are1
shown to be controlled by |π|, which corresponds to a convergence rate of order |π| 2
as in the non-reflected case.
Remark 3.2.4. The representation (3.9) of Z b is of own interest as it provides a naturalMonte-Carlo estimator for the so-called “delta” of Bermudean options. Such a result
was already known in finance for American options, under a uniform ellipticity condition, see [41].
Remark 3.2.5. Extensions to doubly reflected and multivariate reflected BSDEs havebeen considered in Chassagneux [23] and [24].

3.3 BSDEs with random time horizon
By BSDE with random time horizon, we mean the solution (Y, Z) of
Yt= g(Xτ ) +Z 1
t1s<τ h(Xs , Ys , Zs )ds −Z 1
t

Zs′ dWs , t ≤ 1 ,(3.10)

where τ is the first exit time of (t, Xt )t≤T from a cylindrical domain D = [0, 1) × O.
The first step toward the definition of a Backward Euler scheme, consists in choosing
a suitable approximation of the random time horizon τ . In Bouchard and Menozzi [19],the authors consider the first exit time (on the grid) τ̄ of the Euler scheme X πτ̄ := inf{t ∈ π : (t, X̄t ) ∈/ D} .

They then reproduce the Backward scheme (2.4)

i
h
1
π
 Z̄tπ =
|
F
∆W
E
Y
t
t
i
i
t
∆thi
i+1
i
i

 Ytπ = E Ytπ | Fti + h Xtπ , Ytπ , Z̄tπ ∆ti
i

i

i+1

with the terminal condition

Y1π = g(Xτ̄π ) .

i

i

19

Discrete-time approximation of BSDEs

Following the approach of [20] and [74], they show that the discrete-time approximation error can be controlled as:



Z τ
2
2
2
2
Err(π) ≤ C |π| + RY (π) + RZ (π) + E ξ|τ − τ̄ | + 1τ̄ <τ
|Zs | ds
(3.11)
τ̄where C > 0 is independent of π and ξ is a positive random variable satisfying E [ξ p ] ≤Cp′ for all p ≥ 2, for some Cp′ independent of π .
Not surprisingly, the error now depends on an additional term which is related tothe approximation error of τ by τ̄ which has been widely studied under the assumptionthat the coefficients are smooth enough and σ is uniformly elliptic. In particular, Gobet[39], and, Gobet and Menozzi [45] have derived an expansion of the form E [τ − τ̄ ] =1
1
C|π| 2 + o(|π| 2 ).
A similar result can be obtained without uniform ellipticity, whenever O is a C 2
domain with a compact boundary satisfying a uniform non-characteristic condition.
Theorem 3.3.1 ([19]). Assume that the coefficients of X are Lipschitz continuous andthat O is a C 2 domain with a compact boundary satisfying a uniform non-characteristic
condition. Then, for each ε ∈ (0, 1) and each positive random variable ξ satisfying
E [ξ p ] < ∞ for all p ≥ 1,
h 

2 i
≤ Oεξ |π|1−ε ,E E ξ |τ − τ̄ | | Fτ+ ∧τ̄

where τ+ := inf{t ∈ π : τ ≤ t}.

This shows in particular that, for each ε ∈ (0, 1/2),
E [|τ − τ̄ |] ≤1

Oε |π| 2 −εwhich, up to the ε term, is consistent with the result of Gobet and Menozzi [45]. Moreover, it is easily extended to the case where
D1) O is an intersection of C 2 domains with compact boundary, and it satisfies a uniform non-characteristic condition outside the corners and a uniform ellipticity
condition on a neighbourhood of the corners.

The next step toward a convergence rate consiststhe modulus of reg in studyingRτularity R2Y (π) and R2Z (π), and the last term E 1τ̄ <τ τ̄ |Zs |2 ds .

In contrast with the case where O = Rd , the uniform Lipschitz continuity of the(viscosity) solution v of the associated semilinear Cauchy Dirichlet problem
−Lv − h(·, v, σ ′ Dv) = 0 on D , v = g on ∂p D := ([0, 1) × ∂O) ∪ {1} × Ō ,(3.12)
is by no means obvious, and usually requires additional assumptions than thesimple Lipschitz continuity of the coefficients. It is shown in [19], by adapting
standard barrier techniques, when g is smooth enough:20

B. Bouchard and R. Elie and N. Touzi

Hg) g ∈ Cb1,2 ([0, 1] × Rd ),and the domain satisfies a uniform exterior sphere condition as well as a uniform
truncated interior cone condition:
D2) For all x ∈ ∂O, there is y(x) ∈ Oc , r(x) ∈ [L−1 , L], L > 0, and δ(x) ∈ B(0, 1)
such that B̄(y(x), r(x)) ∩ Ō = {x} and {x′ ∈ B(x, L−1 ) : hx′ − x, δ(x)i ≥(1 − L−1 )|x′ − x|} ⊂ Ō .
Here B(z, r) denotes the open ball of centre z and radius r, B̄(z, r) is its closure.
Theorem 3.3.2 ([19]). Assume that all the coefficients are Lipschitz continuous andthat the conditions Hg), D1) and D2) hold. Then, v is uniformly 1/2-Hölder continuous in time and Lipschitz continuous in space.
In particular,
this allows to show that R2Y (π) is controlled by |π| and the last termRτ1
E 1τ̄ <τ τ̄ |Zs |2 ds is controlled by E [|τ − τ̄ |] ≤ C ε |π| 2 −ε , for ε > 0.

It thus remains to study the modulus of continuity of Z . As in the reflected case, it
is not possible to apply the techniques leading to (2.10) since it would involve taking
the Malliavin derivative of g(Xτ ), which has no reason
R ·to be well-defined. Instead, we
can remark that, if v is smooth enough, Dv(·, X)η + 0 Dx h (Θs ) ηs ds is a martingale
on [0, τ ], recall the definition of η in (2.9). This implies that (d = 1 for notational
simplicity)Z τ
Zt = Dv(t, Xt )σ(Xt ) = E Dv(τ, Xτ )ητ +Dx h (Θs ) ηs ds | Ft ηt−1 σ(Xt ) .
tThis formula is very close to (2.10) except that Dg is replaced by Dv . But, since Dv
is bounded, one can actually repeat the martingale-type argument of Step 3 of Section
2.3. The smoothness of v can be obtained for smooth coefficients and domains, undera uniform ellipticity condition of σ , which can then be relaxed by standard approximation arguments.
Theorem 3.3.3 ([19]). Assume that all the coefficients are Lipschitz continuous andthat the conditions Hg), D1) and D2) hold. Then,
R2Y (π) + R2Z (π)≤

O(|π|) .

1

This provides the following |π| 4 −ε convergence rate.
Corollary 3.3.4 ([19]). Assume that all the coefficients are Lipschitz continuous andthat the conditions Hg), D1) and D2) hold. Then, for all ε ∈ (0, 1), 1 
Err(π) ≤ Oε |π| 2 −ε .
As a matter of facts, the global error is driven by the approximation error of the exit
time which propagates backward thanks to the Lipschitz continuity of v .

Remark 3.3.5. It is shown in [19] that one can actually achieve a convergence rate of1order |π| 2 −ε if one computes the approximation error only up to τ̄ ∧ τ .

21Discrete-time approximation of BSDEs

4

Extension to fully nonlinear PDEs

The discretisation of backward stochastic differential equations can be viewed as a probabilistic approach for the numerical resolution of semilinear partial differential equations in parabolic or elliptic form, and possibly with a nonlocal integral term. This restriction on the nonlinearity of the corresponding PDEs excludes the important generalclass of Hamilton–Jacobi–Bellman equations which appear naturally as the infinitesimal version of the dynamic programming principle for stochastic control problems.
An extension of the theory of backward SDEs to the second order case was initiatedby [25] and is still in progress. An alternative approach for such an extension in thecontext of a specific form of the nonlinearity was also introduced by Peng in various
papers, see e.g. [69], and is based on the new concept of G-Brownian motion.
We shall not elaborate more on this aspect. Instead, our objective is to concentrateon the probabilistic numerical implications as discussed in Section 5 of [25]. Following [37], we provide a natural presentation of the discrete-time approximation without
appealing to the theory of backward stochastic differential equations. This presentation shows in particular the important connection between the discretisation schemesof this paper and the finite differences approximations of solutions of PDEs. Finally,
we were not able to extend the methods of proofs outlined in the previous sections.
The convergence results and the rates of convergence are derived by means of the PDE
approach for monotonic schemes.

4.1 Discretisation
Consider the fully nonlinear parabolic PDE:

−Lv − F ·, v, Dv, D2 v = 0, on [0, 1) × Rd ,v(1, ·) =dg, on ∈ R ,

(4.1)
(4.2)

where, given bounded and continuous maps b and σ from R+ × Rd to Md and Rd , thelinear second order operator L is defined by1 ∂v+ hb, Dvi + Tr σσ ′ D2 v ,
Lv :=
∂t
2
and the nonlinearity is isolated in the map F : (t, x, r, p, γ) ∈ R+ ×Rd ×R×Rd ×Sd 7−→
F (t, x, r, p, γ) ∈ R.
For a positive integer n, set h := 1/n, ti := ih, i = 1, . . . , n, and define:X̂ht,x := x + b(t, x)h + σ(t, x)(Wt+h − Wt ),(4.3)

which is the one-step Euler discretisation of the diffusion corresponding to the operator
L.
Assuming that the PDE (4.1) has a classical solution, it follows from Itô’s formulathatZ ti+1Eti ,x v ti+1 , Xti+1
Lv(t, Xt )dt= v (ti , x) + Eti ,xti

22

B. Bouchard and R. Elie and N. Touzi

where we ignored the difficulties related to local martingale part, and Eti ,x :=E[·|Xti = x] denotes the expectation operator conditional on {Xti = x}. Since v
solves the PDE (4.1), this provides

Eti ,x v ti+1 , Xti+1=

v(ti , x) − Eti ,xZ ti+12F (·, v, Dv, D v)(t, Xt )dt .

ti

By approximating the Riemann integral, and replacing the process X by its Euler discretisation, this suggests the following approximation of the value function vv h (1, .) = g

and v h (ti , x) = Th [v h ](ti , x)

(4.4)

where, for a function ψ : R+ × Rd −→ R with exponential growth, we denote:hiTh [ψ](ti , x) := E ψ(ti+1 , X̂hti ,x ) + hF (·, Dh ψ) (ti , x),Dh ψ := Dh0 ψ, Dh1 ψ, Dh2 ψ ,and for k = 0, 1, 2, it follows from an easy integration by parts argument that:
Dhk ψ(ti , x):= E[Dk ψ(ti+1 , X̂hti ,x )] = E[ψ(X̂hti ,x )Hkh (ti , x)],

(4.5)

′
Wh
−1 Wh Wh − h1d
σ −1 .
, H2h = (σ ′ )
h
h2

(4.6)

where
H0h = 1, H1h = (σ ′ )−1

Remark 4.1.1. In the semilinear case, i.e. F is independent of the argument γ , thenumerical scheme coincides with the discretisation of BSDEs as studied in the previous
sections of this paper.
Remark 4.1.2. Other choices can be made for the above integration by parts. For
instance, we also have:
"
#
′
W
W
h/2
h/2
t,x
D2h ϕ(t, x) = E ϕ(X̂h )(σ ′ )−1
σ −1 ,(h/2) (h/2)
which shows that the backward scheme (4.4) is very similar to the probabilistic numerical algorithm suggested in [25].
Remark 4.1.3. So far, the choice of the drift and the diffusion coefficients b and σ inthe nonlinear PDE (4.1) is arbitrary, and was only used to fix the underlying diffusion
X . Our convergence result will however place some restrictions on this choice.
Remark 4.1.4. The integration by parts in (4.5) is intimately related to finite differences for numerical schemes on deterministic grids. We formally justify this connection for d = 1, b ≡ 0, and σ = 1 for simplicity:

23

Discrete-time approximation of BSDEs•


Let {wj , j ≥ 1} be independent r.v. distributed as 21 δ√h + δ−√h . Then, thePbinomial random walk approximation of the Brownian motion Ŵtk := kj=1 wj ,tk := kh, k ≥ 1, suggests the following approximation:√√ψ(t, x + h) − ψ(t, x − h)
t,x
1
h
√
≈
Dh ψ(t, x) := E ψ(t + h, Xh )H1,
2 hwhich is the centred finite differences approximation of the gradient.
√√• Let {w , j ≥ 1} be independent r.v. distributed as 1j
6 δ{ 3h} + 4δ{0} + δ{− 3h} .

Then, the trinomial random walk approximation of the Brownian motion Ŵtk :=Pk
j=1 wj , tk := kh, k ≥ 1, suggests the following approximation:

Dh2 ψ(t, x) := E ψ(t, x +ψ(t + h, Xht,x )H2h ≈√3h)−2ψ(t, x)+ψ(t, x −
3h

√

3h)

,

which is the centred finite differences approximation of the Hessian.
In view of the above interpretation, the numerical scheme studied in this paper can beviewed as a mixed Monte Carlo–Finite Differences algorithm. The Monte Carlo component of the scheme consists in the choice of an underlying diffusion process X . The
finite differences component of the scheme consists in approximating the remaining
nonlinearity by means of the integration-by-parts (4.5).

4.2 Convergence of the discretisation
Definition 4.2.1. We say that (4.1) has strong comparison of bounded solutions if forany bounded upper semicontinuous viscosity supersolution v and any bounded lower
semicontinuous subsolution v on [0, 1) × Rd , satisfying••

either v(1, ·) ≥ v(1, ·),or the viscosity property of v and v holds true on [0, 1] × Rd ,we have v ≥ v .
The strong comparison principle is an important notion in the theory of viscosity
solutions which allows to handle situations where the boundary condition g is not compatible with the equation.
In the sequel, we denote by Fr , Fp and Fγ the partial gradients of F with respect tor, p and γ , respectively. We recall that any Lipschitz function is differentiable a.e.
Assumption F. The nonlinearity F is Lipschitz-continuous with respect to (r, p, γ)
uniformly in (t, x) and sup(t,x)∈[0,T ]×Rd |F (t, x, 0, 0, 0)| < ∞. Moreover, F is uniformly elliptic and dominated by the diffusion of the linear operator L, i.e.
εId ≤ Fγ ≤ a on Rd × R × Rd × Sd

for some ε > 0.

(4.7)

24

B. Bouchard and R. Elie and N. Touzi

Theorem 4.2.2. Let Assumption F hold true, and assume that the fully nonlinear PDE(4.1) has strong comparison of bounded solutions. Then for every bounded function g ,there exists a bounded function v so thatv h −→ v

locally uniformly.

In addition, v is the unique bounded viscosity solution of the relaxed boundary problem(4.1)–(4.2).
Remark 4.2.3. The restriction to bounded terminal data g in the above Theorem 4.2.2
can be relaxed by an immediate change of variable, thanks to the boundedness conditions on the coefficients b and σ , see [37].
Remark 4.2.4. Theorem 4.2.2 states that the right hand-side inequality of (4.7) (i.e.
diffusion must dominate the nonlinearity in γ ) is sufficient for the convergence of theMonte Carlo–Finite Differences scheme. We do not know whether this condition isnecessary. However, in [37], it is shown that this condition is not sharp in the simple
linear case, while the numerical experiments reveal that the method may have a poor
performance in the absence of this condition.
The rest of this section is dedicated to the proof of Theorem 4.2.2 which is basedon the convergence of monotonic schemes of Barles and Souganidis [9]. Since theprobabilistic community is not familiar with this result we report a complete proof.
Properties of the discretisation scheme
properties of the scheme.

The approach of [9] is based on three main

1 The first ingredient is the so-called consistency property, which states that for every
smooth function ϕ with bounded derivatives and (t, x) ∈ [0, T ] × Rd :lim(t′ , x′ ) → (t, x)(h, c) → (0, 0)
t′ + h ≤ 1

(c + ϕ − Th [c + ϕ]) (t′ , x′ )= − LX ϕ + F (·, ϕ, Dϕ, D2 ϕ) (t, x).
h(4.8)
The verification of this property follows from an immediate application of Itô’s formula.
2 The second ingredient is the so-called stability property:
the family (v h )h is L∞ − bounded, uniformly in h.

(4.9)

This property is implied by the boundedness of g and F (t, x, 0, 0, 0) (Assumption F),
together with the stronger property that for every L∞ -bounded functions ϕ, ψ : [0, 1] ×
Rd −→ R:|Th [ϕ] − Th [ψ]|∞ ≤ |ϕ − ψ|∞ (1 + Ch) for some C > 0,

(4.10)

25

Discrete-time approximation of BSDEs

For ease of presentation, we report the proof of (4.10) in the one-dimensional case; thegeneral multi-dimensional case follows the same line of argument. Set f := ϕ − ψ . By
the mean value theorem there exists some θ = (t, x, r̄, p̄, γ̄) that:
W2 − hWh
Fp (θ) + h 2
Fγ (θ)(Th [ϕ] − Th [ψ]) (t, x) = E f (X̂ht,x ) 1 + hFr (θ) +
σ
σ h



2
Fγ (θ)
Fp (θ)
t,x
2= E f (X̂h ) 1 −+ hFr − h
+ Ah (4.11)σ24Fγ (θ)
where
Ah

:=

Wh
√
σ h

√
q
Fp (θ) h,Fγ (θ) + p
2 Fγ (θ)(4.12)

Since 1 − σ −2 Fγ ≥ 0 and |Fγ−1 |∞ < ∞ by (4.7) of Assumption F, it follows from the
Lipschitz property of F that
|Th [ϕ] − Th [ψ]|∞ ≤ |f |∞ 1 − σ −2 Fγ (θ) + E[|Ah |2 ] + Ch
= |f |∞ 1 + h(4Fγ )−1 (θ)(Fp )2 (θ) + Ch ≤ |f |∞ (1 + Ch).
3 The final ingredient – and the most important – is the monotonicity. In the context
of our scheme (4.4), we have for all functions ϕ, ψ : [0, 1]×Rd −→ R with exponential
growth:ϕ≤ψ=⇒ Th [ϕ](t, x) ≤ Th [ψ](t, x) + Ch E[(ψ − ϕ)(t + h, X̂ht,x )] (4.13)

for some constant C > 0. To prove (4.13), we proceed as in the proof of the stability
result to arrive at (4.11). Since f := ψ − ϕ ≥ 0 in the present context, and Fγ ≤ σ 2 by
(4.7) of Assumption F, we deduce that:h Fp (θ)2t,x,(4.14)(Th [ψ] − Th [ϕ]) (t, x) ≥ E f (X̂h ) hFr (θ) −4 Fγ (θ)
and the required result follows from the Lipschitz property of F and the fact that|Fγ−1 |∞ < ∞ by (4.7).
The Barles–Souganidis monotone scheme convergence argument Given the above
consistency, stability and monotonicity properties, the convergence of the family (v h )htowards some function v which is the unique viscosity solution of the PDE (4.1) follows from [9]. We report the full argument for completeness.
From the stability property the semi-relaxed limitsV∗ (t, x) :=lim inf(t′ ,x′ ,h)→(t,x,0)v h (t′ , x′ )and V ∗ (t, x) :=lim supv h (t′ , x′ )(t′ ,x′ ,h)→(t,x,0)

are finite lower-semicontinuous and upper-semicontinuous functions, respectively. Also,
they obviously inherit the boundedness of the family (v h )h . The key-point is to show

26

B. Bouchard and R. Elie and N. Touzi

that V∗ and V ∗ are respectively viscosity supersolution and subsolution of the PDE(4.1)–(4.2), which implies by the strong comparison assumption that V∗ ≤ V ∗ . Since
the converse inequality is trivial, this shows that V∗ = V ∗ is the unique bounded solution of (4.1), thus completing the proof of convergence.
In the rest of this section, we show that V∗ is a viscosity supersolution of the PDE(4.1). A symmetric argument applies to prove the corresponding property for V ∗ . Let(t0 , x0 ) ∈ [0, 1) × Rd and ϕ ∈ C 2 [0, 1] × Rd be such that
0 = (V∗ − ϕ)(t0 , x0 )= min(V∗ − ϕ)Without loss of generality, we may assume that (t0 , x0 ) is a strict minimiser of (V∗ −ϕ).
Let (hn , tn , xn ) be a sequence in (0, 1] × B1 (t0 , x0 ) with(hn , tn , xn ) −→ (0, t0 , x0 )and v hn (tn , xn ) −→ V∗ (t0 , x0 ),and define (t′n , x′n ) by
cn := (v∗hn − ϕ)(t′n , x′n ) =min(v hn − ϕ),B0

where B0 ⊂ [0, 1] × Rd is a closed ball containing (t0 , x0 ) in its interior. By definition,
we have cn + ϕ(t′n , x′n ) = v hn (t′n , x′n ) = Thn [v hn ](t′n , x′n ). Since v hn ≥ cn + ϕ, itfollows from the monotonicity property (4.13) that:cn + ϕ(t′n , x′n ) ≥ Thn [cn + v hn ](t′n , x′n ).

(4.15)

Since the sequence (t′n , x′n )n is bounded, it converges to some (t1 , x1 ) after possibly
passing to a subsequence. Observe that(V∗ −ϕ)(t0 , x0 ) = lim (v hn −ϕ)(tn , xn ) ≥ lim inf (v hn −ϕ)(t′n , x′n ) ≥ (V∗ −ϕ)(t1 , x1 ).
n→∞

n→∞

Since (t0 , x0 ) is a strict minimiser of the difference V∗ − ϕ, this shows that (t1 , x1 ) =(t0 , x0 ), and cn = (v∗hn − ϕ)(t′n , x′n ) −→ 0. We now go back to (4.15), normalise byhn and send n to infinity. By the consistency property (4.8), this provides the required
result:−Lϕ(t0 , x0 ) − F ., ϕ, Dϕ, D2 ϕ (t0 , x0 ) ≥ 0.

Remark 4.2.5. In the context of nonlinear PDEs of the HJB type, see Assumption HJB
below, Bonnans and Zidani [14] introduced finite differences approximations whichobey to the monotonicity condition. The choice of the discretisation turns out to bepretty involved and is specific to each choice of control. Our probabilistic scheme has
the nice property to be automatically monotonic.

4.3 Rate of convergence of the discretisation of HJB equations
We next provide an error estimate for our probabilistic numerical scheme. For ease ofpresentation, we assume that the nonlinearity F satisfies the condition
1Fr − Fp′ Fγ−1 Fp
4≥ 0,(4.16)

Discrete-time approximation of BSDEs

27

which implies that the monotonicity property (4.13) is strengthened to:ϕ≤ψ=⇒ Th [ϕ] ≤ Th [ψ],(4.17)see (4.14). In [37], it is shown that our main results hold without Condition (4.16).
The subsequent argument for the derivation of the error estimate is crucially basedon the following comparison principle satisfied by the scheme.
Lemma 4.3.1. Let Assumption F holds true, and consider two arbitrary bounded functions ϕ and ψ satisfying:h−1 (ϕ − Th [ϕ]) ≤ g1

and h−1 (ψ − Th [ψ]) ≥ g2

(4.18)

for some bounded functions g1 and g2 . Then, for every i = 0, . . . , n:(ϕ − ψ)(ti , x) ≤ eβ(T −ti ) |(ϕ − ψ)+ (1, ·)|∞ + (1 − h)eβ(1−ti ) |(g1 − g2 )+ |∞ (4.19)

for some parameter β > |Fr |∞ .
Before stating precise results, let us explain the key idea as introduced by Krylov[52, 53, 54]. We also refer to the lecture notes by Bonnans [13] which provide a clear
summary in the context of infinite horizon stochastic control problems.
Key-idea for the lower bound Given the solution v of the nonlinear PDE (4.1),
suppose that there is a function uε satisfyinguε

is a classical subsolution of (4.1) and v − C(ε) ≤ uε ≤ v(4.20)

for some function C(ε). Let
Rh [uε ] :=uε − Th [uε ]+ Luε + F ., uε , Duε , D2 uε ,hand let R(h, ε) be a bound on Rh [uε ]:|Rh [uε ](t, x)| ≤ R(h, ε) for every(t, x) ∈ [0, T ] × Rn .

(4.21)

Since uε is a subsolution of the nonlinear PDE (4.1), it follows thathiuε (t, x) ≤ E uε (t + h, X̂ht,x ) + h [F (., Duε )(t, x) + R(h, ε)] .

We next introduce the function U ε defined by U ε (1, x) = uε (1, x) andhi
U ε (ti , x) = E U ε (ti + h, X̂hti ,x ) + h [F (., DU ε )(ti , x) + R(h, ε)] .
Then, it follows from the comparison property of Lemma 4.3.1 that
uε

≤

U ε.

(4.22)

28

B. Bouchard and R. Elie and N. Touzi

Moreover, arguing as in the proof of (4.10), we see thathi
U ε − v h ≤ (1 + Ch)Et,x U ε − v h (t + h, X̂ht,x ) + hR(h, ε),which provides the estimate:U ε − vh

≤ R(h, ε).

(4.23)We then deduce from (4.20), (4.22) and (4.23) thatv − v h = v − uε + uε − v ≤ v − uε + U ε − v h ≤ C(ε) + R(h, ε),and thereforev − vhinf (C(ε) + R(h, ε)) .

≤ε>0(4.24)

The approximating subsolution The following construction of the function uε satisfying (4.20) requires that the nonlinearity F (t, x, r, p, γ) be concave in (r, p, γ) withsome additional conditions in order to ensure some regularity of the solution. Notice that the case where F (t, x, r, p, γ) is convex in (r, p, γ) can also be dealt with bya symmetric argument (inverting the roles of supersolutions and subsolutions). Fora bounded function ψ(t, x) Lipschitz in x and 1/2–Hölder continuous in t, we denote|ψ|1

Assumption HJB
F (t, x, r, p, γ):= |ψ|∞ +ψ(t, x) − ψ(t′ , x′ )′′ 1/2([0,1]×Rd )2 |x − x | + |t − t |supThe nonlinearity F is of the Hamilton–Jacobi–Bellman type:=Lα (t, x, r, p, γ) :=inf {Lα (t, x, r, p, γ)}α∈A1T r[σ α σ αT (t, x)γ] + bα (t, x)p + cα (t, x)r + f α (t, x)2

where b, σ , σ α , bα , cα and f α satisfy:|b|1 + |σ|1 + sup (|σ α |1 + |bα |1 + |cα |1 + |f α |1 ) < ∞.
α∈AThe above Assumption implies in particular that our nonlinear PDE satisfies a strong
comparison result for bounded functions. Let wε be the unique viscosity solution of
the nonlinear PDE with shaked coefficients (in the terminology of Krylov):

ε
ε
ε
2 ε
− ∂w
∂t − inf |e|≤1 F̄ t + εe, x + εe, w (t, x), Dw (t, x), D w (t, x) = 0,(4.25)wε (1, .) = g,where:F̄ (t, x, r, p, γ) :=1
hb(t, x), pi + Tr[σσ ′ γ] + F (t, x, r, p, γ).
229Discrete-time approximation of BSDEs

The existence of wε follows from a direct identification of this PDE as the HJB equationof a corresponding stochastic control problem. By the definition of wε together witha strong comparison argument, we then have thatwε is a subsolution of (4.1) and |wε − v| ≤ Cε,(4.26)

for some constant C > 0. It also follows from classical estimates that for a Lipschitz–
continuous final condition g :wε is bounded, Lipschitz in x, and 1/2–Hölder continuous in t.

(4.27)

Let ρ(t, x) be a C ∞ positive function supported in {(t, x) : t ∈ [0, 1], |x| ≤ 1} with unit
mass, and define


1
t x
ε
ε
ε
ε
(4.28)u (t, x) := w ∗ ρ where ρ (t, x) := d+2 ρ 2 ,εε εso that it follows from (4.27) thatuε∂tβ0 Dβ uε ≤ Cε1−2β0 −|β|1

is C ∞ , and(4.29)

P
for any (β0 , β) ∈ N × Nd \ {0}, where |β|1 := di=1 βi , and C > 0 is some constant.
Then, it follows from two successive applications of Itô’s formula thatRh [uε ]≤ R(h, ε) := Chε−3(4.30)

for some constant C . Finally, by the concavity of F , the function uε inherits (4.26):
uεis a viscosity subsolution of and |uε − v| ≤ C(ε) := Cε.

(4.31)εSince u satisfies all the requirements of the previous section, we deduce thatv − v h ≤ inf C ε + hε−3 ∼ Ch1/4 .
ε>0(4.32)The upper bound In the context of our scheme, we shall obtain an upper bound on
the error by inverting the roles of v and v h in the above lower bound argument. Let vεhbe the solution of the scheme with shaked coefficients:vεh (T, .) = g, vεh = Thε [vεh ], Thε [ψ] = inf Th,e [ψ], Th,e [ψ] := De0 [ψ] + hF (., De [ψ])
|e|≤ε

where
De ψ(t, x) = E [ψ(t + h, x + µ(t + e, x + e)h + σ(t + e, x + e)Wh )Hh (t + e, x + e)] .

Then, 0 ≤ v h − vεh by the comparison result for the scheme. When the volatility σ(.)
is chosen to be constant, we also have v h − vεh ≤ Cε for some constant C . Moreover,
it is proved in [37] that vεh is Lipschitz-continuous in x and 1/2–Hölder continuous int, uniformly in h. Then, the mollification vεh ∗ ρε of vεh satisfies the estimates (4.29).
Collecting all these ingredients, we obtain a lower bound on the error:v − vh≥ −C inf (ε + R(h, ε)) = −Ch1/4 .
ε30The main result

B. Bouchard and R. Elie and N. Touzi

Summing up the above results, we have

Theorem 4.3.2. Let the nonlinearity F be as in Assumption HJB, and let σ be constant.
Then, for any bounded Lipschitz final condition g :Ch1/10 ≤ v − v h ≤ Ch1/4

Remark 4.3.3. In the PDE Finite Differences literature, the rate of convergence isusually stated in terms of the discretisation in the space variable |∆x|. In our context of
stochastic differential equation, notice that |∆x| is or the order of h1/2 . Therefore, theabove upper bound on the rate of convergence corresponds to the classical rate |∆x|1/2 .

Bibliography[1] F. Antonelli and A. Kohatsu-Higa, Filtration stability of backward SDE’s, Stochastic Analysis
and Its Applications 18 (2000), pp. 11–37.
[2] V. Bally, An approximation scheme for BSDEs and applications to control and nonlinear
PDE’s, Pitman Research Notes in Mathematics Series, Longman 364 (1997).
[3] V. Bally and G. Pagès, A quantization algorithm for solving discrete time multidimensional
optimal stopping problems, Bernoulli 9 (2002), pp. 1003–1049.
[4] V. Bally, G. Pagès, and J. Printems, A quantization method for pricing and hedging multidimensional American style options, Mathematical Finance 15 (2005), pp. 119–168.
[5] G. Barles, R. Buckdahn, and E. Pardoux, Backward stochastic differential equations andintegral-partial differential equations, Stochastics and Stochastics Reports 60 (1997), pp. 57–
83.
[6] G. Barles and E. R. Jakobsen, On the convergence rate of approximation schemes for
Hamilton-Jacobi-Bellman equations, Mathematical Modelling and Numerical Analysis,
ESAIM, M2AM 36 (2002), pp. 33–54.
[7] G. Barles and E. R. Jakobsen, Error bounds For monotone approximation schemes for
Hamilton-Jacobi-Bellman equations, SIAM J. Numer. Anal. 43 (2005), pp. 540–558.
[8] G. Barles and E. R. Jakobsen, Error bounds For monotone approximation schemes forparabolic Hamilton-Jacobi-Bellman equations, Math. Comp. 76 (2007), pp. 1861–1893.
[9] G. Barles and P. E. Souganidis, Convergence of approximation schemes for fully non-linear
second order equation, Asymptotic Analysis 4 (1991), pp. 271–283.
[10] D. Becherer, Bounded solutions to backward SDE’s with jumps for utility optimization andindifference hedging, preprint, Imperial College, 2005.
[11] C. Bender and J. Zhang, Time discretization and markovian iteration for coupled FBSDEs,
Ann. Appl. Probab. 18 (2008), pp. 143–177 .
[12] J.M. Bismut, Conjugate convex functions in optimal stochastic control, J. Math. Anal. Appl.
44 (1973), pp. 384–404.
[13] F. Bonnans, Quelques aspects numériques de la commande optimale stochastique, Lecture
Notes, 2009.

Discrete-time approximation of BSDEs

31[14] F. Bonnans and H. Zidani, Consistency of generalized finite difference schemes for thestochastic HJB equation, SIAM J. Numerical Analysis 41 (2003),pp. 1008–1021.
[15] B. Bouchard, A stochastic target formulation for optimal switching problems in finite horizon,
to appear in Stochastics.
[16] B. Bouchard and J.-F. Chassagneux, Discrete time approximation for continuously and discretely reflected BSDEs, Stochastic Processes and their Applications 118 (2008), pp. 2269–
2293.
[17] B. Bouchard, I. Ekeland, and N. Touzi, On the Malliavin approach to Monte Carlo approximation of conditional expectations, Finance and Stochastics 8 (2004), pp. 45–71.
[18] B. Bouchard and R. Elie, Discrete time approximation of decoupled forward-backward SDEwith jumps, Stochastic Processes And Their Applications 118 (2008), pp. 53–75[19] B. Bouchard and S. Menozzi, Strong approximation of BSDEs in a domain, to appear inBernoulli.
[20] B. Bouchard and N. Touzi, Discrete-time approximation and Monte-Carlo simulation of backward stochastic differential equations, Stochastic Processes And Their Applications 111(2004), pp. 175–206.
[21] P. Briand, B. Delyon, and J. Mémin, Donsker-type theorem for BSDE’s, Electronic Communications In Probability 6 (2001), pp. 1–14.
[22] J. F. Carrière, Valuation of the early-exercise price for options using simulations and nonparametric regression, Insurance : Mathematics And Economics 19 (1996), pp. 19–30.
[23] J.-F. Chassagneux, Discrete-time approximation of doubly reflected BSDEs, Advances in Applied Probability 41 (2009).
[24] J.-F. Chassagneux, Processus réfléchis en finance et probabilités numérique, PhD Thesis,
University Paris-Diderot, 2008.
[25] P. Cheridito, H. M. Soner, N. Touzi, and N. Victoir, Second order backward stochastic differential equations and fully non-linear parabolic PDEs, Communications On Pure And Applied
Mathematics 60 (2007), pp. 1081–1110.
[26] D. Chevance, Numerical methods for backward stochastic differential equations, in Numerical methods in finance 1997, Edt L.C.G. Rogers and D. Talay, Cambridge University Press,
pp. 232–244.
[27] E. Clément, D. Lamberton, and P. Protter, An analysis of a least squares regression method
for American option pricing, Finance And Stochastics 6 (2002), pp. 449–472.
[28] F. Coquet, J. Mémin, and L. Slominski, On weak convergence of filtrations, Séminaire De
Probabilités, XXXV, pp. 306–328. Lecture Notes in Math. 1755 (2001) Springer, Berlin.
[29] F. Coquet, V. Mackevičius, and J. Mémin, Stability in D of martingales and backward equations under discretization of filtration, Stochastic Processes And Their Applications 75(1998), pp. 235–248.
[30] R. W. R. Darling and E. Pardoux, BSDE with random terminal time, Annals Of Probability
25 (1997), pp. 1135–1159.
[31] F. Delarue and S. Menozzi, A forward backward algorithm for quasi-linear PDEs, Annals
Of Applied Probability 16 (2006), pp. 140–184.
[32] F. Delarue and S. Menozzi, An interpolated stochastic algorithm for quasi-linear PDEs,
Mathematics Of Computation 261 (2008), pp. 125–158.

32

B. Bouchard and R. Elie and N. Touzi[33] H. Dong, N. V. Krylov, On the rate of convergence of finite-difference approximations forbellman’s equations with constant coefficients, St. Petersburg Math. J. 17 (2005), pp. 108–132.
[34] J. Jr. Douglas, J. Ma, and P. Protter, Numerical methods for forward- backward stochastic
differential equations, Annals Of Applied Probability 6 (1996), pp. 940–968.
[35] R. Elie, Contrôle stochastique et méthodes numériques en finance mathématique, PHD thesis,
University Paris-Dauphine, 2006.
[36] A. Eyraud-Loisel, Backward stochastic differential equations with enlarged filtration. option
hedging of an insider trader in a financial market with jumps, Stochastic Processes And Their
Applications 115 (2005), pp. 1745–1763,[37] A. Fahim, N. Touzi, and X. Warin, A probabilistic numerical scheme for fully nonlinear PDEs,
preprint, 2009.
[38] E. Fournié, J.-M. Lasry, J. Lebuchoux, P.-L. Lions, and N. Touzi, Applications of Malliavin
calculus to Monte Carlo methods in finance, Finance And Stochastics 3 (1999), pp. 391–412.
[39] E. Gobet, Schéma d’Euler pour diffusions tues. Application aux options barrire, PhD Thesis,
University Paris-Diderot, 1998.
[40] E. Gobet, Weak approximation of killed diffusion using Euler schemes, Stochastic Processes
And Their Applications 87 (2000), pp. 167–197.
[41] E. Gobet, Revisiting the Greeks for European and American options, Proceedings of the ”International Symposium on Stochastic Processes and Mathematical Finance” at Ritsumeikan
University, Kusatsu, Japan, March 2003. Edited by J. Akahori, S. Ogawa, S. Watanabe. World
Scientific, pp. 53–71, 2004.
[42] E. Gobet and C. Labart, Error expansion for the discretization of backward stochastic differential equations, Stochastic Processes And Applications 117 (2007), pp. 803–829.
[43] E. Gobet, J.P. Lemor, and X. Warin, Rate of convergence of empirical regression method forsolving generalized BSDE, Bernoulli 12 (2006), pp. 889–916.
[44] E. Gobet and A. Makhlouf, L2 -time regularity of BSDEs with irregular terminal functions,preprint 2008.
[45] E. Gobet and S. Menozzi, Stopped diffusion processes: overshoots and boundary correction,
preprint PMA, University Paris-Diderot, 2007.
[46] Y. Hu, P. Imkeller, and M. Müller, Utility maximization in incomplete markets, Ann. Appl.
Probab. 15 (2005), pp. 1691–1712.
[47] N. El Karoui, C. Kapoudjan, E. Pardoux, S. Peng, and M.C. Quenez, Reflected solutions ofbackward stochastic differential equations and related obstacle problems for PDE’s, Annalsof Probability 25 (1997), pp. 702–737.
[48] N. El Karoui, S. Pend and M.-C. Quenez, Backward stochastic differential equations in finance, Mathematical Finance 7 (1997), pp. 1–71.
[49] N. El Karoui, S. Peng, and M.-C. Quenez, A dynamic maximum principle for the optimization of recursive utilities under constraints, The Annals Of Applied Probability 11 (2001),
pp. 664–693[50] N. El Karoui and R. Rouge, Pricing via utility maximization and entropy, Mathematical Finance 10 (2000), pp. 259–276.
[51] I. Kharroubi, H. Pham, J. Ma, and J. Zhang, Backward SDEs with constrained jumps andquasi-variational inequalities, preprint 2008.

Discrete-time approximation of BSDEs

33[52] N. V. Krylov, On the rate of convergence of finite-difference approximations for bellman’sequations, St. Petersburg Math. J. 9 (1997), pp. 245–256.
[53] N. V. Krylov, On the rate of convergence of finite-difference approximations for bellman’sequations with variable coefficients, Probability Theory And Related Fields 117 (2000),
pp. 1–16.
[54] N. V. Krylov, On the rate of convergence of finite-difference approximations for Bellman’sequations with Lipschitz coefficients, Applied Math. And Optimization 52 (2005), pp. 365–399.
[55] P.-L. Lions and H. Regnier, Calcul du prix et des sensibilités d’une option américaine par uneméthode de Monte Carlo, preprint 2001.
[56] F. A. Longstaff and R. S. Schwartz, Valuing american options by simulation : a simple leastsquare approach, Review Of Financial Studies 14 (2001), pp. 113–147.
[57] J. Ma, P. Protter, J. San Martin, and S. Torres, Numerical method for backward stochastic
differential equations, Annals Of Applied Probability 12 (2002), pp. 302–316.
[58] J. Ma, P. Protter, and J. Yong, Solving forward-backward stochastic differential equations
explicitly - a four step scheme, Probability Theory And Related Fields 98 (1994), pp. 339–
359.
[59] J. Ma and J. Zhang, Path regularity of solutions to backward stochastic differential equations,
Probability Theory And Related Fields 122 (2002), pp. 163–190.
[60] J. Ma and J. Zhang, Representations and regularities for solutions to BSDEs with reflections,
Stochastic Processes And Their Applications 115 (2005), pp. 539–569.
[61] D. Nualart, The Malliavin Calculus and Related Topics, Springer, Berlin, 1995.
[62] E. Pardoux, Backward stochastic differential equations and viscosity solutions of semilinear
parabolic and elliptic PDE’s of second order, In Stochastic Analysis and Related Topics:
The Geilo Workshop 1996. L. Decreusefond, J. Gjerd, B. Oksendal, and A.S. Ustnel (eds.),
Birkhäuser, pp. 79–127, 1998.
[63] E. Pardoux and S. Peng, Adapted solution of a backward stochastic differential equation,
Systems and control letters 14 (1990), pp. 55–61.
[64] E. Pardoux and S. Peng, Backward stochastic differential equations and quasilinear parabolic
partial differential equations, Lecture Notes In Control And Inform. Sci. 176 (1992), pp. 200–217.
[65] E. Pardoux, F. Pradeilles, and Z. Rao, Probabilistic interpretation for a system of semilinear
parabolic partial differential equations, Ann. Inst. H. Poincare 33 (1997), pp. 467–490.
[66] S. Peng, Probabilistic interpretation for systems of quasilinear parabolic partial differential
equations, Stochastics And Stochastics Reports 37 (1991), pp. 61–74.
[67] S. Peng, Backward stochastic differential equations and applications to optimal control, Appl
Math Optim 27 (1993), pp. 125–144.
[68] S. Peng, Backward stochastic differential equations, nonlinear expectations, nonlinear evaluations and risk measures, Lecture Notes in Chinese Summer School in Mathematics Weihai,
July 19–August 14, 2004.
[69] S. Peng, G-Brownian motion and dynamic risk measure under volatility uncertainty,
arXiv:0711.2834v1, 2007.
[70] S. Rubenthaler, Numerical simulation of the solution of a stochastic differential equation
driven by a Levy process, Stochastic Processes And Their Applications 103 (2003), pp. 311–
349.

34

B. Bouchard and R. Elie and N. Touzi[71] A. B. Sow and E. Pardoux, Probabilistic interpretation of a system of quasilinear parabolic
PDEs, Stochastics And Stochastics Reports 76 (2004), pp. 429–477.
[72] S. Tang and X. Li, Necessary conditions for optimal control of stochastic systems with random
jumps, SIAM J. Control Optim. 32 (1994), pp. 1447–1475.
[73] J. Zhang, Some fine properties of backward stochastic differential equations, PhD thesis, Purdue University, 2001.
[74] J. Zhang, A numerical scheme for BSDEs, Annals of Applied Probability 14 (2004), pp. 459–
488.

Author information
Bruno Bouchard, CEREMADE, Université Paris Dauphine and CREST-ENSAE, France.
Email: bouchard@ceremade.dauphine.fr
Romuald Elie, CEREMADE, Université Paris Dauphine and CREST-ENSAE, France.
Email: elie@ceremade.dauphine.
Nizar Touzi, CMAP, Ecole Polytechnique Paris, France.
Email: nizar.touzi@polytechnique.edu