A Forward Scheme for Backward SDEs

Universit¨ at Konstanz

A Forward Scheme for Backward SDEs

Christian Bender Robert Denk

Konstanzer Schriften in Mathematik und Informatik Nr. 209, August 2005

ISSN 1430–3558

c

Fachbereich Mathematik und Statistik c

Fachbereich Informatik und Informationswissenschaft Universit¨at Konstanz

Fach D 188, 78457 Konstanz, Germany Email: preprints@informatik.uni–konstanz.de

Konstanzer Online-Publikations-System (KOPS) URL: http://www.ub.uni-konstanz.de/kops/volltexte/2006/2229/

URN: http://nbn-resolving.de/urn:nbn:de:bsz:352-opus-22297

A Forward Scheme for Backward SDEs ¹

Christian Bender and Robert Denk WIAS Berlin and University of Konstanz

Abstract

We introduce a forward scheme to simulate backward SDEs. Com- pared to existing schemes, we avoid high order nestings of conditional expectations backwards in time. In this way the error, when approxi- mating the conditional expectation, in dependence of the time partition is significantly reduced. Besides this generic result, we present an im- plementable algorithm and provide an error analysis for it. Finally, we demonstrate the strength of the new algorithm by solving some financial problems numerically.

AMS classifications:Primary 65C05; Secondary 65C30, 91B28.

Keywords: BSDE, Numerics, Monte-Carlo simulation, Finance.

1 Introduction

The study of nonlinear backward stochastic differential equations (BSDEs) was initiated by Pardoux and Peng (1990). Mainly motivated by financial problems (see e.g. the survey article by El Karoui et al. (1997)) the theory of BSDEs was developed at high speed during the 1990s. Comparably slow progress has been made on the numerics of BSDEs.

Up to now basically two types of schemes have been considered. Based on the theoretical 4-step-scheme from Ma et al. (1994), numerical algorithms for BSDEs have been developed by Douglas et al. (1996) and more recently by Milstein and Tretyakov (2004). The main focus of these algorithms is the numerical solution of a parabolic PDE which is related to the BSDE.

A second type of algorithms works backwards through time and tries to tackle the stochastic problem directly. Bally (1997) and Chevance (1997) were the first to study this type of algorithm with a (hardly implementable) random time partition respectively under strong regularity assumptions. The work of Ma et al. (2002) is in the same spirit, replacing, however, the Brownian motion by a binary random walk in the approximative equation. See also Briand et al.

(2001) for the binary random walk approach. Only recently, a new notion of L²-regularity on the control part of the solution was introduced in Zhang (2004), which allowed to prove convergence of this backward approach with determin- istic partitions under rather weak regularity assumptions, see Zhang (2004), Bouchard and Touzi (2004), and Gobet et al. (2004) for different algorithms.

1C. Bender is supported by the DFG Research CenterMatheon‘Mathematics for key tech- nologies’ in Berlin. R. Denk is partially supported by the AFF grant 28/04 of the University of Konstanz. Some parts of the present paper are contained in the preprint “Forward Simula- tion of Financial Problems via BSDEs”,Konstanzer Schriften in Mathematik und Informatik 207, 2005.

A main drawback of the backward schemes is, that nestings of conditional expectations backwards through the time steps have to been evaluated. For a practical implementation the conditional expectations must be replaced by some estimator. A generic result of Bouchard and Touzi (2004) shows that the error due to the approximation of the conditional expectation grows with order 1/2, as the number of time steps goes to infinity. This leads to high computational costs, when a fine mesh of the time discretization is required.

In this paper we propose a new forward scheme, which avoids nestings of con- ditional expectations backwards through the time steps. Instead it mimics the Picard type iteration for BSDEs and, consequently, has nestings of conditional expectation along the iterations.

We show that the additional error due to the iteration converges to zero at geometric rate (Theorem 2.6). At this cost the error, when approximating the conditional expectations by a generic estimator, in dependence of the time par- tition is reduced by order 1/2 compared to existing backward schemes (Theorem 3.1). In fact, in our scheme this error does neither explode when the number of time steps nor when the number of iterations tends to infinity. We believe that this is a striking advantage compared to the backward scheme.

Besides this generic results, we develop a practically implementable numeri- cal scheme. In particular, we use the regression-based least squares Monte-Carlo method to approximate the conditional expectation as was suggested by Gobet et al. (2004) in the context of the backward scheme. We analyze the error, when replacing the conditional expectation by the orthogonal projections on subspaces (Theorem 4.1), and also provide rates of convergence when the projection co- efficients are substituted by their simulation-based analogues (Theorem 4.9).

Again we have an error reduction of order 1/2 in the mesh size of the time partition compared to the results in Gobet et al. (2004). Depending on the number Lof simulated paths, the best expected rate of L^−1/2can be achieved for appropriate projection spaces (Theorem 4.9).

Finally, we present some simulations related to financial problems (Section 5). We consider the hedging problem under different interest rates for invest- ing and borrowing and the superhedging problem under borrowing constraints, which lead to nonlinear BSDEs.

2 A Discretization of the Picard Type Iteration

In this section we introduce a discretized Picard iteration and prove its conver- gence for the following type of BSDE:

dXt = b(t, Xt)dt+σ(t, Xt)dWt

dYt = f(t, Xt, Yt, Zt)dt+ZtdWt

X0 = x YT = ξ.

Here Wt= (W1,t, . . . , WD,t)^∗, (the star denoting matrix transposition), is aD- dimensional Brownian motion on [0, T] andZt= (Z1,t, . . . , ZD,t). The process X is R^M-valued and the process Y is R-valued. Throughout the paper we assume:

Assumption 2.1. There is a constantK such that

|b(t, x)−b(t⁰, x⁰)|+|σ(t, x)−σ(t⁰, x⁰)|+|f(t, x, y, z)−f(t⁰, x⁰, y⁰, z⁰)|

≤ K(p

|t−t⁰|+|x−x⁰|+|y−y⁰|+|z−z⁰|) for all(t, x, y, z), (t⁰, x⁰, y⁰, z⁰)∈[0, T]×R^M ×R×R^D,

ξ= Φ(X)

where Φ is a functional on the space of R^M-valued RCLL-functions on [0, T] satisfying theL^∞-Lipschitz condition,

|Φ(x)−Φ(x⁰)| ≤K sup

0≤t≤T|x(t)−x⁰(t)| for all RCLL-functions x, x⁰. Moreover,

sup

0≤t≤T

(|b(t,0)|+|σ(t,0)|+|f(t,0,0,0)|) +|Φ(0)| ≤K where0denotes the constant function taking value 0 on [0, T].

Note, that we do neither assume that the matrixσis quadratic nor thatσσ^∗ is invertible.

Remark 2.2. We shall say that a constant depends on the data, if it depends on K, T, x0 and the dimensions M and D only. Throughout the paper C denotes a generic constant depending on the data which may vary from line to line.

Theoretically, the backward part (Y, Z) can be obtained as the limit of a Picard type iteration (Y⁽ⁿ⁾, Z⁽ⁿ⁾), see e.g. Yong and Zhou (2000), theorem 7.3.4. Here (Y⁽⁰⁾, Z⁽⁰⁾)≡(0,0), and (Y⁽ⁿ⁾, Z⁽ⁿ⁾) is the solution of the simple BSDE

dY_t⁽ⁿ⁾ = f(t, Xt, Y_t⁽ⁿ⁻¹⁾, Z_t⁽ⁿ⁻¹⁾)dt+Z_t⁽ⁿ⁾dWt

Y_T⁽ⁿ⁾ = ξ withX as above.

The solution is given by Y_t⁽ⁿ⁾=E

"

ξ− Z T

t

f(s, Xs, Y_s⁽ⁿ⁻¹⁾, Z_s⁽ⁿ⁻¹⁾)ds

¯¯

¯¯

¯Ft

#

and Z⁽ⁿ⁾ is obtained via the martingale representation theorem. As is empha- sized in Yong and Zhou (2000), ch. 7, the above Picard iteration is still implicit due to the use of the martingale representation theorem.

We will now introduce a time discretization of the above iteration, which is explicit but for the occurrence of conditional expectations.

Suppose a partition π = {t0, t1, . . . , tN} of [0, T] with mesh size |π| :=

maxi|ti+1−ti| is given and a corresponding discretization X^(π) of X as well

as some approximation ξ^(π) of ξ. Let (Y^(0,π), Z^(0,π)) ≡ (0,0). Then define iteratively fori= 0,1, . . . , N, with ∆i=ti+1−tiand ∆Wd,i=Wd,ti+1−Wd,ti,

Y_t^(n,π)_i = E



ξ^(π)−

NX−1 j=i

f(tj, X_t^(π)_j , Y_t^(n−1,π)_j , Z_t^(n−1,π)_j )∆j

¯¯

¯¯

¯¯F^ti



,

Z_d,t^(n,π)_i = E



∆Wd,i

∆i



ξ^(π)−

N−1X

j=i+1

f(tj, X_t^(π)_j , Y_t^(n−1,π)_j , Z_t^(n−1,π)_j )∆j





¯¯

¯¯

¯¯Fti



,

d= 1, . . . , D. (Here we used the convention ∆Wd,N = 0). The processesY^(n,π) andZ^(n,π)are extended to RCLL processes by constant interpolation. Note that the discretized Picard-type iteration has no nestings of conditional expectations backward in time, but forward in the number of Picard iterations. This turns out to be an advantage from the numerical point of view (see section 3 below).

We can now state convergence of the discretized Picard-type iteration:

Theorem 2.3. Suppose Assumption 2.1 holds, and for some constant C de- pending on the data

sup

0≤t≤T

Eh

|Xt−X_t^(π)|²i

≤ C|π|, sup

|π|≤1

Eh

|ξ^(π)|²i

≤ C.

Then there is a constant C depending on the data such that sup

0≤t≤T

E·¯¯¯Yt−Y_t^(n,π)¯¯¯

2¸ +E

Z T

0 |Zt−Z_t^(n,π)|²dt

≤ C

µ

|π|+E[|ξ−ξ^(π)|²] + µ1

2+C|π|

¶n¶

provided|π| is sufficiently small.

Remark 2.4. (i) Note, the condition on the discretization X^(π) of X is, for instance, satisfied by the Euler scheme.

(ii) The condition onξ^(π) is satisfied, whenever for|π| ≤1 E[|ξ−ξ^(π)|²]≤C|π|^α

with some constantC depending on the data and someα >0. Indeed, E[|ξ−ξ^(π)|²]≤2E[|ξ|²] + 2E[|ξ−ξ^(π)|²],

and, thanks to theL^∞-Lipschitz condition and a classical estimate for SDEs, E[|ξ|²] ≤ 2K²E[ sup

0≤t≤T|Xt|²] + 2|Φ(0)|²

≤ C

Ã

|x|²+ Z T

0 |b(t,0)|²+|σ(t,0)|²dt

!

+ 2K²≤C.

The proof of theorem 2.3 is split into two parts. Given the partition π and a corresponding discretizationX^(π)ofX we define (Y^(∞,π), Z^(∞,π)) as the solution of

Y_t^(∞,π)_N = ξ^(π), Z_d,t^(∞,π)_i = E

·∆Wd,i

∆i

Y_t^(∞,π)_i+1

¯¯

¯¯Fti

¸ ,

Y_t^(∞,π)_i = E[Y_t^(∞,π)_i+1 |Fti]−f(ti, X_t^(π)_i , Y_t^(∞,π)_i , Z_t^(∞,π)_i )∆i.

It exists, when the mesh |π| of the partition π is sufficiently fine. Again, the processesY^(∞,π) and Z^(∞,π) are extended to RCLL processes by constant in- terpolation. Note, (Y^(∞,π), Z^(∞,π)) is (up to the interpolation of the Z-part) the backward scheme considered in Bouchard and Touzi (2004). We remark that this backward scheme is still implicit, and inner iterations are required for numerical implementation.

We shall separately consider the convergence of (Y^(n,π), Z^(n,π)) to (Y^(∞,π), Z^(∞,π)) and of (Y^(∞,π), Z^(∞,π)) to (Y, Z).

Concerning the backward scheme we need an extension of the results by Bouchard and Touzi (2004). The following variant of theorem 3.1 in Bouchard and Touzi (2004) is a slight generalization concerning the assumptions on the coefficients. Moreover, it allows for path-depending terminal data and the ap- proximating processes are piecewise constant.

Theorem 2.5. Suppose Assumption 2.1 holds, and the discretization X^(π) of X satisfies

sup

0≤t≤T

Eh

|Xt−X_t^(π)|²i

≤C|π| (1)

for some constantCdepending on the data. Then there is a constantC depend- ing on the data such that

sup

0≤t≤T

E·¯¯¯Yt−Y_t^(∞,π)¯¯¯

2¸ +E

Z T

0 |Zt−Z_t^(∞,π)|²dt

≤ C³

|π|+E[|ξ−ξ^(π)|²]´ provided|π| is sufficiently small.

The proof combines ideas of Bouchard and Touzi (2004) and Zhang (2004), who suggests a different time discretization. For the reader’s convenience we sketch the proof of Theorem 2.5 in the Appendix.

We now investigate the iteration for a fixed partition. Our aim is to derive rates of convergence uniform inπ.

Theorem 2.6. Under the assumptions of theorem 2.3 there are constants C1

andC2 depending on the data such that

0≤i≤Nmax E·¯¯¯Y_t^(∞,π)_i −Y_t^(n,π)_i ¯¯¯

2¸ +

NX−1 i=0

E·¯¯¯Z_t^(∞,π)_i −Z_t^(n,π)_i ¯¯¯

2¸

∆i

≤ C1

µ1

2+C2|π|

¶n

provided|π| is sufficiently small.

Clearly, Theorem 2.3 follows from a straightforward combination of Theo- rems 2.5 and 2.6.

Remark 2.7. LetK denote the Lipschitz constant off. Then it follows from the proof below that Theorem 2.6 holds, for instance, for|π| ≤Γ with

C2 = Γ 4, where

Γ = 16T(T+ 1)²D²K⁴+ 4K(T+ 1)K².

We prepare the proof of Theorem 2.6 with some a priori estimates.

Lemma 2.8. SupposeΓandγ are positive real numbers,y˜^(ι),z˜^(ι),ι= 1,2 are adapted processes and

Y˜_t^(ι)_i = E



ξ^(π)−

N−1X

j=i

f(tj, X_t^(π)_j ,y˜_t^(ι)_j ,z˜^(ι)_tj )∆j

¯¯

¯¯

¯¯F^ti



,

Z˜_d,t^(ι)_i = E



∆Wd,i

∆i



ξ^(π)−

NX−1 j=i+1

f(tj, X_t^(π)_j ,y˜^(ι)_tj ,z˜_t^(ι)_j )∆j





¯¯

¯¯

¯¯Fti



.

Moreover, assume thatf is Lipschitz in(y, z) uniformly in(t, x) with constant K. Then:

0≤i≤Nmax λiEh

|Y˜_t⁽¹⁾_i −Y˜_t⁽²⁾_i |²i +

NX−1 i=0

λiEh

|Z˜_t⁽¹⁾_i −Z˜_t⁽²⁾_i |²i

∆i

≤ K²(T+ 1) µ¡

|π|+ Γ⁻¹¢

(γDT+ 1) + D γ

¶

× Ã1

T

NX−1 i=0

λiEh

|y˜_t⁽¹⁾_i −y˜⁽²⁾_ti |²i

∆i+

NX−1 i=0

λiEh

|z˜_t⁽¹⁾_i −z˜_t⁽²⁾_i |²i

∆i

! ,

whereλ0= 1 andλi= (1 + Γ∆i−1)λi−1.

Proof. The proof goes through several steps. For notational convenience let us introduce

yt_i = y˜⁽¹⁾_ti −y˜_t⁽²⁾_i , zti = z˜_t⁽¹⁾_i −z˜_t⁽²⁾_i ,

∆fi = f(ti, X_t^(π)_i ,y˜_t⁽¹⁾_i ,z˜_t⁽¹⁾_i )−f(ti, X_t^(π)_i ,y˜⁽²⁾_ti ,z˜_t⁽²⁾_i ).

First note that

Y˜_t^(ι)_i =E[ ˜Y_t^(ι)_i+1|Fti]−f(ti, X_t^(π)_i ,y˜_t^(ι)_i ,z˜_t^(ι)_i )∆i (2) and, for thedth component of ˜Z^(ι),

Z˜_d,t^(ι)_i =E

·∆Wd,i

∆i

Y˜_t^(ι)_i+1

¯¯

¯¯F^ti

¸

. (3)

Step 1: We prove that for any 1≤d≤D

NX−1 i=0

λiE·¯¯¯Z˜_d,t⁽¹⁾_i−Z˜_d,t⁽²⁾_i¯¯¯

2¸

∆i

≤ γ

NX−1 i=0

λiE·¯¯¯Y˜_t⁽¹⁾_i −Y˜_t⁽²⁾_i ¯¯¯

2¸

∆i+(1 +T)K² γ

NX−1 i=0

λiEh

|zti|²i

∆i

+(1 +T)K² T γ

NX−1 i=0

λiEh

|yt_i|²i

∆i. (4)

First note that by (3) and H¨older’s inequality, Z˜_d,t⁽¹⁾

i−Z˜_d,t⁽²⁾

i = E

·∆Wd,i

∆i

³Y˜_t⁽¹⁾_i+1−Y˜_t⁽²⁾_i+1´¯¯¯¯Fti

¸

= E

·∆Wd,i

∆i

³Y˜_t⁽¹⁾_i+1−Y˜_t⁽²⁾_i+1−E[ ˜Y_t⁽¹⁾_i+1−Y˜_t⁽²⁾_i+1|Fti]´¯¯¯¯Fti

¸

≤ r 1

∆i

E· ³

Y˜_t⁽¹⁾_i+1−Y˜_t⁽²⁾_i+1−E[ ˜Y_t⁽¹⁾_i+1−Y˜_t⁽²⁾_i+1|Fti]´2¯¯¯¯Fti

¸1/2

. Thus, by (2),

Eh

|Z˜_d,t⁽¹⁾_i−Z˜_d,t⁽²⁾_i|²i

≤ 1

∆i

Eh

|Y˜_t⁽¹⁾_i+1−Y˜_t⁽²⁾_i+1|²−E[ ˜Y_t⁽¹⁾_i+1−Y˜_t⁽²⁾_i+1|Fti]²i

= 1

∆i

Eh

|Y˜_t⁽¹⁾_i+1−Y˜_t⁽²⁾_i+1|²− |Y˜_t⁽¹⁾_i −Y˜_t⁽²⁾_i + ∆fi∆i|²i

≤ 1

∆iEh

|Y˜_t⁽¹⁾_i+1−Y˜_t⁽²⁾_i+1|²− |Y˜_t⁽¹⁾_i −Y˜_t⁽²⁾_i |²−2( ˜Y_t⁽¹⁾_i −Y˜_t⁽²⁾_i )∆fi∆i

i . Multiplying both sides with the weights λi∆i and summing from 0 to N−1 yields forγ >0,

N−1X

i=0

λiEh

|Z˜_t⁽¹⁾_i −Z˜_t⁽²⁾_i |²i

∆i+λ0Eh

|Y˜_t⁽¹⁾₀ −Y˜_t⁽²⁾₀ |²i

≤ λNEh

|Y˜_t⁽¹⁾_N −Y˜_t⁽²⁾_N |²i

−2

NX−1 i=0

λiEh

( ˜Y_t⁽¹⁾_i −Y˜_t⁽²⁾_i )∆fi

i∆i

≤ γ

N−1X

i=0

λiEh

|Y˜_t⁽¹⁾_i −Y˜_t⁽²⁾_i |²i

∆i+K² γ

N−1X

i=0

λiEh

(|yti|+|zti|)²i

∆i. Here we used ˜Y_t⁽¹⁾_N −Y˜_t⁽²⁾_N = 0 and Young’s inequality. (4) may now be obtained by another application of Young’s inequality.

Step 2: We show

0≤i≤Nmax λiEh

|Y˜_t⁽¹⁾_i −Y˜_t⁽²⁾_i |²i

≤ K²(T+ 1) µ

|π|+ 1 Γ

¶ Ã^NX⁻¹

i=0

λiE£

|zti|²∆i¤ +1

T

NX−1 i=0

λiE£

|yti|²∆i¤! (5)

By (2), Jensen’s inequality, and Young’s inequality we get Eh

|Y˜_t⁽¹⁾_j −Y˜_t⁽²⁾_j |²i

≤ (1 + Γ∆j)Eh

|Y˜_t⁽¹⁾_j+1−Y˜_t⁽²⁾_j+1|²i

+ (∆j+ Γ⁻¹)E[(∆fj)²]∆j

≤ (1 + Γ∆j)Eh

|Y˜_t⁽¹⁾_j+1−Y˜_t⁽²⁾_j+1|²i +¡

|π|+ Γ⁻¹¢

K²(T + 1)E[|zt_j|²]∆j

+¡

|π|+ Γ⁻¹¢

K²T + 1

T E[|yt_j|²]∆j.

Multiplying with λj and summing from j =ito N −1 easily yields (5), since Y˜_t⁽¹⁾_N −Y˜_t⁽²⁾_N = 0.

Final Step: The assertion follows from a straightforward combination of (4) and (5).

Proof of theorem 2.6. Denote,

y_t^(n+1,π)_i = Y_t^(n+1,π)_i −Y_t^(n,π)_i , z_t^(n+1,π)_i = Z_t^(n+1,π)_i −Z_t^(n,π)_i . By Lemma 2.8,

0≤i≤Nmax λiEh

|y^(n+1,π)_ti |²i +

N−1X

i=0

λiEh

|z_t^(n+1,π)_i |²i

∆i

≤ K²(T+ 1) µ¡

|π|+ Γ⁻¹¢

(γDT + 1) +D γ

¶

× Ã

0≤i≤Nmax λiEh

|y_t^(n,π)_i |²i +

N−1X

i=0

λiEh

|z_t^(n,π)_i |²i

∆i

! .

We now choose γ= 4DK²(T+ 1) and Γ = 4K²(T + 1)(γDT + 1) and iterate the above inequality to obtain,

0≤i≤Nmax λiEh

|y_t^(n+1,π)_i |²i +

NX−1 i=0

λiEh

|z^(n+1,π)_ti |²i

∆i

≤

µΓ|π| 4 +1

2

¶nÃ

0≤i≤Nmax λiEh

|Y_t^(1,π)_i |²i +

N−1X

i=0

λiEh

|Z_t^(1,π)_i |²i

∆i

! .

Recalling the definition ofλi from Lemma 2.8 we have,

0≤i≤Nmax Eh

|y^(n+1,π)_ti |²i +

NX−1 i=0

Eh

|z_t^(n+1,π)_i |²i

∆i

≤ e^ΓT µΓ|π|

4 +1 2

¶nÃ

0≤i≤Nmax Eh

|Y_t^(1,π)_i |²i +

NX−1 i=0

Eh

|Z_t^(1,π)_i |²i

∆i

! .

Denote the square root of the right-hand side by A(π, n). Clearly the se- ries P

nA(π, n) converges, when |π| is sufficiently small. This shows, that

(Y^(n,π), Z^(n,π)) is a Cauchy sequence and thus converges to (Y^(∞,π), Z^(∞,π)) (when|π|is sufficiently small) by means of (2)–(3). Moreover, forn∈N,

0≤i≤Nmax E·¯¯¯Y_t^(∞,π)_i −Y_t^(n,π)_i ¯¯¯²

¸ +

N−1X

i=0

E·¯¯¯Z_t^(∞,π)_i −Z_t^(n,π)_i ¯¯¯²

¸

∆i

≤ Ã_∞

X

ν=n

A(π, ν)

!2

≤ e^ΓT Ã

0≤i≤Nmax Eh

|Y_t^(1,π)_i |²i +

N−1X

i=0

Eh

|Z_t^(1,π)_i |²i

∆i

! Ã 1−

rΓ|π| 4 +1

2

!−2

× µΓ|π|

4 +1 2

¶n

.

It remains to prove a uniform bound for Ã

0≤i≤Nmax Eh

|Y_t^(1,π)_i |²i +

NX−1 i=0

Eh

|Z_t^(1,π)_i |²i

∆i

! ,

which is given in the following lemma.

Lemma 2.9. Under the assumptions of theorem 2.3, there is a constant C depending on the data only such that

0≤i≤Nmax Eh

|Y_t^(1,π)_i |²i +

NX−1 i=0

Eh

|Z_t^(1,π)_i |²i

∆i≤C provided|π| ≤1 .

Proof. By Young’s and H¨older’s inequality we have

0≤i≤Nmax Eh

|Y_t^(1,π)_i |²i

≤2E[|ξ^(π)|²] + 2T

N−1X

j=0

Eh

|f(tj, X_t^(π)_j ,0,0)|²i

∆j

The first term on the right hand side is bounded by a constant depending on the data for|π| ≤1 by assumption. For the second term we observe

Eh

|f(tj, X_t^(π)_j ,0,0)|²i

≤ 2Eh

|f(tj, X_t^(π)_j ,0,0)−f(tj,0,0,0)|²i

+ 2|f(tj,0,0,0)|²

≤ 2K² µ

sup

0≤t≤T

E[|X_t^(π)|²] + 1

¶ . Now, by assumption and a classical result on SDEs

sup

0≤t≤T

E[|X_t^(π)|²]≤2 sup

0≤t≤T

E[|X_t^(π)−Xt|²] + 2 sup

0≤t≤T

E[|Xt|²]

≤ C|π|+C Ã

|x|²+ Z T

0 |b(t,0)|²+|σ(t,0)|²dt

!

≤C(1 +|π|).

We have thus shown that for |π| ≤1,

0≤i≤Nmax Eh

|Y_t^(1,π)_i |²i

+ max

0≤i≤NEh

|f(tj, X_t^(π)_j ,0,0)|²i

≤C. (6)

Analogously to step 1 in Lemma 2.8 we obtain, Eh

|Z_d,t^(1,π)

i |²i2

≤ 1

∆i

Eh

|Y_t^(1,π)_i+1 |²− |Y_t^(1,π)_i |²−2Y_t^(1,π)_i f(ti, X_t^(π)_i ,0,0)∆i

i.

Multiplying with ∆i and summingifrom 0 toN−1 easily gives the L²-bound forZ^(1,π)in view of (6).

As a corollary we obtain a uniform bound for theL²-norms:

Corollary 2.10. Under the assumptions of Theorem 2.3, there is a constant C depending on the data only such that

0≤i≤Nmax Eh

|Y_t^(n,π)_i |²i +

NX−1 i=0

Eh

|Z_t^(n,π)_i |²i

∆i≤C provided|π| is sufficiently small.

Proof. With the notation from the proof of theorem 2.6 we get for sufficiently small|π|,

0≤i≤Nmax Eh

|Y_t^(n,π)_i |²i +

NX−1 i=0

Eh

|Z_t^(n,π)_i |²i

∆i

≤ max

0≤i≤N

Xn ν=1

Ã

E·¯¯¯y_t^(n,π)_i ¯¯¯

2¸ +

N−1X

i=0

E·¯¯¯z_t^(n,π)_i ¯¯¯

2¸

∆i

!

≤ Ã_∞

X

ν=1

A(π, ν)

!2

≤ C

Ã

0≤i≤Nmax Eh

|Y_t^(1,π)_i |²i +

NX−1 i=0

Eh

|Z_t^(1,π)_i |²i

∆i

!

with a constantC depending on the data only. Lemma 2.9 concludes.

3 Generic Analysis of the Error Propagation

For numerical implementation of the iteration proposed in the previous section, one has to approximate the conditional expectations. This section is devoted to an analysis of the error due to the replacement of the conditional expectation by a generic estimator. It turns out that the error grows moderately when the mesh of the partition goes to zero and the number of Picard iterations tends to infinity. We believe, this is an important advantage over the backward scheme, where the error explodes when the mesh tends to zero.

Suppose a generic estimatorEb^π[·|F^t] of the conditional expectation is given.

We consider first the corresponding approximation of the backward scheme of Bouchard and Touzi (2004), namely

Yb_t^(∞,π)_N = ξ^(π), Zb_d,t^(∞,π)_i = Eb^π

·∆Wd,i

∆i

Yb_t^(∞,π)_i+1

¯¯

¯¯Fti

¸ ,

Yb_t^(∞,π)_i = Eb^π[Yb_t^(∞,π)_i+1 |F^ti]−f(ti, X_t^(π)_i ,Yb_t^(∞,π)_i ,Zb_t^(∞,π)_i )∆i. (7) Bouchard and Touzi (2004), Theorem 4.1, prove, under slightly stronger as- sumptions than Assumption 2.1, that

0≤i≤Nmax E[|Yb_t^(∞,π)_i −Y_t^(∞,π)_i |²]

≤ C

|π| max

0≤j≤NE Ã

|Eb^π[Yb_t^(∞,π)_i+1 |F^ti]−E[Yb_t^(∞,π)_i+1 |F^ti]|²

+

¯¯

¯¯Eb^π

·Wt_i+1−Wt_i

ti+1−ti

Yb_t^(∞,π)_i+1

¯¯

¯¯F^ti

¸

−E

·Wt_i+1−Wt_i

ti+1−ti

Yb_t^(∞,π)_i+1

¯¯

¯¯F^ti

¸¯¯¯¯

2!

for some constant Cdepending on the data.

This means, given the same accuracy of the conditional expectation estima- tor the error due to the approximation of the conditional expectation explodes when the mesh of the partition tends to zero. Put differently, due to the nu- merical approximation of the conditional expectation by a Monte-Carlo based estimator one has to simulate the more paths the finer the partition. This increases the computational costs. This effect is particularly unfavorable when the constant in Theorem 2.5 is large (e.g. due to a large Lipschitz constant or time horizon) and, thus, a fine mesh is needed forY_t^(∞,π) to be a good approx- imation of Yt. We note that the described effect has also been observed in the numerical examples by Gobet et al. (2004).

We shall now show that the error due to the approximation of the conditional expectation by its generic estimator does not explode for the discretized Picard iteration. We define

bb^(n,π)_i = ξ^(π)−

N−1X

j=i

f(tj, X_t^(π)_j ,Yb_t^(n−1,π)_j ,Zb_t^(n−1,π)_j )∆j, Yb_t^(n,π)_i = Eb^π[bb^(n,π)_i |F^ti],

Zb_d,t^(n,π)

i = Eb^π

·∆Wd,i

∆i bb^(n,π)_i+1

¯¯

¯¯Fti

¸ ,

initialized at (Yb^(0,π),Zb^(0,π)) = (0,0).

Theorem 3.1. Under Assumption 2.1 there is a constantC depending on the

data such that for any sufficiently fine partition π,

0≤i≤Nmax E[|Yb_t^(n,π)_i −Y_t^(n,π)_i |²] +

N−1X

i=0

E[|Zb_t^(n,π)_i −Z_t^(n,π)_i |²]∆i

≤ C max

1≤ν≤n

Ã

0≤i≤Nmax Eh

|Eb^π[bb^(ν,π)_i |F^ti]−E[bb^(ν,π)_i |F^ti]|²i

+E

NX−1 i=0

¯¯

¯¯Eb^π

·∆Wi

∆i

bb^(ν,π)_i+1

¯¯

¯¯F^ti

¸

−E

·∆Wi

∆i

bb^(ν,π)_i+1

¯¯

¯¯F^ti

¸¯¯¯¯

2

∆i

!

holds for all n∈N. Proof. Define,

b^(n,π)_i = ξ^(π)−

NX−1 j=i

f(tj, X_t^(π)_j , Y_t^(n−1,π)_j , Z_t^(n−1,π)_j )∆j. Then, by Young’s inequality, and with the notation from Lemma 2.9,

0≤i≤Nmax λiE[|Yb_t^(n,π)_i −Y_t^(n,π)_i |²] +

N−1X

i=0

λiE[|Zb_t^(n,π)_i −Z_t^(n,π)_i |²]∆i

≤ 2 Ã

0≤i≤Nmax λiEh

|Eb^π[bb^(n,π)_i |Fti]−E[bb^(n,π)_i |Fti]|²i

+E

N−1X

i=0

λi

¯¯

¯¯Eb

·∆Wi

∆i

bb^(n,π)_i+1

¯¯

¯¯Fti

¸

−E

·∆Wi

∆i

bb^(n,π)_i+1

¯¯

¯¯Fti

¸¯¯¯¯

2

∆i

!

+2 Ã

0≤i≤Nmax λiEh

|E[bb^(n,π)_i −b^(n,π)_i |Fti]|²i

+

N−1X

i=0

λiE

"¯¯¯¯E

·∆Wi

∆i

bb^(n,π)_i+1 −∆Wi

∆i

b^(n,π)_i+1

¯¯

¯¯Fti

¸¯¯¯¯

2#

∆i

! .

Lemma 2.9 can be applied to the second term. Hence, with a suitable choice of Γ and γ,

0≤i≤Nmax λiE[|Yb_t^(n,π)_i −Y_t^(n,π)_i |²] +

NX−1 i=0

λiE[|Zb_t^(n,π)_i −Z_t^(n,π)_i |²]∆i

≤ 2 Ã

0≤i≤Nmax λiEh

|Eb^π[bb^(n,π)_i |Fti]−E[bb^(n,π)_i |Fti]|²i

+E

NX−1 i=0

λi

¯¯

¯¯Eb^π

·∆Wi

∆i

bb^(n,π)_i+1

¯¯

¯¯Fti

¸

−E

·∆Wi

∆i

bb^(n,π)_i+1

¯¯

¯¯Fti

¸¯¯¯¯

2

∆i

!

+ µ1

4+ Γ|π|

¶ µ

0≤i≤Nmax λiE[|Yb_t^(n−1,π)_i −Y_t^(n−1,π)_i |²] +

NX−1 i=0

λiE[|Zb_t^(n−1,π)_i −Z_t^(n−1,π)_i |²]∆i

¶ .

Now for|π|sufficiently small (e.g. less or equal (4Γ)⁻¹) the above estimate can be iterated to obtain the theorem. Note, 1 ≤λi ≤e^ΓT. Thus, we can choose C= 2e^ΓT∨Γ.

4 A Numerical Forward Scheme

In this section we specify an estimator for the conditional expectation. We shall utilize the so-called least-squares Monte-Carlo regression method, which was introduced in Longstaff and Schwartz (2001) in the context of American options and is also applied to the backward scheme in Gobet et al. (2004). The approximation takes place in two steps. First, the conditional expectation is replaced by an orthogonal projection on finite dimensional subspaces. Then, the coefficients of the orthogonal projections are estimated from a sample of independent simulations by the least squares method. Convergence of these two steps will be analyzed in the following subsections. Subsection 4.3 summarizes the results in a Markovian setting relevant for the practical implementation of the numerical scheme.

4.1 Orthogonal Projection on Subspaces of L

( F

)

We will first replace the conditional expectationsE[·|Fti] by orthogonal projec- tions on subspaces ofL²(Fti). Precisely, we fixD+1 subspaces Λd,i, 0≤d≤D, ofL²(Fti) for each 0≤i≤k. The orthogonal projection on Λd,iis denoted by Pd,i.

We now consider the algorithm Yb_t^(n,π)_i = P0,i



ξ^(π)−

N−1X

j=i

f(tj, X_t^(π)_j ,Yb_t^(n−1,π)_j ,Zb_t^(n−1,π)_j )∆j



,

Zb_d,t^(n,π)_i = Pd,i



∆Wd,i

∆i



ξ^(π)−

N−1X

j=i+1

f(tj, X_t^(π)_j ,Yb_t^(n−1,π)_j ,Zb_t^(n−1,π)_j )∆j







,

initialized at (Yb^(0,π),Zb^(0,π)) = 0.

Our aim is to analyze the error of (Yb^(n,π),Zb^(n,π)) as compared to (Y^(n,π), Z^(n,π)) in terms of the projection errors |Y_t^(n,π)_i −P0,i[Y_t^(n,π)_i ]| and

|Z_d,t^(n,π)_i −Pd,i[Z_d,t^(n,π)_i ]|. The main feature of the algorithm – as can be expected in view of Theorem 3.1 – is that the error does not propagate backwards in time. Neither does it explode, when the number of iteration tends to infinity.

This is an important advantage compared to the scheme proposed in Gobet et al. (2004) where the projection errors sum up over the time steps. Roughly speaking, in the Gobet et al. (2004)-scheme the L²-error is bounded by √

N times a constant times the worst L²-projection error (see their Theorem 2).

The following theorem states that in our scheme theL²-error is bounded by a constant times the worstL²-projection error.

Theorem 4.1. Supposef is Lipschitz in(y, z)uniformly in(t, x)with constant

K. Then there is a constantC depending on the data such that

0≤i≤Nmax Eh

|Yb_t^(n,π)_i −Y_t^(n,π)_i |²i +

NX−1 i=0

Eh

|Zb_t^(n,π)_i −Z_t^(n,π)_i |²i

∆i

≤ C

Xn ν=0

µ1 2 +C|π|

¶n−νÃ_N−1 X

i=0

Eh

|Y_t^(ν,π)_i −P0,i[Y_t^(ν,π)_i ]|²i

∆i

+ XD d=1

N−1X

i=0

Eh

|Z_d,t^(ν,π)_i −Pd,i[Z_d,t^(ν,π)_i ]|²i

∆i

!

for sufficiently small|π|. In particular, with a possibly different constant C,

0≤i≤Nmax Eh

|Yb_t^(n,π)_i −Y_t^(n,π)_i |²i +

NX−1 i=0

Eh

|Zb_t^(n,π)_i −Z_t^(n,π)_i |²i

∆i

≤ C max

0≤ν≤n max

0≤i≤N

à Eh

|Y_t^(ν,π)_i −P0,i[Y_t^(ν,π)_i ]|²i

+ XD d=1

Eh

|Z_d,t^(ν,π)

i −Pd,i[Z_d,t^(ν,π)

i ]|²i! .

Proof. We define Y^(n,π)_ti = E



ξ^(π)−

N−1X

j=i

f(tj, X_t^(π)_j ,Yb_t^(n−1,π)_j ,Zb_t^(n−1,π)_j )∆j

¯¯

¯¯

¯¯F^ti



,

Z^(n,π)_d,ti = E



∆Wd,i

∆i



ξ^(π)−

NX−1 j=i

f(tj, X_t^(π)_j ,Yb_t^(n−1,π)_j ,Zb_t^(n−1,π)_j )∆j





¯¯

¯¯

¯¯F^ti



.

Notice, that P0,i

³Y^(n,π)_ti −Y_t^(n,π)_i ´

= Yb_t^(n,π)_i −P0,i

³Y_t^(n,π)_i ´ , Pd,i

³Z^(n,π)_d,ti −Z_d,t^(n,π)_i ´

= Zb_d,t^(n,π)_i −Pd,i

³Z_d,t^(n,π)_i ´ .

Since the orthogonal projection has norm 1 and applying Lemma 2.8 with ˜Y⁽¹⁾= Y^(n,π), ˜Z⁽¹⁾=Z^(n,π), ˜Y⁽²⁾=Y^(n,π), and ˜Z⁽²⁾=Z^(n,π), we obtain:

0≤i≤Nmax λiEh

|Yb_t^(n,π)_i −P0,i(Y_t^(n,π)_i )|²i +

XD d=1

NX−1 i=0

λiEh

|Zb_d,t^(n,π)_i −Pd,i(Z_d,t^(n,π)_i )|²i

∆i

≤ max

0≤i≤NλiEh

|Y^(n,π)_ti −Y_t^(n,π)_i |²i +

NX−1 i=0

λiEh

|Z^(n,π)_ti −Z_t^(n,π)_i |²i

∆i

≤ K²(T+ 1) µ¡

|π|+ Γ⁻¹¢

(γDT+ 1) + D γ

¶

× Ã1

T

N−1X

i=0

λiEh

|Yb_t^(n−1,π)_i −Y_t^(n−1,π)_i |²i +

NX−1 i=0

λiEh

|Zb_t^(n−1,π)_i −Z_t^(n−1,π)_i |²i

∆i

!

for anyγ, Γ>0 withλ0= 1 andλi= (1 + Γ∆i−1)λi−1. The rest of the proof now follows the same lines as the proof of Theorem 3.1 taking into account that, due to the orthogonality of the orthogonal projection,

Eh

|Yb_t^(ν,π)_i −Y_t^(ν,π)_i |²i

=Eh

|Yb_t^(ν,π)_i −P0,i[Y_t^(ν,π)_i ]|²i +Eh

|Y_t^(ν,π)_i −P0,i[Y_t^(ν,π)_i ]|²i .

We also get uniformL²-bounds forYb^(n,π)andZb^(n,π).

Corollary 4.2. Under the assumptions of Theorem 2.3, there is a constant C depending on the data only such that

0≤i≤Nmax Eh

|Yb_t^(n,π)_i |²i +

NX−1 i=0

Eh

|Zb_t^(n,π)_i |²i

∆i≤C provided|π| is sufficiently small.

Proof. This assertion directly follows from Corollary 2.10 and Theorem 4.1, because the orthogonal projection has norm 1.

4.2 A Monte-Carlo Least-Squares Method to Approxi- mate Conditional Expectations

In a next step we replace the projection on subspaces by a simulation based least-squares estimator.

To avoid an overload in notation and since the generalization is straightfor- ward, we shall consider the caseD= 1 only.

We now assume that the projection spaces from the previous section are all finite-dimensional and denote by

{η₁ⁱ, . . . , ηⁱ_K(i)}, resp. {η˜ⁱ₁, . . . ,η˜ⁱ_K(i)˜ }

a basis of Λ0,i and Λ1,i, respectively. The inner-product-matrices associated to these bases are denoted by

Bi=¡

E[η_kⁱηⁱ_l]¢

k,l=0,···K(i), resp. Bei =¡

E[˜ηⁱ_kη˜_lⁱ]¢

k,l=0,···K(i)˜ . In this situation the processes Yb^(n,π)andZb^(n,π) may be rewritten as

Yb_t^(n,π)_i =

K(i)X

k=1

α^(n,π)_i,k η_kⁱ, (8)

Zb_t^(n,π)_i =

K(i)˜

X

k=1

e α^(n,π)_i,k η˜_kⁱ, where, e.g. with ηⁱ= (ηⁱ₁, . . . , η_K(i)ⁱ )^∗,

α^(n,π)_i,· = B⁻¹i E



ηⁱ



ξ^(π)−

N−1X

j=i

f(tj, X_t^(π)_j ,Yb_t^(n−1,π)_j ,Zb_t^(n−1,π)_j )∆j







, (9)

e

α^(n,π)_i,· = Beⁱ⁻¹E



˜ηⁱ∆Wi

∆i



ξ^(π)−

N−1X

j=i+1

f(tj, X_t^(π)_j ,Yb_t^(n−1,π)_j ,Zb_t^(n−1,π)_j )∆j







.