A stable numerical scheme for stochastic differential equations with multiplicative noise

# A stable numerical scheme for stochastic differential equations with multiplicative noise

C. M. Mora 111Departamento de Ingeniería Matemática, Facultad de Ciencias Físicas y Matemáticas, Universidad de Concepción, Casilla 160 C, Concepción, Chile. e-mails: cmora@ing-mat.udec.cl, hmardones@ing-mat.udec.cl, mselva@ing-mat.udec.cl 444CIMA    H. A. Mardones 111Departamento de Ingeniería Matemática, Facultad de Ciencias Físicas y Matemáticas, Universidad de Concepción, Casilla 160 C, Concepción, Chile. e-mails: cmora@ing-mat.udec.cl, hmardones@ing-mat.udec.cl, mselva@ing-mat.udec.cl 444CIMA    J.C. Jimenez 222Departamento de Matemática Interdisciplinaria, Instituto de Cibernética, Matemática y Física. La Habana, Cuba. e-mail: jcarlos@icimaf.cu    M. Selva 111Departamento de Ingeniería Matemática, Facultad de Ciencias Físicas y Matemáticas, Universidad de Concepción, Casilla 160 C, Concepción, Chile. e-mails: cmora@ing-mat.udec.cl, hmardones@ing-mat.udec.cl, mselva@ing-mat.udec.cl    R. Biscay 333Departamento de Probabilidad y Estadística, Centro de Investigación en Matemática, Guanajuato, México. e-mail: rolando.biscay@cimat.mx
###### Abstract

We introduce a new approach for designing numerical schemes for stochastic differential equations (SDEs). The approach, which we have called direction and norm decomposition method, proposes to approximate the required solution by integrating the system of coupled SDEs that describes the evolution of the norm of and its projection on the unit sphere. This allows us to develop an explicit scheme for stiff SDEs with multiplicative noise that shows a solid performance in various numerical experiments. Under general conditions, the new integrator preserves the almost sure stability of the solutions for any step-size, as well as the property of being distant from . The scheme also has linear rate of weak convergence for a general class of SDEs with locally Lipschitz coefficients, and one-half strong order of convergence.

55footnotetext: CMM, HAM and JCJ was partially supported by FONDECYT Grant 1110787. In addition, CMM and HAM thank the founding of BASAL Grant PFB-03 and CONICYT Grant 21090691, respectively.

Key words. stochastic differential equation, stable numerical scheme, weak error, mean-square convergence, rate of convergence, bilinear SDEs, unstable equilibrium point, locally Lipschitz SDEs.

AMS subject classifications. 60H10, 60H35, 65C20, 65C30, 65C05.

## 1 Introduction

This paper deals with the numerical solution of stiff stochastic differential equations (SDEs) with multiplicative noise. More precisely, we develop an almost sure stable explicit scheme for the Itô SDE

 Xt=X0+∫t0b(Xs)ds+m∑k=1∫t0σk(Xs)dWks, \hb@xt@.01(1.1)

where are independent Wiener processes on a filtered complete probability space , is an adapted -valued stochastic process, and have continuous first-order partial derivatives. In order to do this, we introduce a new method for designing numerical schemes for SDEs with non-constant diffusion coefficients . We are mainly interested in the computation of , where has at most polynomial growth at infinity.

In many cases the stochastic theta-methods, like the backward Euler scheme (LABEL:scheme:BackwardEuler), preserve dynamical properties of (LABEL:eq:1.1) provided that the step-size of the time discretization is small enough (see, e.g., [28, 39, 44, 67, 66]). This does not prevent that the so-called drift-implicit methods (see, e.g., [46, 48]) have poor numerical performance in situations where, for example, some partial derivatives of the diffusion coefficients are not small (see, e.g., [46, 48]). Though a variety of numerical methods for (LABEL:eq:1.1) have been developed in recent times, the schemes for SDEs with multiplicative noise suffer from step-size restrictions due, for instance, to stability issues.

Milstein, Platen and Schurz [46] introduced the general formulation of the balanced implicit methods, a class of fully implicit schemes for (LABEL:eq:1.1) whose implementation depends on the choice of certain weights (see, e.g., [5, 6, 46, 48, 57, 68]), which is a complex problem [58]. The reported balanced schemes exhibit low rate of weak convergence, except incipient progress achieved by [41, 57]. Hutzenthaler, Jentzen and Kloeden [31] designed the tamed Euler scheme for solving (LABEL:eq:1.1) in case satisfies a one-sided Lipschitz condition and are globally Lipschitz continuous (see, e.g., [29, 55] for subsequent developments). Abdulle and Cirilli [2] extended Chebyshev’s methods to solve mean-square stable stiff SDEs (see also, e.g., [3]). Truncating the Brownian motion increments, Milstein, Repin and Tretyakov [47] constructed a class of fully implicit mean-square schemes (see also [34, 47]). Multistep, composite and splitting-step methods have been develop, for instance, in [4, 7, 18, 43, 48, 51]. Using the local linearization method, [14] introduced an exponential scheme for (LABEL:eq:1.1) with . The article [50] develops an integrator of Euler-exponential type for multidimensional SDEs with multiplicative noise (see also, e.g., [22, 35, 37, 59, 60]), and [32] provides a numerical method based on the computation of the conditional mean and the square root of the conditional covariance matrix of a local linearization approximation to (LABEL:eq:1.1). Schemes adapted to specific SDEs are given, for instance, in [10, 13, 15, 50].

To the best of our knowledge, the current numerical methods for (LABEL:eq:1.1), some of them showed up after the first version of this paper [42], have to use small step-sizes in many cases where the diffusion coefficients play an essential role in the dynamics of (see, e.g., Section LABEL:Sec:Numerical_Experiments). This motivates the introduction in Subsection LABEL:subsec:EquilibriumAt0 of a new technique for constructing almost sure stable methods for (LABEL:eq:1.1) with equilibrium point at . In case , we divide the numerical approximation of into the computation of and . Since and satisfy the system (LABEL:eq:2.3) and (LABEL:eq:2.2) given below, we propose to simulate by solving numerically these two coupled SDEs with smooth coefficients. We call this approach Direction and Norm Decomposition method (DND). We take advantage of the unit norm property of and the one-dimensionality of (LABEL:eq:2.3). Namely, this article computes by projecting the Euler-Maruyama scheme applied to (LABEL:eq:2.2) onto the unit sphere, and the norm of is obtained by applying an exponential scheme to the scalar SDE (LABEL:eq:2.3). This yields Scheme LABEL:scheme:EulerStable given below. Subsection LABEL:subsec:GeneralSDEs provides a way to extend Scheme LABEL:scheme:EulerStable to the framework where , may be different from . By adding an auxiliary function, we transform (LABEL:eq:1.1) into a -stochastic differential equation with equilibrium point at . Then, applying Scheme LABEL:scheme:EulerStable, along with suitable approximations, we get Scheme LABEL:scheme:EulerStableG. It is worth pointing out that Scheme LABEL:scheme:EulerStableG becomes Scheme LABEL:scheme:EulerStable whenever . Section LABEL:Sec:Numerical_Experiments presents various numerical experiments that illustrate the very good performance of the new scheme even for large step-sizes.

Suppose for the moment that is an equilibrium point of (LABEL:eq:1.1). In case , have at most linear growth and

 −λ:=supx≠0⎛⎜⎝⟨x,b(x)⟩+12∑mk=1∥∥σk(x)∥∥2∥x∥2−∑mk=1⟨x,σk(x)⟩2∥x∥4⎞⎟⎠<0, \hb@xt@.01(1.2)

Higham, Mao and Yuan [28] prove the almost sure exponential stability of the Euler-Maruyama method for small enough step-sizes (see also [39]). In this situation, (LABEL:eq:1.1) is almost sure exponential stable. If the linear growth condition of is replaced by the one-sided Lipschitz condition and (LABEL:eq:StabCond) is substituted by a slightly stronger requirement, then the backward Euler method

 ¯En+1=¯En+b(¯En+1)Δ+m∑k=1σk(¯En)(WkTn+1−WkTn) \hb@xt@.01(1.3)

is almost sure exponential stable for sufficiently small step-sizes (see [28]). Here and subsequently, for any , where . The almost sure asymptotic stability of the stochastic theta-methods applied to test SDEs has been studied, for instance, in [12, 17, 27, 39, 54, 56, 58]. Similarly, the almost sure asymptotic stability of Balanced schemes has been tested, for example, in [6, 41, 57, 58, 61]. Under (LABEL:eq:StabCond) and the linear growth condition of , we obtain the almost sure exponentially stability of Scheme LABEL:scheme:EulerStable for any step-size (see Theorem LABEL:th:Stable of Subsection LABEL:subsec:LongTime). To the best of our knowledge, Scheme LABEL:scheme:EulerStable is the first numerical method that preserves, in a large class of SDEs, the almost sure asymptotic stability of (LABEL:eq:1.1) no matter the value of the step-size. In case are globally Lipschitz functions, Mao [39] proved that (LABEL:eq:1.1) is small-moment exponentially stable iff the stochastic theta methods applied to (LABEL:eq:1.1) are small-moment exponentially stable for sufficiently small step sizes (see also [28]). We recall that (LABEL:eq:1.1) is small-moment exponentially stable iff (LABEL:eq:1.1) is th moment exponentially stable for a sufficiently small . In Subsection LABEL:subsec:LongTime, we show that Scheme LABEL:scheme:EulerStable preserves the small-moment exponential stability of (LABEL:eq:1.1) for any step-size whenever (LABEL:eq:StabCond) holds (see Theorem LABEL:th:MomentStable).

It is important that the numerical solution of (LABEL:eq:1.1) captures the behavior of when is a non-stable fixed point of (LABEL:eq:1.1). In this direction, we prove that Scheme LABEL:scheme:EulerStable preserves the non-stability of the origin for any step-size provided that , satisfy a general criterion for being a non-stable equilibrium point of (LABEL:eq:1.1) (see Theorem LABEL:th:NoStable below). Previously, [12] has verified that the stochastic theta methods applied to a test SDE reproduce the almost sure instability of when the step-sizes are sufficiently small. On the other hand, Scheme LABEL:scheme:EulerStable also keeps intact the sign of in case , which is an interesting property (see, e.g., [15, 51]).

Many applications deal with the computation of , with . This motivates the study of the weak errors, i.e., the difference between and the expectation of the approximate value of . Using the Kolmorogov equation, Talay [62, 63] and Milstein [45] developed a methodology for obtaining the rate of weak convergence of the numerical schemes for (LABEL:eq:1.1) (see also, e.g., [25, 48, 64] and [20]). Thus, [45, 62, 63] got the linear weak convergence rate of the Euler Maruyama scheme under the global Lipschitz condition. Few weak convergence results are available for SDEs with non-globally Lipschitz continuous coefficients, class of SDEs that appears in important applications. In case are locally Lipschitz functions, Milstein and Tretyakov [49] proposed to discard the numerical trajectories leaving a sufficiently large sphere and studied the weak error involved in this procedure. Hairer, Hutzenthaler and Jentzen [26] showed the existence of a locally Lipschitz SDE with smooth coefficients for which the Euler-Maruyama converges in the weak sense (also in the strong one) without any arbitrarily small polynomial rate of convergence (see also [29, 30]). Bossy and Diop [15] obtained that the symmetrized Euler scheme attains the order of weak convergence when it is applied to (LABEL:eq:1.1) with , and globally Lipschitz continuous, where . In the ergodic case, Talay [66] addressed the computation of integrals with respect to the invariant probability law of a kind of stochastic Hamiltonian system, with additive noise and locally Lipschitz continuous; Talay [66] showed that the discretization error of the backward Euler scheme has the same expansion as in the globally Lipschitz case [67]. We prove that Scheme LABEL:scheme:EulerStableG converges weakly with order under a global coercivity condition and the smoothness of the solution of the Kolmogorov equation associated to (LABEL:eq:1.1) (see Theorem LABEL:th:RateConvergence below), which is a general class of SDEs with multiplicative noise. To this end, we derive a fundamental weak convergence theorem (see Theorem LABEL:th:GeneralWeakConvergence).

Tretyakov and Zhang [68] gave a fundamental mean-square convergence theorem for SDEs with non-global Lipschitz coefficients that satisfy a global monotonicity condition. Moreover, [68] introduced a particular balanced scheme which has rate of strong convergence in a non-global Lipschitz setting (see, e.g., [26, 29, 40, 68] for a recent account of strong convergence results for SDEs with non-global Lipschitz coefficients). We prove that Scheme LABEL:scheme:EulerStableG converge in with rate , where , under the assumptions of the fundamental mean-square convergence theorem proved in [68] (see Theorem LABEL:th:StrongRateConvergence below),

The paper is organized as follows. In Section LABEL:sec:2 we introduce the direction and norm decomposition method (DND). Section LABEL:sec:Convergence is devoted to the stability and convergence properties of the DND scheme. Section LABEL:Sec:Numerical_Experiments provides numerical experiments. All proofs are deferred to Section LABEL:Sec:Proofs, and Section LABEL:sec:Conclusions presents our conclusions.

### 1.1 Notation

For simplicity, we consider the equidistant time discretization , where and We will use the same symbol (resp. and ) for different non-negative real numbers (resp. natural numbers and non-negative increasing functions) that have the common property to be independent of . We set , and denotes the identity operator. For any we define , and we write whenever . Let be the set of all such that for any , with , (i) is continuous; and (ii) for all and . Here and below, and stand for the norm and the dot product (the usual Euclidean scalar product) on , respectively. We say that is in the class if belongs to . The Jacobian matrix of is denoted by .

## 2 Direction and norm decomposition method

Since and are locally Lipschitz, (LABEL:eq:1.1) has a unique continuous strong solution up to an explosion time (see, e.g., [53]), which we assume to be a.s. This happens, for instance, when for all (see, e.g., [24, 38]).

### 2.1 SDEs with equilibrium at 0

Suppose that . Then, there is no loss of generality in assuming a.s., and so, almost surely, for all . In this case, we propose to divide the computation of into the numerical approximations of and the norm of .

We start by obtaining the SDEs describing the evolution of and . Applying Itô’s formula to we obtain

 ∥∥Xt∧˜τj∥∥ =∥X0∥+m∑k=1∫t∧˜τj0⟨Xs,σk(Xs)⟩∥Xs∥dWks +∫t∧˜τj0⎛⎜⎝⟨Xs,b(Xs)⟩+12∑mk=1∥∥σk(Xs)∥∥2∥Xs∥−12m∑k=1⟨Xs,σk(Xs)⟩2∥Xs∥3⎞⎟⎠ds,

where . Since for all , , and so taking limit as gives

 ∥Xt∥ =∥X0∥+m∑k=1∫t0⟨Xs,σk(Xs)⟩∥Xs∥dWks \hb@xt@.01(2.1) +∫t0⎛⎜⎝⟨Xs,b(Xs)⟩+12∑mk=1∥∥σk(Xs)∥∥2∥Xs∥−12m∑k=1⟨Xs,σk(Xs)⟩2∥Xs∥3⎞⎟⎠ds.

We rewrite (LABEL:eq:2.1) as

 ∥Xt∥=∥X0∥+∫t0β(∥Xs∥,Zs)∥Xs∥ds+m∑k=1∫t0⟨Zs,¯σk(∥Xs∥,Zs)⟩∥Xs∥dWks, \hb@xt@.01(2.2)

where for any and we define

 \hb@xt@.01(2.3)

Applying Itô’s formula to we get after a long calculation that

 Xt∧˜τj∥∥Xt∧˜τj∥∥ =X0∥X0∥+∫t∧˜τj0(b(Xs)∥Xs∥−⟨Xs∥Xs∥,b(Xs)∥Xs∥⟩Xs∥Xs∥)ds +12m∑k=1∫t∧˜τj0⎛⎝3⟨Xs∥Xs∥,σk(Xs)∥Xs∥⟩2−⟨σk(Xs)∥Xs∥,σk(Xs)∥Xs∥⟩⎞⎠Xs∥Xs∥ds +m∑k=1∫t∧˜τj0(σk(Xs)∥Xs∥−⟨Xs∥Xs∥,σk(Xs)∥Xs∥⟩Xs∥Xs∥)dWks.

This implies

 Zt =Z0+∫t0(¯b(∥Xs∥,Zs)−⟨Zs,¯b(∥Xs∥,Zs)⟩Zs+Ψ(∥Xs∥,Zs))ds \hb@xt@.01(2.4)

where was defined by and

 \hb@xt@.01(2.5)

Since , the main idea of this paper is to compute by solving the system of SDEs formed by (LABEL:eq:2.3) and (LABEL:eq:2.2), which defines the Direction and Norm Decomposition method (DND). Taking advantage of the unit norm property of and the one-dimensionality of (LABEL:eq:2.3), we next design a simple numerical scheme based on our DND method.

Suppose that the -measurable random variable simulates the initial condition . Then, we set and . In order to compute and , we next generate recursively the pairs of -measurable random variables such that takes values in and lies on the unit sphere of for any . Here, and will approximate and , respectively. Fix and , which satisfy and . From (LABEL:eq:2.1) it follows that for all ,

Freezing , over at the values , , and replacing the pair , by , , we obtain the linear scalar SDE

 ηt=¯ηn+∫tTnβ(¯ηn,^Zn)ηsds+m∑k=1∫tTn⟨^Zn,¯σk(¯ηn,^Zn)⟩ηsdWks,

whose solution at time is

 ηTn+1=¯ηnexp((β(¯ηn,^Zn)−12m∑k=1⟨^Zn,¯σk(¯ηn,^Zn)⟩2)Δ \hb@xt@.01(2.6) +m∑k=1⟨^Zn,¯σk(¯ηn,^Zn)⟩(WkTn+1−WkTn)).

We can solve (LABEL:eq:2.2) using various iterative schemes. For simplicity, we apply the Euler approximation to

 Zt +m∑k=1∫tTn(¯σk(∥Xs∥,Zs)−⟨Zs,¯σk(∥Xs∥,Zs)⟩Zs)dWks,

which yields

 ZTn+1 ≈ZTn+(¯b(∥XTn∥,ZTn)−⟨ZTn,¯b(∥XTn∥,ZTn)⟩ZTn+Ψ(∥XTn∥,ZTn))Δ +m∑k=1(¯σk(∥XTn∥,ZTn)−⟨ZTn,¯σk(∥XTn∥,ZTn)⟩ZTn)(WkTn+1−WkTn).

Hence, substituting and by and , respectively, we obtain

 ˜Zn+1 \hb@xt@.01(2.7) +m∑k=1(¯σk(¯ηn,^Zn)−⟨^Zn,¯σk(¯ηn,^Zn)⟩^Zn)(WkTn+1−WkTn).

As we are interested in simulating the distribution of , in (LABEL:eq:2.22) and (LABEL:eq:2.23) we replace by , where are independent and identically distributed (i.i.d.) -measurable random variables with symmetric law and variance , which are also independent of . Then and become

 ¯ηn+1:=¯ηnexp(⟨^Zn,¯b(¯ηn,^Zn)⟩Δ+m∑k=1⟨^Zn,¯σk(¯ηn,^Zn)⟩√Δ^Wkn+1

and

 ¯Zn+1 +m∑k=1(¯σk(¯ηn,^Zn)−⟨^Zn,¯σk(¯ηn,^Zn)⟩^Zn)√Δ^Wkn+1.

Finally, we take , and so we have defined the new pair . The fact that allows us to improve the accuracy of by projecting on the unit sphere, a normalization procedure that has been used with success in the numerical solution of the non-linear Schödinger equations (see, e.g., [50, 52]) and the computation of Lyapunov exponents (see, e.g., [19, 65]). Since approximates the law of , we can expect that is not close to , and so the weak approximation of should reproduce efficiently the unit-norm property of without incurring round-off errors. In summary, we have designed the following numerical scheme.

###### Scheme 1

Let and be random variables satisfying and . Consider the i.i.d. symmetric random variables with variance that are independent of . For any , we define recursively the pair by

 \hb@xt@.01(2.8)

and where

 ¯Zn+1 \hb@xt@.01(2.9) +m∑k=1(¯σk(¯ηn,^Zn)−⟨^Zn,¯σk(¯ηn,^Zn)⟩^Zn)√Δ^Wkn+1,

is given by (LABEL:eq:2.21), the functions , are described by (LABEL:eq:3.5), and

 μn=⟨^Zn,¯b(¯ηn,^Zn)⟩+12m∑k=1∥∥ ¯σk(¯ηn,^Zn)∥∥2−m∑k=1⟨^Zn,¯σk(¯ηn,^Zn)⟩2.

Here, approximates the solution of the SDE (LABEL:eq:1.1) with for all .

###### Remark 2.1

Since could be approximately , we implement Scheme LABEL:scheme:EulerStable by computing and rather than , which avoids possible round-off errors (see Subsection LABEL:subsec:R_Errors).

If (LABEL:eq:1.1) reduces to the bilinear SDE

 Xt=X0+∫t0BXs ds+m∑k=1∫t0σkXs dWks \hb@xt@.01(2.10)

with , then Scheme LABEL:scheme:EulerStable becomes

###### Scheme 2

Define recursively where

 ¯Zn+1=^Zn+Bn^ZnΔ+m∑k=1(σk−⟨^Zn,σk^Zn⟩)^Zn√Δ^Wkn+1 \hb@xt@.01(2.11)

with as in Scheme LABEL:scheme:EulerStable, and

 Bn=B−⟨^Zn,B^Zn⟩+m∑k=1(32⟨^Zn,σk^Zn⟩2−⟨^Zn,σk^Zn⟩σk−12∥∥σk^Zn∥∥2).

The stochastic process is given by the iterative formula

 ¯ηn+1=¯ηnexp((⟨^Zn,B^Zn⟩+12m∑k=1∥∥σk^Zn∥∥2−m∑k=1⟨^Zn,σk^Zn⟩2)Δ+m∑k=1⟨^Zn,σk^Zn⟩√Δ^Wkn+1).

Thus, approximates the solution of (LABEL:eq:BilinealSDE) for all .

In case , from (LABEL:eq:2.4) we have for all . Therefore, Scheme LABEL:scheme:EulerStable reduces to the follow

###### Scheme 3

Given ,

 ¯Xn+1=¯Xnexp((¯b(¯Xn)−12m∑k=1¯σk(¯Xn)2)Δ+m∑k=1¯σk(¯Xn)√Δ^Wkn+1)

approximates the solution of (LABEL:eq:1.1) with , where and the ’s are i.i.d. symmetric random variables with variance that are independent of .

###### Remark 2.2

A projection technique on the sphere was introduced in [65] to approximate the upper Lyapunov exponent of (LABEL:eq:1.1) with , linear. As a difference with the DND approach considered here, the method of [65] does not involve the numerical solution of the coupled system of the SDEs (LABEL:eq:2.3) and (LABEL:eq:2.2), which describes the evolution of the norm of and the projection of on the unit sphere, respectively.

### 2.2 General SDEs

Suppose that at least one of the vectors , is different from . Consider the constant function , where . From (LABEL:eq:1.1) it follows that

 \hb@xt@.01(2.12)

with and for all and . As , we can compute by applying Scheme LABEL:scheme:EulerStable to (LABEL:eq:3.4). Since , a better alternative is to approximate and in by the solution of

 (˜Xt˜Vt)=(¯Xnα)+∫tTnf((˜Xs˜Vs))ds+m∑k=1∫tTngk((˜Xs˜Vs))dWks, \hb@xt@.01(2.13)

with . Then , where and are given by one iteration of Scheme LABEL:scheme:EulerStable applied to (LABEL:eq:3.1). Hence, we compute by projecting onto its first coordinates. After a short algebraic manipulation we get Scheme LABEL:scheme:EulerStableG with (see Remark LABEL:rem:ObtencionEulerStableG for details).

###### Scheme 4

Let be a random variable with values in such that in case . Suppose that are i.i.d. symmetric random variables with variance that are independent of . Choose such that if , and otherwise. Then, for all the solution of (LABEL:eq:1.1) is recursively approximated by

 ¯Xn+1=⎧⎨⎩¯ρn+1¯Un+1/√∥∥¯Un+1∥∥2+(¯Vn+1)2if (¯Un+1,¯Vn+1)≠0,¯ρn+1^Unif (¯U