Mean-square convergence of the BDF2-Maruyama and backward Euler schemes for SDE satisfying a global monotonicity condition

# Mean-square convergence of the BDF2-Maruyama and backward Euler schemes for SDE satisfying a global monotonicity condition

Institut für Mathematik, Secr. MA 5-3
Straße des 17. Juni 136
DE-10623 Berlin
Germany
Tel.: +49 (0) 30 314 - 2 96 80
Fax: +49 (0) 30 314 - 2 89 67
22email: andersson@math.tu-berlin.deRaphael Kruse Technische Universität Berlin
Institut für Mathematik, Secr. MA 5-3
Straße des 17. Juni 136
DE-10623 Berlin
Germany
Tel.: +49 (0) 30 314 - 2 33 54
Fax: +49 (0) 30 314 - 2 89 67
44email: kruse@math.tu-berlin.de
Institut für Mathematik, Secr. MA 5-3
Straße des 17. Juni 136
DE-10623 Berlin
Germany
Tel.: +49 (0) 30 314 - 2 96 80
Fax: +49 (0) 30 314 - 2 89 67
22email: andersson@math.tu-berlin.deRaphael Kruse Technische Universität Berlin
Institut für Mathematik, Secr. MA 5-3
Straße des 17. Juni 136
DE-10623 Berlin
Germany
Tel.: +49 (0) 30 314 - 2 33 54
Fax: +49 (0) 30 314 - 2 89 67
44email: kruse@math.tu-berlin.de
###### Abstract

In this paper the numerical approximation of stochastic differential equations satisfying a global monotonicity condition is studied. The strong rate of convergence with respect to the mean square norm is determined to be for the two-step BDF-Maruyama scheme and for the backward Euler-Maruyama method. In particular, this is the first paper which proves a strong convergence rate for a multi-step method applied to equations with possibly superlinearly growing drift and diffusion coefficient functions. We also present numerical experiments for the -volatility model from finance and a two dimensional problem related to Galerkin approximation of SPDE, which verify our results in practice and indicate that the BDF2-Maruyama method offers advantages over Euler-type methods if the stochastic differential equation is stiff or driven by a noise with small intensity.

###### Keywords:
SODE backward Euler-Maruyama method BDF2-Maruyama method strong convergence rates global monotoncity condition
###### Msc:
65C30 65L06 65L20

## 1 Introduction

Strong convergence rates of numerical approximations to stochastic differential equations (SDEs) are a well studied topic. Under a global Lipschitz condition on the coefficients the picture is rather complete, both for one-step methods kloeden1999 (), milstein2004 () and multi-step methods buckwar2006 (), kruse2011 (). Many important equations in application have coefficients that do not satisfy the global Lipschitz condition, and it is therefore important to study a more general setting. Many convergence results for explicit and implicit one-step methods have also been proven for equations without the global Lipschitz condition, see for instance beynisaakkruse2014 (), higham2002b (), hu1996 (), hutzenthaler2014c (), hutzenthaler2014a (), hutzenthaler2012 (), mao2013a (), sabanis2013b (), tretyakov2013 (). In the present paper we determine the strong rate for the backward Euler-Maruyama method (BEM), in the mean-square norm, which improves mao2013a () in terms of a weaker assumption on the coefficients.

For multi-step schemes, on the other hand, there are no previously known results on strong convergence for equations with coefficients not satisfying a global Lipschitz condition. In this paper we determine the strong rate for the BDF2-Maruyama scheme for equations whose, possibly superlinearly growing, coefficient functions satisfy a global monotonicity condition. Backward difference formulas (BDF) are popular in applied sciences for the approximation of stiff equations, see buckwar2006 () for a list of references to such works.

Let , and be a filtered probability space satisfying the usual conditions, on which an -valued standard -Wiener process is defined. We consider the equation

 X(t)=X0+∫t0f(X(s))ds+∫t0g(X(s))dW(s),t∈[0,T], (1)

with drift and diffusion coefficient function . The functions and are assumed to satisfy a global monotonicity, a coercivity and a local Lipschitz condition in Assumption 2.1 below. The initial condition fulfills with some additional integrability, admitting higher moments of the solution.

For a given equidistant time step size we discretize the exact solution to (1) along the temporal grid . Here is uniquely determined by the inequality . We set for . We consider discretizations by means of the backward Euler-Maruyama method

 Xjh−Xj−1h=hf(Xjh)+g(Xj−1h)ΔhWj,j∈{1,…,Nh}, (2)

with satisfying , and by means of the BDF2-Maruyama scheme from buckwar2006 (). The latter is given by the recursion

 32Xjh−2Xj−1h+12Xj−2h=hf(Xjh)+32g(Xj−1h)ΔhWj−12g(Xj−2h)ΔhWj−1, (3)

for , with initial values where is the same as above and is determined, for instance, by one step of the backward Euler scheme or some other one-step method satisfying . In practice, the implementation of the methods (2) and (3) often requires to solve a nonlinear equation in each time step. In Section 3 we discuss that under our assumptions a solution does indeed always exists provided the step size is small enough. The choice of the root-finding algorithm may depend on the coefficient function and its smoothness. We refer to ortega2000 () for a collection of such methods.

We prove that for , determined either by (2) or (3), and being the solution to (1), there exist a constant such that the following mean-square convergence holds:

 maxj∈{0,…,Nh}∥∥X(tj)−Xjh∥∥L2(Ω;Rm)≤C√h,h∈(0,1). (4)

The precise statements of our convergence results are found in Theorems 4.4 and 5.4. The proofs are based on two elementary identities: for all it holds that

 2(u2−u1,u2)=|u2|2−|u1|2+|u2−u1|2, (5)

and for all it holds that

 4(32u3−2u2+12u1,u3)=|u3|2−|u2|2+|2u3−u2|2−|2u2−u1|2+|u3−2u2+u1|2, (6)

found in emmrich2009 (), which has been derived from results on -stability for linear multi-step methods, see girault1979 (); stuart1996 (). Up to the best of our knowledge (6) has not previously been used in the study of the BDF2 scheme for stochastic differential equations.

The paper is organized as follows: Section 2 contains notation and our precise assumptions on the coefficients and in (1). We cite well known results on existence, uniqueness and moment bounds for the solution under these conditions. A well-posedness result for general implicit stochastic difference equations is proved in Section 3. Sections 4 and 5 contain the analysis of the backward Euler-Maruyama and the BDF2-Maruyama schemes, respectively. Subsections 4.1 and 5.1 contain a priori estimates for the respective schemes, in Sections 4.2 and 5.2 stability results are proved, while Subsections 4.3 and 5.3 are concerned with the consistency of the two schemes. The two main results on the strong mean-square convergence rate are stated in Sections 4.4 and 5.4, respectively. Further, in Subsection 5.5 we have a closer look on the second initial value for the BDF2-Maruyama scheme and it is shown that using one step of the BEM method is a feasible choice. Section 6 contains numerical experiments involving the -volatility model from finance which verify our theoretical results and indicate that the BDF2-Maruyama method performs better than Euler-type methods in case of stiff problems or equations with a small noise intensity.

## 2 Setting and preliminaries

### 2.1 Notation and function spaces

Let and denote the scalar product and norm in and let denote the Hilbert-Schmidt norm on the space of all times matrices, i.e., for .

Let be a filtered probability space satisfying the usual conditions. For and a sub--field we denote by the Banach space of all -fold integrable, -measurable random variables taking values in a Banach space with norm

 ∥Z∥Lp(Ω,G,P;E)=(E[|Z|pE])1p,Z∈Lp(Ω,G,P;E).

If we write . If and we obtain the Hilbert space with inner product and norm

 ⟨X,Y⟩=E[(X,Y)],∥X∥=√⟨X,X⟩,

for all . We denote by the norm in , i.e., for .

We next introduce notation related to the numerical discretizations. Recall from Section 1 the temporal grids , . For and , we denote by

 Pjh:L2(Ω,F,P;Rm)→L2(Ω,Ftj,P;Rm),

the orthogonal projector onto the closed sub-space , which is also known as the conditional expectation. More precisely, for we set . We introduce the spaces of all adapted grid functions, which enjoy the following integrability properties

 G2h:={Z :{0,…,Nh}×Ω→Rm:Zn,f(Zn)∈L2(Ω,Ftn,P;Rm),

These will play an important role in the error analysis.

### 2.2 Setting

Consider the setting introduced in Section 1. We now formulate our assumptions on the initial condition and the coefficient functions and which we work with throughout this paper.

###### Assumption 2.1

There exists such that the initial condition satisfies . Moreover, the mappings and , are continuous and there exist and such that for all it holds

 (f(x1)−f(x2),x1−x2)+η∣∣g(x1)−g(x2)∣∣2HS≤L|x1−x2|2, (7) ∣∣f(x1)−f(x2)∣∣≤L(1+|x1|q−1+|x2|q−1)|x1−x2|, (8)

where is the same as above. Further, it holds for all that

 (f(x),x)+4q−32∣∣g(x)∣∣2HS≤L(1+|x|2). (9)

Assumption 2.1 guarantees the existence of an up to modification unique adapted solution to (1) with continuous sample paths, satisfying

 supt∈[0,T]∥X(t)∥L4q−2(Ω;Rm)<∞, (10)

see, e.g., (mao1997, , Chap. 2). In the proof of Theorem 5.3 on the consistency of the BDF2 scheme, the -moment bound is of importance in order to apply the bounds (23), (24) below.

For later reference we note several consequences of Assumption 2.1. ¿From (8) we deduce the following polynomial growth bound:

 ∣∣f(x)∣∣≤~L(1+|x|q),x∈Rm, (11)

where . Indeed, (8) implies that

 ∣∣f(x)∣∣ ≤∣∣f(x)−f(0)∣∣+∣∣f(0)∣∣≤L(1+|x|q−1)|x|+∣∣f(0)∣∣ ≤(2L+|f(0)|)(1+|x|q),x∈Rm.

Moreover, from (7) followed by a use of (8) it holds, for , that

 ∣∣g(x1)−g(x2)∣∣2HS ≤Lη|x1−x2|2+1η∣∣(f(x1)−f(x2),x1−x2)∣∣ ≤Lη(2+|x1|q−1+|x2|q−1)|x1−x2|2.

This gives the local Lipschitz bound

 ∣∣g(x1)−g(x2)∣∣2HS≤2Lη(1+|x1|q−1+|x2|q−1)|x1−x2|2,x1,x2∈Rm, (12)

and, in the same way as above, the polynomial growth bound

 ∣∣g(x)∣∣HS≤¯L(1+|x|q+12),x∈Rm, (13)

where . Finally, we note for later use that the restriction of the exact solution to the time grid , given by

 [X|τh]j:=X(jh),j∈{0,…,Nh},

is an element of the space for every . This follows directly from (10) and the growth bounds (11) and (13).

### 2.3 Preliminaries

Here we list some basic results that we use in this paper. Frequently, we apply the Young inequality and the weighted Young inequality

 ab≤a22+b22andab≤ν2a2+12νb2, (14)

which holds true for all and . We make use of the following discrete version of Gronwall’s Lemma: If , , then

 ∀n∈{1,…,Nh}: an≤c+bhn−1∑j=1ajimplies∀n∈{1,…,Nh}: an≤cebtn. (15)

Finally we cite a standard result from nonlinear analysis which we use for the well-posedness of the numerical schemes, see for instance (ortega2000, , Chap. 6.4) or (stuart1996, , Thm. C.2):

###### Proposition 2.1

Let be a continuous mapping satisfying for some

 (G(x1)−G(x2),x1−x2)≥c|x1−x2|2,x1,x2∈Rm.

Then is a homeomorphism with Lipschitz continuous inverse. In particular, it holds

 ∣∣G−1(y1)−G−1(y2)∣∣≤1c|y1−y2|

for all .

## 3 A well-posedness result for stochastic difference equations

In this section we prove existence and uniqueness of solutions to general stochastic -step difference equations. This result applies in particular to all implicit linear multi-step schemes for SDE with coefficients satisfying Assumption 2.1 including the backward Euler-Maruyama method, the Crank-Nicolson scheme, the -step BDF-schemes, and the -step Adams-Moulton methods. We refer the reader to buckwar2006 (); kruse2011 () for a thorough treatment of these schemes for stochastic differential equations with Lipschitz continuous coefficients.

###### Theorem 3.1

Let the mappings and satisfy Assumption 2.1 with and , let , , , , and with . Assume that initial values are given with and for all . Then, for every there exists a unique family of adapted random variables satisfying

 k∑ℓ=0αk−ℓUj−ℓh=hk∑ℓ=0βk−ℓf(Uj−ℓh)+k∑ℓ=1γk−ℓg(Uj−ℓh)ΔhWj−ℓ+1, (16)

for . In particular, it holds true that and for all .

###### Proof

Let , , be the mappings defined by

 Fh(x):=x−hβkf(x),x∈Rm;h∈(0,h1].

Note that for every it holds that and from the global monotonicity condition (7) we have that

 (Fh(x1)−Fh(x2),x1−x2) =|x1−x2|2−βkh(f(x1)−f(x2),x1−x2) ≥(1−βkh1L)|x1−x2|2.

Consequently, by Proposition 2.1 the inverse of exists for every and is globally Lipschitz continuous with Lipschitz constant . Using these properties and the fact that we can rewrite (16) as

 Fh(Ujh)=Rjh⟺Ujh=F−1h(Rjh) (17)

for all , where

 Rjh:=−k∑ℓ=1αk−ℓUj−ℓh+hk∑ℓ=1βk−ℓf(Uj−ℓh)+k∑ℓ=1γk−ℓg(Uj−ℓh)ΔhWj−ℓ+1,

for . Therefore, by (17) and the continuity of we have that for every , is an adapted collection of random variables, uniquely determined by the initial values .

In order to prove that for all , by means of an induction argument, we introduce . By the assumptions on the initial values it holds that . For the induction step, we now assume that for some . This assumption and the fact that

 ∥∥g(Uj−ℓh)ΔhWj−ℓ+1∥∥2=h|||g(Uj−ℓh)|||2,

imply immediately that . Thus, from the linear growth of we get .

In addition, we recall from (beynisaakkruse2014, , Cor. 4.2) the fact that under Assumption 2.1 the mapping is also globally Lipschitz continuous with a Lipschitz constant independent of . More precisely, there exists a constant such that for all and all it holds true that

 h∣∣g(F−1h(x1))−g(F−1h(x2))∣∣2HS≤C|x1−x2|2. (18)

Consequently, the mapping is also of linear growth and we conclude as above

In particular, this gives . Finally, from the definition of and (17) we have that

 f(Unh)=1βkh(Unh+βkhf(Unh)−Unh)=1βkh(Unh−Fh(Unh))=1βkh(F−1h(Rnh)−Rnh).

Hence, by the linear growth of and the fact that we conclude that .

Altogether, this proves that and, therefore, by induction for every . By finally noting that , , the proof is complete. ∎

## 4 The backward Euler-Maruyama method

In this section we prove that the backward Euler-Maruyama scheme is mean-square convergent of order under Assumption 2.1. The proof is split over several subsections: First we familiarize ourselves with the connection between the BEM method and the identity (5). This is done by proving an a priori estimate in Subsection 4.1. In Subsection 4.2 we then derive a stability result which gives an estimate of the distance between an arbitrary adapted grid function and the one generated by the BEM method. As it turns out this distance is bounded by the error in the initial value and a local truncation error. The latter is estimated for the restriction of the exact solution to (1) to the temporal grid in Subsection 4.3. Altogether, this will then yield the desired convergence result in Section 4.4.

### 4.1 Basic properties of the backward Euler-Maruyama scheme

Here and in Subsection 4.2 we study , , satisfying

 Uj=Uj−1+hf(Uj)+g(Uj−1)ΔhWj,j∈{1,…,Nh}, (19)

with initial condition such that and . Here is the parameter in Assumption 2.1 and from Theorem 3.1 there exist for every a unique satisfying (19). is not necessarily related to the initial value of (1).

In order to prove the a priori bound of Theorem 4.1 and the stability in Theorem 4.2 the following lemma is used:

###### Lemma 1

For all and with satisfying (19) it holds for all -almost surely that

 |Ej|2−|Ej−1|2+|Ej−Ej−1|2 =2h(f(Uj),Ej)+2(g(Uj−1)ΔhWj,Ej−Ej−1)−2(Vj−Vj−1,Ej)+Zj,

where and are the centered random variables given by

 Zj:=2(g(Uj−1)ΔhWj,Ej−1).
###### Proof

From the identity (5), and since satisfies (19) by assumption the assertion follows directly. Note that is well-defined as a centered real-valued integrable random variable due to the independence of the centered Wiener increment and the square integrable random variables and . ∎

The proof of the next theorem is the first and simplest demonstration of the, in principle, same technique used to prove Theorems 4.2, 5.1 and 5.2 below. This a priori estimate is in fact not needed further in the analysis and it can be deduced from the stability Theorem 4.2, but with larger constants, and for a more narrow range for the parameter . We include it for completeness.

###### Theorem 4.1

Let Assumption 2.1 hold with , . For denote by the unique adapted grid function satisfying (19). Then, for all it holds that

where .

###### Proof

Lemma 1 applied with and taking expectations yields

 ∥Uj∥2−∥Uj−1∥2+∥Uj−Uj−1∥2 =2h⟨f(Uj),Uj⟩+2⟨g(Uj−1)ΔhWj,Uj−Uj−1⟩.

From the coercivity condition (9) and the Young inequality (14) we have that

Summing over from to gives that

 ≤2LT+∥U0∥2+h|||g(U0)|||2+(4−4q)hn−1∑j=1|||g(Uj)|||2+2Lhn−1∑j=1∥Uj∥2.

Since it holds and . By elementary bounds we get

We conclude by a use of the discrete Gronwall Lemma (15). ∎

### 4.2 Stability of the backward Euler-Maruyama scheme

For the formulation of the stability Theorem 4.2 we define for and the local truncation error of given by

 ρBEMh(V):=Nh∑j=1∥∥ϱjh(V)∥∥2+1hNh∑j=1∥∥Pj−1hϱjh(V)∥∥2, (20)

where the local residuals of are defined as

 ϱjh(V):=hf(Vj)+g(Vj−1)ΔhWj−Vj+Vj−1,

for . Note that for every . We also introduce a maximal step size for the stability, which guarantees that the stability constant in Theorem 4.2 does not depend on . It is given by

 hE=1max{4L,2}. (21)

Using the same arguments as in Theorem 4.1 the assertion of Theorem 4.2 stays true for all but with a constant depending on as in Theorem 4.1.

###### Theorem 4.2

Let Assumption 2.1 hold with , . For all , satisfying (19), , and all it holds that

 ∥Un−Vn∥2+h|||g(Un)−g(Vn)|||2

where .

###### Proof

Fix arbitrary and . To ease the notation we suppress the dependence of and and simply write, for instance, . We also write and

 Δfj:=f(Uj)−f(Vj),Δgj:=g(Uj)−g(Vj),j∈{0,…,Nh}. (22)

From Lemma 1 we get after taking expectations that

 ∥Ej∥2−∥Ej−1∥2+∥Ej−Ej−1∥2

In order to treat the residual term we first notice that . Then, by taking the adjoint of the projector and by applying the weighted Young inequality (14) with we obtain

 2⟨ϱj,Ej−1⟩ =2⟨Pj−1hϱj,Ej−1⟩≤1h∥∥Pj−1hϱj∥∥2+h∥∥Ej−1∥∥2.

Moreover, further applications of the Cauchy-Schwarz inequality, the triangle inequality, and the weighted Young inequality (14) with yield

 2⟨Δgj−1ΔWj+ϱj,Ej−Ej−1⟩ ≤∥∥Δgj−1ΔWj+ϱj∥∥2+∥∥Ej−Ej−1∥∥2

Therefore, together with the global monotonicity condition (7) this gives

 ∥Ej∥2−∥Ej−1∥2

Setting gives that . Then, summing over from to and thereby identifying two telescoping sums yields

Since as well as and we obtain after some elementary transformations the inequality

 +2(1+2L)hn−1∑j=1∥∥Ej∥∥2.

The proof is completed by applying the discrete Gronwall Lemma (15). ∎

### 4.3 Consistency of the backward Euler-Maruyama scheme

In this subsection we give an estimate for the local truncation error (20) of the BEM method. For the proof we first recall that the restriction of the exact solution to the temporal grid is an element of the space , see Subsection 2.2. Further, we make use of (beynisaakkruse2014, , Lemma 5.5, Lemma 5.6), which provide estimates for the drift integral

 ∫τ2τ1∥∥f(X(τ))−f(X(s))∥∥ds≤C(1+supt∈[0,T]∥∥X(t)∥∥2q−1L4q−2(Ω;Rm))|τ2−τ1|32, (23)

for all with , and for the stochastic integral

 (24)

for all