Stochastic C-Stab. and B-cons. of Milstein-type schemes

# Stochastic C-stability and B-consistency of explicit and implicit Milstein-type schemes

Wolf-Jürgen Beyn Wolf-Jürgen Beyn
Fakultät für Mathematik
Universität Bielefeld
Postfach 100 131
DE-33501 Bielefeld
Germany
Elena Isaak Elena Isaak
Fakultät für Mathematik
Universität Bielefeld
Postfach 100 131
DE-33501 Bielefeld
Germany
and  Raphael Kruse Raphael Kruse
Technische Universität Berlin
Institut für Mathematik, Secr. MA 5-3
Straße des 17. Juni 136
DE-10623 Berlin
Germany
###### Abstract.

This paper focuses on two variants of the Milstein scheme, namely the split-step backward Milstein method and a newly proposed projected Milstein scheme, applied to stochastic differential equations which satisfy a global monotonicity condition. In particular, our assumptions include equations with super-linearly growing drift and diffusion coefficient functions and we show that both schemes are mean-square convergent of order . Our analysis of the error of convergence with respect to the mean-square norm relies on the notion of stochastic C-stability and B-consistency, which was set up and applied to Euler-type schemes in [Beyn, Isaak, Kruse, J. Sci. Comp., 2015]. As a direct consequence we also obtain strong order convergence results for the split-step backward Euler method and the projected Euler-Maruyama scheme in the case of stochastic differential equations with additive noise. Our theoretical results are illustrated in a series of numerical experiments.

###### Key words and phrases:
stochastic differential equations, global monotonicity condition, split-step backward Milstein method, projected Milstein method, mean-square convergence, strong convergence, C-stability, B-consistency
65C30, 65L20

## 1. Introduction

More than four decades ago Grigori N. Milstein proposed a new numerical method for the approximate integration of stochastic ordinary differential equations (SODEs) in [15] (see [16] for an English translation). This scheme is nowadays called the Milstein method and offers a higher order of accuracy than the classical Euler-Maruyama scheme. In fact, G. N. Milstein showed that his method converges with order to the exact solution with respect to the root mean square norm under suitable conditions on the coefficient functions of the SODE while the Euler-Maruyama scheme is only convergent of order , in general.

In its simplest form, that is for scalar stochastic differential equations driven by a scalar Wiener process , the Milstein method is given by the recursion

 (1)

where denotes the step size, is the stochastic increment, and and are the drift and diffusion coefficient functions of the underlying SODE (Equation (3) below shows the SODE in the full generality considered in this paper).

Since the derivation of the Milstein method in [15] relies on an iterated application of the Itō formula, the error analysis requires the boundedness and continuity of the coefficient functions and and their partial derivatives up to the fourth order. Similar conditions also appear in the standard literature on this topic [9, 17, 18].

In more recent publications these conditions have been relaxed: For instance in [10] it is proved that the strong order result for the scheme (1) stays true if the coefficient functions are only two times continuously differentiable with bounded partial derivatives, provided the exact solution has sufficiently high moments and the mapping is globally Lipschitz continuous for every . On the other hand, from the results in [7] it follows that the explicit Euler-Maruyama method is divergent in the strong and weak sense if the coefficient functions are super-linearly growing. Since the same reasoning also applies to the Milstein scheme (1) it is necessary to consider suitable variants in this situation.

One possibility to treat super-linearly growing coefficient functions is proposed in [23]. Here the authors combine the Milstein scheme with the taming strategy from [8]. This allows to prove the strong convergence rate in the case of SODEs whose drift coefficient functions satisfy a one-sided Lipschitz condition. The same approach is used in [12], where the authors consider SODEs driven by Lèvy noise. However, both papers still require that the diffusion coefficient functions are globally Lipschitz continuous.

This is not needed for the implicit variant of the Milstein scheme considered in [6], where the strong convergence result also applies to certain SODEs with super-linearly growing diffusion coefficient functions. However, the authors only consider scalar SODEs and did not determine the order of convergence. The first result bypassing all these restrictions is found in [25], which deals with an explicit first order method based on a variant of the taming idea. A more recent result based on the taming strategy is also given in [13].

In this paper we propose two further variants of the Milstein scheme which apply to multi-dimensional SODEs of the form (3). First, we follow an idea from [2] and study the projected Milstein method which consists of the standard explicit Milstein scheme together with a nonlinear projection onto a sphere whose radius is expanding with a negative power of the step size. The second scheme is a Milstein-type variant of the split-step backward Euler scheme (see [5]) termed split-step backward Milstein method.

For both schemes we prove the optimal strong convergence rate in the following sense: Let and denote the exact solution and its numerical approximation with corresponding step size . Then, there exists a constant independent of such that

 (2) maxn∈{1,…,N}∥X(tn)−Xh(tn)∥L2(Ω;Rd)≤Ch,

where . For the proof we essentially impose the global monotonicity condition (4) and certain local Lipschitz assumptions on the first order derivatives of the drift and diffusion coefficient functions. For a precise statement of all our assumptions and the two convergence results we refer to Assumption 2.1 and Theorems 2.2 and 2.3 below. Together with the result on the balanced scheme found in [25], these theorems are the first results which determine the optimal strong convergence rate for some Milstein-type schemes without any linear growth or global Lipschitz assumption on the diffusion coefficient functions and for multi-dimensional SODEs.

The remainder of this paper is organized as follows: In Section 2 we introduce the projected Milstein method and the split-step backward Milstein scheme in full detail. We state all assumptions and the convergence results, which are then proved in later sections. Further, we apply the convergence results to SODEs with additive noise for which the Milstein-type schemes coincide with the corresponding Euler-type scheme.

The proofs follow the same steps as the error analysis in [2]. In order to keep this paper as self-contained as possible we briefly recall the notions of C-stability and B-consistency and the abstract convergence theorem from [2] in Section 3. Then, in the following four sections we verify that the two considered Milstein-type schemes are indeed stable and consistent in the sense of Section 3. Finally, in Section 8 we report on a couple of numerical experiments which illustrate our theoretical findings. Note that both examples include non-globally Lipschitz continuous coefficient functions, which are not covered by the standard results found in [9, 17].

## 2. Assumptions and main results

This section contains a detailed description of our assumptions on the stochastic differential equation, under which our strong convergence results hold. Further, we introduce the projected Milstein method and the split-step backward Milstein scheme and we state our main results.

Our starting point is the stochastic ordinary differential equation (3) below. We apply the same notation as in [2] and we fix , , and a filtered probability space satisfying the usual conditions. By we denote a solution to the SODE

 (3) dX(t)=f(t,X(t))dt+m∑r=1gr(t,X(t))dWr(t),t∈[0,T],X(0)=X0.

Here stands for the drift coefficient function, while , , are the diffusion coefficient functions. By , , we denote an independent family of real-valued standard -Brownian motions on . For a sufficiently large the initial condition is assumed to be an element of the space .

Let us fix some further notation: We write and for the Euclidean inner product and the Euclidean norm on , respectively. Further, we denote by the set of all bounded linear operators on endowed with the matrix norm induced by the Euclidean norm. For a sufficiently smooth mapping and a given we denote by the Jacobian matrix of the mapping .

Having established this we formulate the conditions on the drift and the diffusion coefficient functions:

###### Assumption 2.1.

The mappings and , , are continuously differentiable. Further, there exist and such that for all and it holds

 (4) ⟨f(t,x1)−f(t,x2),x1−x2⟩+ηm∑r=1∣∣gr(t,x1)−gr(t,x2)∣∣2 ≤L|x1−x2|2.

In addition, there exists such that

 (5) ∣∣∂f∂x(t,x1)−∂f∂x(t,x2)∣∣L(Rd)≤L(1+|x1|+|x2|)q−2|x1−x2|

and, for every ,

 (6) ∣∣∂gr∂x(t,x1)−∂gr∂x(t,x2)∣∣L(Rd) ≤L(1+|x1|+|x2|)q−32|x1−x2|

for all and . Moreover, it holds

 (7) ∣∣∂f∂t(t,x)∣∣≤L(1+|x|)q,∣∣∂gr∂t(t,x)∣∣ ≤L(1+|x|)q+12,

for all , , and all .

First we note that Assumption 2.1 is slightly weaker than the conditions imposed in [25, Lemma 4.2] in terms of smoothness requirements on the coefficient functions. Further, we recall that Equation (4) is often termed global monotonicity condition in the literature. It is easy to check that Assumption 2.1 is satisfied (with ) if and and all their first order partial derivatives are globally Lipschitz continuous. However, Assumption 2.1 includes several SODEs which cannot be treated by the standard results found in [9, 17]. We refer to Section 8 for two more concrete examples.

For a possibly enlarged the following estimates are an immediate consequence of Assumption 2.1 and the mean value theorem: For all and it holds

 (8) |f(t,x)| ≤L(1+|x|)q, (9) ∣∣∂f∂x(t,x)∣∣L(Rd) ≤L(1+|x|)q−1, (10) |f(t1,x)−f(t2,x)| ≤L(1+|x|)q|t1−t2|, (11) |f(t,x1)−f(t,x2)| ≤L(1+|x1|+|x2|)q−1|x1−x2|,

and, for all ,

 (12) |gr(t,x)| ≤L(1+|x|)q+12, (13) ∣∣∂gr∂x(t,x)∣∣L(Rd) ≤L(1+|x|)q−12, (14) |gr(t1,x)−gr(t2,x)| ≤L(1+|x|)q+12|t1−t2|, (15) |gr(t,x1)−gr(t,x2)| ≤L(1+|x1|+|x2|)q−12|x1−x2|.

Thus, Assumption 2.1 implies [2, Assumption 2.1] and all results of that paper also hold true in the situation considered here. Note that in this paper we use the weights instead of as in [2]. For this makes no difference, however in condition (6) we may have if , so that Lipschitz constants actually decrease at infinity.

In the following it will be convenient to introduce the abbreviation

 (16) gr1,r2(t,x):=∂gr1∂x(t,x)gr2(t,x),t∈[0,T],x∈Rd,

for . As above, one easily verifies under Assumption 2.1 that the mappings satisfy (for a possibly larger ) the polynomial growth bound

 (17) ∣∣gr1,r2(t,x)∣∣ ≤L(1+|x|)q

as well as the local Lipschitz bound

 (18) ∣∣gr1,r2(t,x1)−gr1,r2(t,x2)∣∣ ≤L(1+|x1|+|x2|)q−1|x1−x2|

for all , , and . For this conclusion to hold in case , it is essential to use the modified weight function in (6).

We say that an almost surely continuous and -adapted stochastic process is a solution to (3) if it satisfies -almost surely the integral equation

 (19) X(t)=X0+∫t0f(s,X(s))ds+m∑r=1∫t0gr(s,X(s))dWr(s)

for all . It is well-known that Assumption 2.1 is sufficient to ensure the existence of a unique solution to (3), see for instance [11], [14, Chap. 2.3] or the SODE chapter in [20, Chap. 3].

In addition, the exact solution has finite -th moments, that is

 (20) supt∈[0,T]∥∥X(t)∥∥Lp(Ω;Rd)<∞,

if the following global coercivity condition is satisfied: There exist and such that

 (21) ⟨f(t,x),x⟩+p−12m∑r=1∣∣gr(t,x)∣∣2 ≤C(1+|x|2)

for all , . A proof is found, for example, in [14, Chap. 2.4].

For the formulation of the numerical methods we recall the following terminology from [2]: By we denote an upper step size bound. Then, for every we say that is a vector of (deterministic) step sizes if . Every vector of step sizes induces a set of temporal grid points given by

 Th:={tn:=n∑i=1hi:n=0,…,N},

where . For short we write

 |h|:=maxi∈{1,…,N}hi

for the maximal step size in .

Moreover, we recall from [9, 17] the following notation for the stochastic increments: Let with . Then we define

 (22) Is,t(r):=∫tsdWr(τ),

for and, similarly,

 (23) Is,t(r1,r2):=∫ts∫τ1sdWr1(τ2)dWr2(τ1),

where . The joint family of the iterated stochastic integrals is not easily generated on a computer. Besides special cases such as commutative noise one relies on an additional approximation method from e.g. [3, 21, 24]. We also refer to the corresponding discussion in [9, Chap. 10.3].

The first numerical scheme, which we study in this paper, is an explicit one-step scheme and termed projected Milstein method (PMil). It is the Milstein-type counterpart of the projected Euler-Maruyama method form [2] and consists of the standard Milstein scheme and a projection onto a ball in whose radius is expanding with a negative power of the step size.

To be more precise, let , , be an arbitrary vector of step sizes with upper step size bound . For a given parameter the PMil method is determined by the recursion

 (24) ¯¯¯¯¯XPMilh(ti)=min(1,h−αi∣∣XPMilh(ti−1)∣∣−1)XPMilh(ti−1),XPMilh(ti)=¯¯¯¯¯XPMilh(ti)+hif(ti−1,¯¯¯¯¯XPMilh(ti))+m∑r=1gr(ti−1,¯¯¯¯¯XPMilh(ti))Iti−1,ti(r)+m∑r1,r2=1gr1,r2(ti−1,¯¯¯¯¯XPMilh(ti))Iti−1,ti(r2,r1), % for 1≤i≤N,

where . The results of Section 4 indicate that the parameter value for is optimally chosen by setting in dependence of the growth rate appearing in Assumption 2.1. One aim of this paper is the proof of the following strong convergence result for the PMil method. It follows directly from Theorems 4.4 and 5.1 as well as Theorem 3.5.

###### Theorem 2.2.

Let Assumption 2.1 be satisfied with polynomial growth rate . If the exact solution to (3) satisfies , then the projected Milstein method (24) with parameter value and with arbitrary upper step size bound is strongly convergent of order .

Next, we come to the second numerical scheme, which is called split-step backward Milstein method (SSBM). For a suitable upper step size bound and a given vector of step sizes , , this method is defined by setting and by the recursion

 (25) ¯¯¯¯¯XSSBMh(ti)=XSSBMh(ti−1)+hif(ti,¯¯¯¯¯XSSBMh(ti)),XSSBMh(ti)=¯¯¯¯¯XSSBMh(ti)+m∑r=1gr(ti,¯¯¯¯¯XSSBMh(ti))Iti−1,ti(r)+m∑r1,r2=1gr1,r2(ti,¯¯¯¯¯XSSBMh(ti))Iti−1,ti(r2,r1),

for every .

Let us note that the recursion defining the SSBM method evaluates the diffusion coefficient functions at time in the -th step. This phenomenon was already apparent in the definition of the split-step backward Euler method in [2]. It turns out that by this modification we avoid some technical issues in the proofs as condition (26) is applied to and , , simultaneously at the same point in time. Compare also with the inequality (50) further below.

It is shown in Section 6 that the SSBM scheme is a well-defined stochastic one-step method under Assumption 2.1. The second main result of this paper is the proof of the following strong convergence result:

###### Theorem 2.3.

Let Assumption 2.1 be satisfied with and . In addition, we assume that there exist and such that it holds

 (26) ⟨f(t,x1)−f(t,x2),x1−x2⟩+η1m∑r=1∣∣gr(t,x1)−gr(t,x2)∣∣2+η2m∑r1,r2=1∣∣gr1,r2(t,x1)−gr1,r2(t,x2)∣∣2≤L|x1−x2|2

for all . If the solution to (3) satisfies , then the split-step backward Milstein method (25) with arbitrary upper step size bound is strongly convergent of order .

As we show below this theorem follows directly from Theorem 3.5 together with Theorems 6.3 and 7.1. Note that (26) is more restrictive than the global monotonicity condition (4) if the mappings are not globally Lipschitz continuous for all .

In the remainder of this section we briefly summarize the corresponding convergence results in the case of stochastic differential equations with additive noise, that is if the mappings , , do not depend explicitly on the state of . In this case it is well-known that Milstein-type schemes coincide with their Euler-type counterparts.

To be more precise, we consider the solution to an SODE of the form

 (27) dX(t)=f(t,X(t))dt+m∑r=1gr(t)dWr(t),t∈[0,T],X(0)=X0.

In this case, the conditions on the drift coefficient function and the diffusion coefficient functions , , in Assumption 2.1 simplify to

The coefficient functions and , , are continuously differentiable, and there exist constants , such that for all and the following properties hold

 ⟨f(t,x1)−f(t,x2),x1−x2⟩ ≤L|x1−x2|2, ∣∣∂f∂t(t,x)∣∣ ≤L(1+|x|)q, ∣∣∂f∂x(t,x1)−∂f∂x(t,x2)∣∣L(Rd) ≤L(1+|x1|+|x2|)q−2|x1−x2|.

Under this assumption it directly follows that for all for the coefficient functions defined in (16) . Consequently, the PMil method and the SSBM scheme coincide with the PEM method and the SSBE scheme from [2], respectively.

Let us note that Assumption 2.4 implies the global coercivity condition (21) for every . Consequently, under Assumption 2.4 the exact solution to (27) has finite -th moments for every . From this and Theorems 2.2 and 2.3 we directly obtain the following convergence result:

###### Corollary 2.5.

Let Assumption 2.4 be satisfied with and . Then it holds that

1. the projected Euler-Maruyama method with and arbitrary upper step size bound is strongly convergent of order .

2. the split-step backward Euler method with arbitrary upper step size bound is strongly convergent of order .

## 3. A reminder on stochastic C-stability and B-consistency

In this section we give a brief overview of the notions of stochastic C-stability and B-consistency introduced in [2]. We also state the abstract convergence theorem, which, roughly speaking, can be summarized by

 stoch. C-stability + stoch. B-consistency⇒ Strong convergence.

We first recall some additional notation from [2]: For an arbitrary upper step size bound we define the set to be

 T(¯¯¯h):={(t,δ)∈[0,T)×(0,¯¯¯h]:t+δ≤T}.

Further, for a given vector of step sizes , , we denote by the space of all adapted and square integrable grid functions, that is

 G2(Th):={Z:Th×Ω→Rd:Z(tn)∈L2(Ω,Ftn,P;Rd) for all n=0,1,…,N}.

The next definition describes the abstract class of stochastic one-step methods which we consider in this section.

###### Definition 3.1.

Let be an upper step size bound and be a mapping satisfying the following measurability and integrability condition: For every and it holds

 (28) Ψ(Z,t,δ)∈L2(Ω,Ft+δ,P;Rd).

Then, for every vector of step sizes , , we say that a grid function is generated by the stochastic one-step method with initial condition if

 (29) Xh(ti)=Ψ(Xh(ti−1),ti−1,hi),1≤i≤N,Xh(t0)=ξ.

We call the one-step map of the method.

For the formulation of the next definition we denote by the conditional expectation of a random variable with respect to the sigma-field . Note that if is square integrable, then coincides with the orthogonal projection onto the closed subspace . By we denote the associated projector onto the orthogonal complement.

###### Definition 3.2.

A stochastic one-step method is called stochastically C-stable (with respect to the norm in ) if there exist a constant and a parameter value such that for all and all random variables it holds

 (30)

A first consequence of the notion of stochastic C-stability is the following a priori estimate: Let be a stochastically C-stable one-step method. If there exists a constant such that for all it holds

 (31) ≤C0δ, (32) ≤C0δ12,

then there exists a positive constant with

 maxn∈{0,…,N}∥Xh(tn)∥L2(Ω;Rd)≤eCT(∥ξ∥2L2(Ω;Rd)+CC20T)12,

for every vector of step sizes , , where denotes the grid function generated by with step sizes . A proof for this result is found in [2, Cor. 3.6].

###### Definition 3.3.

A stochastic one-step method is called stochastically B-consistent of order to (3) if there exists a constant such that for every it holds

 (33) ∥∥E[X(t+δ)−Ψ(X(t),t,δ)|Ft]∥∥L2(Ω;Rd)≤Cconsδγ+1

and

 (34)

where denotes the exact solution to (3).

Finally, it remains to give our definition of strong convergence.

###### Definition 3.4.

A stochastic one-step method converges strongly with order to the exact solution of (3) if there exists a constant such that for every vector of step sizes , , it holds

 maxn∈{0,…,N}∥∥Xh(tn)−X(tn)∥∥L2(Ω;Rd)≤C|h|γ.

Here denotes the exact solution to (3) and is the grid function generated by with step sizes .

We close this section with the following abstract convergence theorem, which is proved in [2, Theorem 3.7].

###### Theorem 3.5.

Let the stochastic one-step method be stochastically C-stable and stochastically B-consistent of order . If , then there exists a constant depending on , , , , and such that for every vector of step sizes , , it holds

 maxn∈{0,…,N}∥∥X(tn)−Xh(tn)∥∥L2(Ω;Rd)≤C|h|γ,

where denotes the exact solution to (3) and the grid function generated by with step sizes . In particular, is strongly convergent of order .

## 4. C-stability of the projected Milstein method

In this section we prove that the projected Milstein (PMil) method defined in (24) is stochastically C-stable.

Throughout this section we assume that Assumption 2.1 is satisfied with growth rate . First, we choose an arbitrary upper step size bound and a parameter value . Later it will turn out to be optimal to set in dependence of the growth in Assumption 2.1.

For the definition of the one-step map of the PMil method it is convenient to introduce the following short hand notation: For every , we denote the projection of onto the ball of radius by

 (35) x∘:=min(1,δ−α|x|−1)x.

Then, the one-step map is given by

 (36) ΨPMil(x,t,δ):=x∘+δf(t,x∘)+m∑r=1gr(t,x∘)It,t+δ(r)+m∑r1,r2=1gr1,r2(t,x∘)It,t+δ(r2,r1)

for every and . Recall (22) and (23) for the definition of the stochastic increments.

First, we check that the PMil method is a stochastic one-step method in the sense of Definition 3.1. At the same time we verify that the one step map satisfies conditions (31) and (32).

###### Proposition 4.1.

Let the functions and , , satisfy Assumption 2.1 with and let . For every initial value and for every it holds that is a stochastic one-step method.

In addition, there exists a constant only depending on and such that

 (37) ∥∥E[ΨPMil(0,t,δ)|Ft]∥∥L2(Ω;Rd) ≤C0δ, (38) ∥∥(id−E[⋅|Ft])ΨPMil(0,t,δ)∥∥L2(Ω;Rd) ≤C0δ12

for all .

###### Proof.

We first verify that satisfies (28). For this let us fix arbitrary and . Then, the continuity and boundedness of the mapping yields

 Z∘∈L∞(Ω,Ft,P;Rd).

Consequently, by the smoothness of the coefficient functions and by (8), (12), and (17) it follows that

 f(t,Z∘),gr1(t,Z∘),gr1,r2(t,Z∘)∈L∞(Ω,Ft,P;Rd)

for every . Therefore, is an -measurable random variable satisfying condition (28).

It remains to show (37) and (38). From (8) we get immediately that

 ∥∥E[ΨPMil(0,t,δ)|Ft]∥∥L2(Ω;Rd)=∣∣δf(t,0)∣∣≤Lδ.

Next, recall that the stochastic increments and are pairwise uncorrelated. Therefore, we obtain that

 ∥∥(id−E[⋅|Ft])ΨPMil(0,t,δ)∥∥2L2(Ω;Rd) =∥∥m∑r=1gr(t,0)It,t+δ(r)+m∑r1,r2=1gr1,r2(t,0)It,t+δ(r1,r2)∥∥2L2(Ω;Rd) =δm∑r=1