Random walk approximation of BSDEs with H{\"o}lder continuous terminal condition

In this paper we consider the random walk approximation of the solution of a Markovian BSDE whose terminal condition is a locally H{\"o}lder continuous function of the Brownian motion. We state the rate of the L 2-convergence of the approximated solution to the true one. The proof relies in part on growth and smoothness properties of the solution u of the associated PDE. Here we improve existing results by showing some properties of the second derivative of u in space.


Introduction
Let (Ω, F, P) be a complete probability space carrying the standard Brownian motion B = (B t ) t≥0 and assume (F t ) t≥0 is the augmented natural filtration. We consider the following backward stochastic differential equation (BSDE for short) where f is Lipschitz continuous and g is a locally α-Hölder continuous and polynomially bounded function (see (3)). In this paper we are interested in the L 2 -convergence of the numerical approximation of (1) by using a random walk. First results dealing with the numerical approximation of BSDEs date back to the late 1990s. Bally (see [1]) was the first to consider this problem by introducing random discretization, namely the jump times of a Poisson process. In his PhD thesis, Chevance (see [5]) proposed the following discretization of a generator independent from z. The general case (f depends on z, terminal condition ξ ∈ L 2 ) has been studied by Briand, Delyon and Mémin (see [3]). In that paper the authors define an approximated solution (Y n , Z n ) based on random walk and prove weak convergence to (Y, Z) using convergence of filtrations. We also refer to [11], [13], [14], [15] for other numerical methods for BSDEs which use a random walk approach. As in [3], let us introduce the following approximation of B, based on a random walk: where h = T n (n ∈ N * ) and (ε i ) i=1,2,... is a sequence of i.i.d. Rademacher random variables. Consider the following approximated solution (Y n , Z n ) of (Y, Z) The main result of our paper gives the rate of convergence in Theorem 3.1). Basically, we get that the L 2 -norm of the error on Y is of order √ T −v . The proof of this result is based on several ingredients. In particular, we need some estimates on the bound of the first and second derivatives of the solution of the PDE associated to the BSDE (1). We establish these bounds in the case of a forward backward SDE (FBSDE for short) whose terminal condition satisfies the Hölder continuity condition (3). This result extends Zhang [18,Theorem 3.2].
The rest of the paper is organized as follows. Section 2 introduces notations, assumptions and the representation for Z and Z n based on the Malliavin weights. Section 3 states the rate of convergence of the error on Y and Z in L 2 -norm, which is the main result of the paper. Section 4 presents numerical simulations and Section 5 recalls some properties of Malliavin weights, of the regularity of solutions to FBSDEs with a locally Hölder continuous terminal condition function and states some properties of the solutions to the PDEs associated to these FBSDEs.

Preliminaries
This section is dedicated to notations, assumptions and the representation of Z and Z n using the Malliavin weights. Notation: • · p := · L p (P) for p ≥ 1 and for p = 2 simply · . constant. Assumption 2.1.
• g is locally Hölder continuous with order α ∈ (0, 1] and polynomially bounded (p 0 ≥ 0, C g > 0) in the following sense Notice that (3) implies In the rest of the paper, the study of the error (Y n − Y, Z n − Z) will either rely on (2) or on its integral version: where the backward equation (6) arises from (2) by setting Y n r := Y n tm and Z n r := Z n tm for r ∈ [t m , t m+1 ). For n large enough, (6) has a unique solution (Y n , Z n ), and (Y n tm , Z n tm ) n−1 m=0 is adapted to the filtration (G m ) n−1 m=0 . Let us now introduce the Malliavin representations for Z and Z n . They are the cornerstone of our study of the error on Z.

Representations for Z and Z n
We will use the representation (see Ma and Zhang [12,Theorem 4.2]) where E t [·] = E[·|F t ], and for all s ∈ (t, T ] we have Lemma 2.2. Suppose that Assumption 2.1 holds. Then the process Z n given by (6) has the representation Proof. We multiply equation (2) by ε k+1 and take the conditional expectation with respect to G k .
where the l.h.s. is equal to zero. Indeed, for m ≥ k + 1, we have Similarly, for m ≥ k + 1, we get (using [3, Proposition 5.1], where it is stated that both Y n tm and Z n tm can be represented as functions of t m and B n tm ) It remains to divide (9) by √ h and rearrange.

Main result
This section is devoted to the main result of the paper: the rate of the L 2 -convergence of (Y n , Z n ) to (Y, Z). The proof will rely on the fact that the random walk B n can be constructed from the Brownian motion B by Skorohod embedding. Let τ 0 := 0 and define We will use this random walk for our approximation, i.e. we will require Properties satisfied by τ k and B τ k are stated in Lemma A.1. We will denote by E τ k the conditional expectation w.r.t. F τ k .

Remark 3.2. Theorem 3.1 implies that
Proof of Theorem 3.1. Let u : [0, T ) × R → R be the solution of the PDE associated to (1). Since by Theorem 5.4 We first give some properties satisfied by F .

Lemma 3.3. If Assumption 2.1 holds then F is a Lipschitz continuous and polynomially bounded function in x :
where Ψ(x) is given in (5).
Proof of Lemma 3.3. Thanks to the mean value theorem and Theorem 5.4-(ii-c) and (iii-b) we have for The second inequality can be shown similarly.
For the estimate of E|Y t k − Y n t k | 2 we will use (1) and (2): We frequently express conditional expectations with the help of an independent copy of B denoted byB, for example By (3) and Lemma A.1, where To estimate the other term in (11) we consider the decomposition For D 1 we have by Theorem 5.3 that where the last inequality follows from We bound D 2 using Lemma 3.3 and Lemma A.1. Similar to (12) we conclude (setting For D 3 we apply again Lemma 3.3 and Lemma A.1, Finally, using the estimates for the terms D 1 (s, m), D 2 (m), ..., D 4 (m) we arrive at For Z t k − Z n t k we exploit the representations (7) and (8) and estimate Then, similar to (12), we have for the terminal condition by Lemma A.1 that Here we have used thatẼ Then by the conditional Hölder inequality and by (13) as well as by Lemma 3.3 we have Indeed, We estimate T 2 with the help of Lemma 3.3 and Lemma A.1 as follows : Bs−Bt k s−t k ds one notices that by the conditional Hölder inequality, where the last inequality follows in the same way as in (13). Consequently, we have Lemma A.2 enables to bound the second and third term of the r.h.s. by C . Thus we get Then we use (14) and the above estimate to get If this inequality is iterated, one gets a shape where the Gronwall lemma applies. Indeed, setting a m := ( Y tm − Y n tm + Z tm − Z n tm ) one has to consider the double sum Consequently, which gives the bound on the error on Z. Moreover, (14) yields If v ∈ [t k , t k+1 ), we have by Theorem 5.3 that

Numerical simulations
This section deals with the algorithm used to compute (Y n t k , Z n t k ) k=0,··· ,n and numerical experiments for three different terminal conditions. In each case the exact solution is available and we are able to compute the error (Y n − Y, Z n − Z) in L 2 -norm.

4.1
Simulation of (τ 1 , · · · , τ n ) and B n In order to simulate (τ 1 , · · · , τ n ), we use the fact that τ 0 = 0 and ∀k ≥ 1, where (σ k ) 1≤n is an i.i.d. sequence whose common law σ represents the first exit time of the From the book of Borodin and Salminen [2], we have that the Laplace transform of σ is given by E(e −λσ ) = 1 cosh( √ 2λh) . Let F denote the cumulative distribution function of σ. It holds E(e −λσ ) = λF (λ), whereF is the Laplace transform of F . Then, to obtain F , it remains to inverse numerically its Laplace transform. Once we have F , we simulate the sequence (σ k ) 1≤k≤n by following the steps of Algorithm 1.

Simulation of (Y n , Z n )
Since B n is built using the random walk (15), it can be represented by a recombining binomial tree. Both (Y n t k ) 0≤k≤n and (Z n t k ) 0≤k≤n−1 can then also be represented as a recombining binomial tree. Since Y n tn = g(B n tn ), we solve backward in time the BSDE by following these equalities, ensuing from (2) (Y n t k has been replaced by Y n t k+1 in the generator term, but the error induced by this modification is smaller than the ones we consider)

Study of the error E|Y
In this subsection we assume that we are able to compute the exact solution (Y, Z). We want to study numerically the convergence in n of E|Y n t k − Y t k | 2 and E|Z n t k − Z t k | 2 , where (Y, Z) solves (1) and (Y n , Z n ) solves (6). To do so, we approximate the error E|A n t k − A t k | 2 (A = Y or A = Z) by Monte Carlo: 1. For each Monte Carlo simulation, we pick at random one sequence (ξ 1 , · · · , ξ n ) (which gives the value of (B n t 1 , · · · , B n tn )) and one sequence (τ 1 , · · · , τ n ).
2. From the sequence (ξ 1 , · · · , ξ n ) we get the trajectory of Y n , including Y n t k .
3. From the sequence (B τ 1 , · · · , B τn ) (which is equal to (B n t 1 , · · · , B n tn )), we compute B t k by using the Brownian bridge method. We deduce (Y t k , Z t k ) as functions of B t k .
In the following experiments, we plot the logarithm of the errors E Y and E Z (defined in (16)) w.r.t. log(n). From Theorem 3.1, we get that log(E Y ) and log(E Z ) decrease as − α 2 log(n). By using a linear regression, we compute the slope of the line solving the least square problem and compare it to − α 2 .

Case g(x) = e T +x and f (y, z) = y + z
We consider the BSDE with terminal condition g(x) = e T +x and driver f (y, z) = y + z. In this case, we know that Y t = e T +Bt+     Figure 2) represents log(error on Y) (the error is defined by (16)) (resp. log(error on Z)) with respect to log(n). For the Y case, the slope ensuing from the linear regression is −0.53. Even though g(x) = e T +x does not satisfy (3), g is locally Lipschitz continuous, and the outcome seems to be consistent with Theorem 3.1 for α = 1. For the Z case, we get the slope −0.61.

Case g(x)
= x 2 and f (y, z) = y + z In that case, we know that . We run M = 20000 Monte Carlo simulations.   Figure 4 represents log(error on Z) with respect to log(n). The slope of the linear regression is −0.48. The results are then consistent with Theorem 3.1.

Case g(x)
= |x| and f (y, z) = y + z In that case, we know that Y t = e   Here we notice that the modulus of the slope we get is larger than 1 4 , the upper bound obtained in that case in Theorem 3.1.

Some properties of solutions to PDEs and BSDEs
In the following we recall and prove results for FBSDEs with a general forward process, even though we apply them in the present paper only for the case where the forward process is just the Brownian motion. Restricting ourselves to the case of Brownian motion would not shorten the proofs considerably. Let us consider the following SDE started in (t, x), where b and σ satisfy Assumption 5.1.

Malliavin weights
In this section we recall the Malliavin weights and their properties from [9, Subsection 1.1 and Remark 3].

Regularity of solutions to BSDEs
Let us now consider the FBSDE where X t,x is the process satisfying (17). The following result is taken from [9, Theorem 1]. We reformulate it here for the simple situation where we need it. On the other hand, we will use P t,x and are interested in an estimate for all (t, x) ∈ [0, T ) × R. (i) There exists a constant C y 5.3 > 0 such that for 0 ≤ t < s < T and x ∈ R, (ii) There exists a constant C z 5.3 > 0 such that for 0 ≤ t < s < T and x ∈ R, The constants C y 5.3 and C z 5.3 depend on K f , L f , C g , c 1,2 5.4 , T, p 0 , b, σ, κ q and p.
Step 1: We first assume additionally that f : [0, T ] × R 3 → R is continuously differentiable in x, y, and z with uniformly bounded derivatives as it was assumed for [9,Theorem 1]. To take the dependency on x into consideration which arises since we use P t,x , it suffices to replace everywhere in the proof in [9] the constant c B Θ p,∞ by C(T, C g , σ, b, p, p 0 )Ψ(x). The constant C z 5.3 depends moreover on L f and κ q .
Step 2: Now let f be as in Assumption 5.1. In [9, Theorem 1, proof of (C4 l ) =⇒ (C1 l )] a linear BSDE is used which describes the behaviour of the process Z minus its counterpart where the generator is identically 0. Here the partial derivatives of f x , f y , f z appear but only their uniform bound is needed in the estimates. Hence if f satisfies (4), we can use mollifying as explained in (25) below (one may choose N = ∞).
Since |∂ x f ε (t, x, y, z)|, |∂ y f ε (t, x, y, z)| and |∂ z f ε (t, x, y, z)| are bounded by L f we conclude from Step 1 that for all ε > 0 the process Z ε corresponding to f ε satisfies for p ≥ 2. Especially, the family {|Z ε s − Z ε t | q : ε > 0} is then uniformly integrable provided that q < p. By an a priori estimate (cf. [4, Lemma 3.1]) we have that Fubini's theorem implies that there exists a sequence ε m → 0 and a measurable set N ⊆ [0, T ] of Lebesgue measure zero, such that lim m→∞ E|Z r − Z εm r | 2 = 0 for all r ∈ [0, T ] \ N. Consequently, for any q < p and all t, s ∈ [0, T ] \ N with t < s, The assertion follows for all q ≥ 2 since (20) holds for all p ∈ [2, ∞). Since by Theorem 5.4 (ii) the process Z does have a continuous version, we finally get the assertion for all t < s.

Properties of the associated PDE
We collect in the theorem below properties of the solution to the PDE which are mainly known. The new part concerns ∂ 2 x u. For Lipschitz continuous g, the behaviour of ∂ 2 x u has been studied in [19]. General results related to this topic can be found in [7]. and where c 2 5.4 depends on C g , T, p 0 , κ 2 , L f , K f and on the bounds and Lipschitz constants of b and σ.
(iii) u xx exists, and Proof of (ii) (c): We show the assertion for a generator not depending on X, since the terms arising from that dependency would be easy to treat. Since E t,x (E t,x (g(X T ))N t,1 T ) = 0 we can subtract it from the right hand side of (21) and get It holds and thanks to the Cauchy-Schwarz inequality with Ψ 1 = C g (1 + |X T | p 0 + |X t,Xt T | p 0 ) and equation (3), Relation (18) and the Lipschitz continuity of f imply (1 + |u(r, X r )| + |∂ x u(r, X r )σ(r, X r )|)|N t,1 r |dr. Proof. (a) Since g is locally Hölder continuous in the sense of (3), |g(x)| ≤ C g (1 + |x| p 0 +1 ). Then, we get |g ε,N (x)| ≤ C g (1 + (N + 1 + ε) p 0 +1 ), and for f being Lipschitz continuous in y and z, uniformly in time, the same type of result applies.
(b) Since φ is a C ∞ function and f and g are of polynomial growth, we get the result.
(e) We simply have to apply the Lipschitz property of f to get the result.