Adaptive Euler-Maruyama method for SDEs with non-globally Lipschitz drift: Part II, infinite time interval

TL;DR

This paper introduces an adaptive Euler-Maruyama method for ergodic SDEs with non-globally Lipschitz drift over infinite time, ensuring stability, uniform moment bounds, and strong convergence, enabling efficient multilevel Monte Carlo simulations.

Contribution

It develops an adaptive timestep scheme for ergodic SDEs with non-globally Lipschitz drift, providing uniform bounds and convergence results over infinite time intervals.

Findings

Numerical experiments confirm theoretical stability and convergence.

Adaptive method achieves uniform moment bounds over infinite time.

Supports efficient multilevel Monte Carlo simulations.

Abstract

This paper proposes an adaptive timestep construction for an Euler-Maruyama approximation of the ergodic SDEs with a drift which is not globally Lipschitz over an infinite time interval. If the timestep is bounded appropriately, we show not only the stability of the numerical solution and the standard strong convergence order, but also that the bound for moments and strong error of the numerical solution are uniform in T, which allow us to introduce the adaptive multilevel Monte Carlo. Numerical experiments support our analysis.

Equations

$$\mathrm{d}{X}_t=f(X_t)\mathrm{d}t+g(X_t)\mathrm{d}{W}_t,$$

$$⟨x,f(x)⟩\le -\alpha \Vert x{\Vert }^2+\beta ,$$

$$\pi (\phi ):=\int \phi (x)\mathrm{d}\pi (x)=\underset{t\to \infty }{\lim }\mathbb{E}\left[\phi (X_t)\right],\phi \in L^1(\pi ).$$

$$\underset{T\to \infty }{\lim }\frac{1}{T}{\int }_0^T\phi (X_t)\mathrm{d}t=\pi (\phi ),a.s.,$$

$$\frac{1}{N}\sum _{n=1}^N\phi ({\hat{X}}_{nh}),$$

$$\underset{N\to \infty }{\lim }\frac{1}{N}\sum _{n=1}^N\phi ({\hat{X}}_{nh})=\int \phi (x)\mathrm{d}\pi (x)+O(h).$$

$${\hat{X}}_{(n+1)h}={\hat{X}}_{nh}+f({\hat{X}}_{nh})h+\Delta W_n$$

$${\hat{X}}_{nh}^{\ast }$$

$${\hat{X}}_{(n+1)h}$$

$${\hat{X}}_{(n+1)h}={\hat{X}}_{nh}+f({\hat{X}}_{(n+1)h})h+g({\hat{X}}_{nh})\Delta W_n.$$

$$\mathrm{d}{\tilde{X}}_t=(f(x)+\nabla f(x)({\tilde{X}}_t-x))\mathrm{d}t+\mathrm{d}{W}_t,$$

$${\hat{X}}_{t_n}^{\ast }$$

$${\hat{X}}_{t_{n+1}}$$

$$h_n\le \min (2h_{n-1},h_{\max }),k_n=\min (k\in \mathbb{Z}:|f({\hat{X}}_{t_n}^{\ast })-f({\hat{X}}_{t_n})|\le \tau ).$$

$$|\mathbb{E}[\phi (X_T)]-\pi (\phi )|\le \epsilon .$$

$$t_{n+1}=t_n+h_n,{\hat{X}}_{t_{n+1}}={\hat{X}}_{t_n}+f({\hat{X}}_{t_n})h_n+g({\hat{X}}_{t_n}){\Delta W}_n,$$

$${\hat{X}}_t={\hat{X}}_{\underline{t}}+f({\hat{X}}_{\underline{t}})(t-\underline{t})+g({\hat{X}}_{\underline{t}})(W_t-W_{\underline{t}}),$$

$$\mathrm{d}{\hat{X}}_t=f({\hat{X}}_{\underline{t}})\mathrm{d}t+g({\hat{X}}_{\underline{t}})\mathrm{d}{W}_t=f({\overline{X}}_t)\mathrm{d}t+g({\overline{X}}_t)\mathrm{d}{W}_t.$$

$$\Vert f(x)-f(y)\Vert +\Vert g(x)-g(y)\Vert \le C_R\Vert x-y\Vert$$

$$⟨x,f(x)⟩\le -\alpha \Vert x{\Vert }^2+\beta ,$$

$$\Vert g(x){\Vert }^2\le \beta .$$

$$\mathbb{E}\left[\Vert X_t{\Vert }^p\right]\le C_p.$$

$${\mathrm{e}}^{p\alpha t/2}\Vert X_t{\Vert }^p-\Vert X_0{\Vert }^p$$

$$\mathbb{E}\left[{\mathrm{e}}^{p\alpha t/2}\Vert X_t{\Vert }^p\right]-\mathbb{E}\left[\Vert X_0{\Vert }^p\right]$$

$$\frac{p(p+1)\beta }{2}{\mathrm{e}}^{p\alpha s/2}\Vert X_s{\Vert }^{p-2}\le \frac{p\alpha }{2}{\mathrm{e}}^{p\alpha s/2}\Vert X_s{\Vert }^p+c_p{\mathrm{e}}^{p\alpha s/2}$$

$$\mathbb{E}\left[{\mathrm{e}}^{p\alpha t/2}\Vert X_t{\Vert }^p\right]-\mathbb{E}\left[\Vert X_0{\Vert }^p\right]$$

$$\mathbb{E}\left[\Vert X_t{\Vert }^p\right]\le \frac{2c_p}{p\alpha }+{\mathrm{e}}^{-p\alpha t/2}\mathbb{E}\left[\Vert X_0{\Vert }^p\right]\le \frac{2c_p}{p\alpha }+\Vert x_0{\Vert }^p\triangleq C_p.$$

$$⟨x,f(x)⟩+\frac{1}{2}h(x)\Vert f(x){\Vert }^2\le -\alpha \Vert x{\Vert }^2+\beta .$$

$$\mathbb{E}\left[\Vert {\hat{X}}_t{\Vert }^p\right]<C_p,\mathbb{E}\left[\Vert {\overline{X}}_t{\Vert }^p\right]<C_p.$$

$$h(x)\ge {\left(\xi \Vert x{\Vert }^q+\zeta \right)}^{-1}.$$

$$\mathbb{E}\left[(N_T-1{)}^p\right]<C_{h,p}{T}^p.$$

Full text

60H10,

60H35,

65C30,

SDE,

Euler-Maruyama,

strong convergence,

adaptive timestep,

arXiv:0000.0000

\startlocaldefs

\endlocaldefs

and

1 Introduction

In this paper we consider an $m$-dimensional stochastic differential equation (SDE) driven by a $d$-dimensional Brownian motion:

[TABLE]

which has a non-globally Lipschitz drift $f:{\mathbb{R}}^{m}\!\rightarrow\!{\mathbb{R}}^{m}$ satisfying the dissipativity condition: for some $\alpha,\beta>0,$

[TABLE]

and a bounded volatility $g:{\mathbb{R}}^{m}\!\rightarrow\!{\mathbb{R}}^{m\times d}.$ In particular, we focus on a class of SDEs which are ergodic in nature and converge exponentially to some invariant measure $\pi$. Evaluating the expectation of some function $\varphi(x)$ with respect to that invariant measure $\pi$ is of great interest in mathematical biology, physics and Bayesian inference in statistics:

[TABLE]

Several different methodologies have been developed to estimate it.

First, we can compute the probability density function $\rho(x)$ of $\pi$ by solving the corresponding stationary Fokker-Planck equation, see [14]. However, the stationary Fokker-Planck equation is a partial differential equation (PDE) and its numerical solution becomes extremely expensive when the dimension of the PDEs becomes large.

The second approach is based on the ergodicity of the SDEs:

[TABLE]

where the limit does not depend on initial value $x_{0}.$ This approach uses discretized numerical schemes to approximate the SDEs and requires the numerical solution $\widehat{X}_{t}$ to preserve the ergodicity. In practice, we can choose a sufficiently large $N$ and compute

[TABLE]

where ${\widehat{X}}_{nh}$ is the numerical solution at the $n$th discretized time point using an ergodic method with a uniform timestep $h.$ Under the dissipativity condition (2) together with the Lipschitz condition for $f$, Talay [15] shows the standard weak convergence order for the Milstein method:

[TABLE]

Roberts $\&$ Tweedie in [13] analyse the ergodicity of the unadjusted Langevin algorithm for the Langevin equation, which has uniform volatility and satisfies dissipativity condition (2). This scheme corresponds to the standard Euler-Maruyama method:

[TABLE]

using a uniform timestep of size $h$ with Brownian increments $\Delta W_{n}.$ The paper shows that the numerical solution is not ergodic when $f$ has a polynomial degree larger than 1. Metropolis-adjusted Langevin algorithm (MALA) is introduced but the numerical solutions are still not exponentially ergodic for non-linear drift $f.$

For globally Lipschitz Langevin SDEs satisfying the dissipativity condition (2), the standard Euler-Maruyama method is shown in [10] to inherit ergodicity provided the timesteps are sufficiently small. However, the standard Euler-Maruyama method and Milstein method fail to be stable for non-globally Lipschitz SDEs. The Split-step backward Euler method

[TABLE]

and the drift-implicit Backward Euler method:

[TABLE]

are proved to be ergodic in [10].

Under the same conditions, Hansen in [6] considers the local linearization of the drift coefficient, that is the first-order Taylor approximation of $\tilde{X}_{t}$: for $t\in[nh,(n+1)h],$ given $\tilde{X}_{nh}=x,$

[TABLE]

which is an Ornstein-Uhlenbeck diffusion and we can calculate the analytical conditional distribution of $\tilde{X}_{(n+1)h}.$ However, this analytical treatment only applies when the diffusion coefficient is uniform.

An adaptive timestepping algorithm proposed by Lamba, Mattingly $\&$ Stuart in [8] chooses the step size by halving or doubling based on the local error estimation and a user-input tolerance $\tau$. More precisely,

[TABLE]

where $h_{n}=2^{-k_{n}}{h_{\rm max}}$ satisfies that

[TABLE]

This scheme is proved to preserve the ergodicity of original SDEs under the dissipativity condition (2) and boundedness and invertibility of the diffusion coefficient $g.$

Lemaire [9] considers an infinite time interval under the dissipativity condition generated by a general Lyapunov function using a timestep with an upper bound which decreases towards zero over time, and proves convergence of the empirical distribution to the invariant distribution of the SDE.

Finally, without requiring the ergodicity of the schemes, for exponentially ergodic SDEs, we can choose a sufficiently large $T$ such that

[TABLE]

Then, for this fixed $T,$ we can use of all the methods mentioned in Part I of this pair of articles [1] to estimate $\mathbb{E}[\varphi(X_{T})].$ Milstein $\&$ Tretyakov [12] analyse the error of this kind of approach based on their quasi-symplectic method. In practice, a suitable choice of initial data is important because the transition period to a sufficient proximity of the equilibrium can be rather long. Therefore, running a small number of pioneer paths can be employed to obtain a good initial distribution for the overwhelming majority of simulations.

We should remark here that the PDE approach is far too expensive in high dimensions. The time-averaging approach (3) requires the numerical methods to preserve the ergodicity but the third approach does not. The length of the time interval $[0,T]$ used in the time-averaging approach is much longer than the third approach, for it needs not only to ensure that the distribution of ${\widehat{X}}_{T}$ is sufficiently close to the invariant measure $\pi,$ but also to guarantee a small variance for the average. Therefore, the second approach needs to simulate a single long path but the third one needs to simulate a lot of relative short paths, which allows multilevel Monte Carlo (MLMC) and parallel computing techniques to be employed. One important concern about the third approach is that in the numerical analysis of existing algorithms in the finite time interval $[0,T],$ the strong error increases exponentially as $T$ increases. This is not acceptable when we need to simulate a much larger $T$ and it will be a key concern in this paper. The final issue about the second and third approaches is how to choose a good $T$ such that the weak error is bounded appropriately.

In this paper, we propose instead to use the standard explicit Euler-Maruyama method, but with an adaptive timestep $h_{n}$ which is a function of the current approximate solution ${\widehat{X}}_{t_{n}}$. By setting a suitable condition for $h,$ we can show that, instead of an exponential bound, the numerical solution has a uniform bound with respect to $T$ for both moments and the strong error. Then, MLMC methodology [2, 3] is employed and non-nested timestepping is used to construct an adaptive MLMC [4]. Following the idea of Glynn and Rhee [5] to estimate the invariant measure of some Markov chains, we introduce an adaptive MLMC algorithm for the infinite interval, in which each level $\ell$ has a different time interval length $T_{\ell},$ to achieve a better computational performance. Note that using different time interval lengths allows us not to worry how to choose an appropriate $T$ before simulation. The MLMC algorithm will automatically terminate at a level $L$ with a sufficiently large $T_{L}.$

The rest of the paper is organised as follows. Section 2 states the main theorems and proves some minor lemmas. Section 3 introduces the MLMC schemes, and the relevant numerical experiments are provided in section 4. The proofs of the three main theorems are deferred to section 5, and finally, section 6 has some conclusions and discusses future extensions.

In this paper we consider the infinite time interval $[0,\infty)$ and let $(\Omega,\mathcal{F},{\mathbb{P}})$ be a probability space with normal filtration $(\mathcal{F}_{t})_{t\in[0,\infty)}$ for section 2 and $(\mathcal{F}_{t})_{t\in(-\infty,0]}$ for section 3 corresponding to a $d$-dimensional standard Brownian motion $W_{t}=(W^{(1)},W^{(2)},\ldots,W^{(d)})_{t}.$ We denote the vector norm by $\|v\|\triangleq(|v_{1}|^{2}+|v_{2}|^{2}+\ldots+|v_{m}|^{2})^{\frac{1}{2}}$, the inner product of vectors $v$ and $w$ by $\langle v,w\rangle\triangleq v_{1}w_{1}+v_{2}w_{2}+\ldots+v_{m}w_{m}$, for any $v,w\in{\mathbb{R}}^{m}$ and the Frobenius matrix norm by $\|A\|\triangleq\sqrt{\sum_{i,j}A_{i,j}^{2}}$ for any $A\in{\mathbb{R}}^{m\times d}.$

2 Adaptive algorithm and theoretical results

2.1 Adaptive Euler-Maruyama method

The adaptive Euler-Maruyama discretisation is

[TABLE]

where $h_{n}\triangleq h({\widehat{X}}_{t_{n}})$ and ${\Delta W}_{n}\triangleq W_{t_{n+1}}\!-\!W_{t_{n}}$, and there is fixed initial data $t_{0}\!=\!0,\ {\widehat{X}}_{0}\!=\!x_{0}$.

We use the notation ${\underline{t}}\triangleq\max\{t_{n}:t_{n}\!\leq\!t\},\ n_{t}\triangleq\max\{n:t_{n}\!\leq\!t\}$ for the nearest time point before time $t$, and its index.

We define the piecewise constant interpolant process ${\overline{X}}_{t}={\widehat{X}}_{\underline{t}}$ and also define the standard continuous interpolant [7] as

[TABLE]

so that ${\widehat{X}}_{t}$ is the solution of the SDE

[TABLE]

In the following subsections, we state some preliminary lemmas, the key results on stability and strong convergence, and related results on the number of timesteps, introducing various assumptions as required for each. The main proofs are deferred to Section 5.

2.2 Stability

Assumption 1** (Local Lipschitz and linear growth).**

$f$* and $g$ are both locally Lipschitz, so that for any $R\!>\!0$ there is a constant $C_{R}$ such that*

[TABLE]

for all $x,y\in{\mathbb{R}}^{m}$ with $\|x\|,\|y\|\leq R$. Furthermore, there exist constants $\alpha,\beta>0$ such that for all $x\in{\mathbb{R}}^{m}$, $f$ satisfies the dissipativity condition:

[TABLE]

and $g$ is globally bounded and non-degenerate:

[TABLE]

Lemma 1** (SDE stability).**

If the SDE satisfies Assumption 1 with $X_{0}=x_{0},$ then for all $p\in(0,\infty),$ there is a constant $C_{p}$ which only depends on $x_{0}$ and $p$ such that, $\forall\,t\geq 0,$

[TABLE]

Proof.

The proof is for $p\!\geq\!2$; the result for $0\!<\!p\!<\!2$ follows from Hölder’s inequality.

By Itô’s formula, we have for any $t\geq 0,$

[TABLE]

and then by taking expectations on both sides and using the dissipativity property (7) and boundedness (8), we obtain

[TABLE]

Young inequality (33) implies that

[TABLE]

where $c_{p}=\left(\frac{p-2}{p\alpha}\right)^{p/2-1}(\beta(p+1))^{p/2}.$ Therefore, we obtain

[TABLE]

and we can conclude that

[TABLE]

∎

We now specify the critical assumption about the adaptive timestep.

Assumption 2** (Adaptive timestep).**

The adaptive timestep function $h:{\mathbb{R}}^{m}\rightarrow(0,{h_{\rm max}}]$ is continuous and bounded, with $0<{h_{\rm max}}<\infty$, and there exist constants $\alpha,\beta>0$ such that for all $x\in{\mathbb{R}}^{m}$, $h$ satisfies the inequality

[TABLE]

Note that if another timestep function $h^{\delta}(x)$ is smaller than $h(x)$, then $h^{\delta}(x)$ also satisfies this Assumption. Note also that the form of (9), which is motivated by the requirements of the proof of the next theorem, is very similar to (7). Indeed, if (9) is satisfied then (7) is also true for the same values of $\alpha$ and $\beta$. Compared with the condition in the finite time analysis [1], we need the additional ${h_{\rm max}}$ upper bound to achieve the following result.

Theorem 1** (Infinite time stability).**

If the SDE satisfies Assumption 1, and the timestep function $h$ satisfies Assumption 2, then for all $p\in(0,\infty)$ there exists a constant $C_{p}$ which depends solely on $p,$ $x_{0}$, ${h_{\rm max}}$ and the constants $\alpha,\beta$ in Assumption 2 such that, $\forall t\geq 0,$

[TABLE]

Proof.

The proof is deferred to Section 5. ∎

To bound the expected number of timesteps, we require an assumption on how quickly $h(x)$ can approach zero as $\|x\|\rightarrow\infty$.

Assumption 3** (Timestep lower bound).**

There exist constants $\xi,\zeta,q\!>\!0$, such that the adaptive timestep function satisfies the inequality

[TABLE]

Given this assumption, we obtain the following lemma.

Lemma 2** (Bounded timestep moments).**

If the SDE satisfies Assumption 1, and the timestep function $h$ satisfies Assumptions 2 and 3, then for all $T,p\in(0,\infty)$ there exists a constant $C_{h,p}$ which depends on $p$ and $C_{p}$ in Theorem 1 such that

[TABLE]

where $\displaystyle N_{T}=\min\{n:t_{n}\geq T\}$ is the number of timesteps required by a path approximation.

Proof.

By Assumption 3, we have

[TABLE]

Therefore, by Jensen’s inequality, we obtain

[TABLE]

and the result is then an immediate consequence of Theorem 1. ∎

2.3 Strong convergence

Standard strong convergence analysis for an approximation with a uniform timestep $h$ considers the limit $h\!\rightarrow\!0$. This clearly needs to be modified when using an adaptive timestep, and we will instead consider a timestep function $h^{\delta}$ controlled by a scalar parameter $0\!<\!\delta\!\leq\!1$, and consider the limit $\delta\!\rightarrow\!0$.

In our analysis, we will make the following assumption.

Assumption 4**.**

The timestep function $h^{\delta}$ satisfies the inequalities

[TABLE]

and $h$ satisfies Assumption 2.

To prove an order of strong convergence requires new assumptions on $f$ and $g$:

Assumption 5** (Contractive Lipschitz properties).**

For some fixed $p^{*}\in[2,\infty),$ there exist constants $\lambda,\eta\!>\!0$ such that for all $x,y\in{\mathbb{R}}^{m}$, $f$ and $g$ satisfy the contractive Lipschitz condition:

[TABLE]

and $g$ satisfies the Lipschitz condition:

[TABLE]

In addition, $f$ satisfies the local polynomial growth Lipschitz condition

[TABLE]

for some $\gamma,\mu,q>0$.

Note that we will prove that this Assumption ensures that two solutions to this SDE starting from different places but driven by the same Brownian motion, will come together exponentially. That means the error made on previous time steps will decay exponentially and then we can prove a uniform bound for the strong error. If the drift and volatility are differentiable, the following assumption is equivalent to Assumption 5, and usually easier to check in practice.

Assumption 6** (Contractive Lipschitz properties).**

For some fixed $p^{*}\in[2,\infty),$ there exists a constant $\lambda\!>\!0$ such that for all $x,e\in{\mathbb{R}}^{m}$ with $\|e\|\!=\!1$, $f$ and $g$ are differentiable and satisfy the contractive Lipschitz condition:

[TABLE]

and $g$ satisfies the Lipschitz condition:

[TABLE]

and in addition $f$ satisfies the local polynomial growth Lipschitz condition

[TABLE]

for some $\gamma,\mu,q>0$.

Lemma 3** (SDE contractivity).**

If the SDE satisfies Assumption 5, then for $p\in[2,p^{*}]$ any two solutions to the SDE: $X_{t}$ and $Y_{t},$ driven by the same Brownian motion but starting from $x_{0}$ and $y_{0}$, where $x_{0}\neq y_{0}$, satisfy, $\forall\ t>0,$

[TABLE]

Proof.

First, we can define $e_{t}\triangleq X_{t}-Y_{t},$ and since $X_{t}$ and $Y_{t}$ are driven by the same Brownian motion, we get

[TABLE]

By Itô’s formula, we have for any $0<t\leq T,$

[TABLE]

Therefore, by taking expectations on both sides and using the contractive Lipschitz property (11), we obtain that

[TABLE]

∎

Theorem 2** (Strong convergence order).**

If the SDE satisfies Assumption 5, and the timestep function $h^{\delta}$ satisfies Assumption 4, then for all $p\in(0,p^{*}]$ there exists a constant $C_{p}$ such that, $\forall t\geq 0,$

[TABLE]

Proof.

The proof is deferred to Section 5. ∎

Note that this theorem implies the half-order weak convergence, but the numerical results shows the standard first order weak convergence. We do not provide the proof for the first order weak convergence since the strong convergence result is sufficient to extend this scheme to MLMC.

Lemma 4** (Number of timesteps).**

If the SDE satisfies Assumption 1, and the timestep function $h^{\delta}$ satisfies Assumption 4 with $h(x)$ satisfying Assumption 3, then for all $T,p\in(0,\infty)$ there exists a constant $C_{h,p}$ same as in Lemma 2 such that

[TABLE]

where $N_{T}$ is again the number of timesteps required by a path approximation.

Proof.

The proof is very similar to the proof of Lemma 2, noting that

[TABLE]

due to Assumptions 3 and 4. ∎

First order strong convergence is achievable for Langevin SDEs in which $m\!=\!d$ and $g$ is the identity matrix $I_{m}$, but this requires stronger assumptions on the drift $f$.

Assumption 7** (Enhanced contractive Lipschitz properties).**

There exists a constant $\lambda\!>\!0$ such that for all $x,y\in{\mathbb{R}}^{m}$, $f$ satisfies the contractive one-sided Lipschitz condition:

[TABLE]

In addition, $f$ is differentiable, and $f$ and $\nabla f$ satisfy the local polynomial growth Lipschitz condition

[TABLE]

for some $\gamma,\mu,q>0$.

Lemma 5**.**

*If $f$ satisfies Assumption 7, then for any $x,y,v\in{\mathbb{R}}^{m}$

[TABLE]

where the remainder term has the bound

[TABLE]

Proof.

If we define the scalar function $u(\lambda)$ for $0\!\leq\!\lambda\!\leq\!1$ by

[TABLE]

then $u(\lambda)$ is continuously differentiable, and by the Mean Value Theorem $u(1)\!-\!u(0)=u^{\prime}(\lambda^{*})$ for some $0\!<\!\lambda^{*}\!<\!1$, which implies that

[TABLE]

The final result then follows from the Lipschitz property of $\nabla f$. ∎

We now state the theorem on improved strong convergence.

Theorem 3** (Strong convergence for Langevin SDEs).**

If $m\!=\!d$, $g\equiv I_{m}$, $f$ satisfies Assumption 7, and the timestep function $h^{\delta}$ satisfies Assumption 4, then for all $p\in(0,\infty)$ there exists a constant $C_{p}$ such that, $\forall\ t\geq 0,$

[TABLE]

Proof.

The proof is deferred to Section 5. ∎

3 Adaptive Multilevel Monte Carlo for invariant distributions

We are interested in the problem of approximating:

[TABLE]

where $\pi$ is the invariant measure of the SDE (1). Numerically, we can approximate this quantity by simulating ${\mathbb{E}}\left[\varphi(X_{T})\right]$ for a sufficiently large $T.$ In the following subsections, we will introduce our adaptive multilevel Monte Carlo algorithm and its numerical analysis.

3.1 Algorithm

To estimate ${\mathbb{E}}\left[\varphi(X_{T})\right],$ the simplest Monte Carlo estimator is

[TABLE]

where ${\widehat{X}}_{T}^{(n)}$ is the terminal value of the $n$th numerical path in the time interval $[0,T]$ using a suitable adaptive function $h^{\delta}.$ It can be extended to Multilevel Monte Carlo by using non-nested timesteps as explained in [4]. Consider the identity

[TABLE]

where $\varphi_{\ell}:=\varphi({\widehat{X}}_{T}^{\ell})$ with ${\widehat{X}}_{T}^{\ell}$ being the numerical estimator of $X_{T},$ which uses adaptive function $h^{\delta}$ with $\delta\!=\!M^{-\ell}$ for some fixed $M\!>\!1.$ Then the standard MLMC estimator is the following telescoping sum:

[TABLE]

where ${\widehat{X}}_{T}^{(n,\ell)}$ is the terminal value of the $n$th numerical path in the time interval $[0,T]$ using a suitable adaptive function $h^{\delta}$ with $\delta=M^{-\ell}.$

Different from the standard MLMC with fixed time interval $[0,T],$ we now allow different levels to have a different length of time interval $T_{\ell},$ satisfying $0<T_{0}<T_{1}<...<T_{\ell}<...<T_{L}=T,$ which means that as level $\ell$ increases, we obtain a better approximation not only by using smaller timesteps but also by simulating a longer time interval. However, the difficulty is how to construct a good coupling on each level $\ell$ since the fine path and coarse path have different lengths of time interval $T_{\ell}$ and $T_{\ell-1}.$

Following the idea of Glynn and Rhee [5] to estimate the invariant measure of some Markov chains, we perform the coupling by starting a level $\ell$ fine path simulation at time $t^{f}_{0}=-T_{\ell}$ and a coarse path simulation at time $t^{c}_{0}=-T_{\ell-1}$ and terminate both paths at $t=0.$ Since the drift $f$ and volatility $g$ do not depend explicitly on time $t,$ the distribution of the numerical solution simulated on the time interval $[-T_{\ell},0]$ is same as one simulated on $[0,T_{\ell}].$ The key point here is that the fine path and coarse path share the same driving Brownian motion during the overlap time interval $[-T_{\ell-1},0].$ Owing to the result of Lemma 3, two solutions to the SDE satisfying Assumption 5, starting from different initial points and driven by the same Brownian motion will converge exponentially. Therefore, the fact that different levels terminate at the same time is crucial to the variance reduction of the multilevel scheme.

Our new multilevel scheme still has the identity (19) but with $\varphi_{\ell}=\varphi({\widehat{X}}_{0}^{\ell})$ with ${\widehat{X}}^{\ell}_{0}$ being the terminal value of the numerical path approximation on the time interval $[-T_{\ell},0]$ using adaptive function $h^{\delta}$ with $\delta\!=\!M^{-\ell}.$ The corresponding new MLMC estimator is

[TABLE]

where ${\widehat{X}}_{0}^{(n,\ell)}$ is the terminal value of the $n$th numerical path through time interval $[-T_{\ell},0]$ using adaptive function $h^{\delta}$ with $\delta=M^{-\ell}.$ Figure 1 and Algorithm 1 illustrate the detailed implementation of a single adaptive MLMC sample using a non-nested adaptive timestep on level $\ell$ with $M=2.$

3.2 Numerical analysis

First, we state the exponential convergence to the invariant measure of the original SDEs, which can help us to measure the approximation error caused by truncating the infinite time interval.

Lemma 6** (Exponential convergence).**

If the SDE satisfies Assumption 1 and $\varphi$ satisfies the Lipschitz condition: there exists a constant $\kappa>0$ such that

[TABLE]

then this SDE is ergodic and has a unique invariant measure $\pi$ and there exist constants $\mu^{*},\,\lambda^{*}>0$ such that

[TABLE]

If the SDE additionally satisfies Assumption 5, then there exists a constant $\mu>0$ depending on $x_{0},$ $\kappa$ and $C_{1}$ in Lemma 1 such that

[TABLE]

Proof.

The proof of existence and uniqueness of the invariant measure and convergence (22) is given in Theorem 4.4 in [10] and Theorem 6.1 in [11].

Next, we can define a new random variable $Y_{0}$ which follows the invariant measure $\pi,$ then the solution $Y_{t}$ to the SDE with the initial value $Y_{0}$ will also follows the invariant measure for any $t\!>\!0.$ Therefore, by Lipschitz property of $\varphi$ and Lemma 1 and 3, there exists a constant $\mu\!>\!0$ such that

[TABLE]

∎

Note that $\lambda$ is easier to estimate than $\lambda^{*}$ through Assumption 6.

Lemma 7** (Variance of MLMC corrections for bounded volatility).**

If $\varphi$ satisfies the Lipschitz condition (21), the SDE satisfies Assumption 5 and the timestep function $h^{\delta}$ satisfies Assumption 4 with $\delta=M^{-\ell}$ for each level, then for each level $\ell,$ there exist constants $c_{1}$ and $c_{2}$ such that the variance of correction $V_{\ell}:=\mathbb{V}\left[\varphi({\widehat{X}}_{0}^{\ell})-\varphi({\widehat{X}}_{0}^{\ell-1})\right]$ satisfies

[TABLE]

Proof.

Lipschitz condition (21) implies

[TABLE]

${\widehat{X}}^{\ell}$ and ${\widehat{X}}^{\ell-1}$ share the same driving Brownian motion from $-T_{\ell-1}$ to $0.$ We can define the corresponding solution to the SDE (1) starting from $x_{0}$ and driven by the same Brownian motion as ${\widehat{X}}^{\ell-1}$ through time interval $[-T_{\ell-1},0]$ by $X^{c},$ and the solution starting from $x_{0}$ driven by the same Brownian motion as ${\widehat{X}}^{\ell}$ through time interval $[-T_{\ell},0]$ by $X^{f}.$

Then, by Jensen’s inequality, we obtain that

[TABLE]

where

[TABLE]

Theorem 2 implies that there exist a constant $C_{2}$ which does not depend on $T_{\ell}$ such that

[TABLE]

and Lemma 1 and Lemma 3 imply that there exists a constant $C$ depending on $x_{0}$ and $C_{2}$ in Lemma 1 such that

[TABLE]

Finally, we obtain the desired result:

[TABLE]

∎

Note that if we set

[TABLE]

then

[TABLE]

which has the same order of magnitude as the variance bound for the standard finite time interval MLMC considered in Part I, [1].

We define the computational cost of a path simulation to be equal to the number of timesteps. Hence, due to Lemma 2 and (25), there exists a constant $C_{0}$ such that the expected cost of a single MLMC sample on level $\ell$ is bounded by $C_{0}(\ell\!+\!1)M^{\ell}$. Given this, we obtain the following theorem for the complexity of the MLMC algorithm to achieve a specified Mean Square Error accuracy.

Theorem 4** (MLMC for invariant measure).**

If $\varphi$ satisfies the Lipschitz condition (21), the SDE satisfies Assumption 5 and the timestep function $h^{\delta}$ satisfies Assumption 4 with $\delta\!=\!M^{-\ell}$ for each level, then by choosing suitable values for $L$ and $T_{\ell},\,N_{\ell}$ for each level $\ell,$ there exists a constant $c_{3}$ such that the MLMC estimator (20) has a mean square error (MSE) with bound

[TABLE]

and an expected computational cost $C$ with bound

[TABLE]

Proof.

By Jensen’s inequality, the mean square error can be decomposed into three parts:

[TABLE]

which enables us to achieve the MSE bound by bounding each part by $\varepsilon^{2}/3$.

Lemma 23 implies that

[TABLE]

provided we ensure that

[TABLE]

If we set $T_{\ell}$ according to (25), this is achieved by requiring

[TABLE]

By Theorem 2 and the Lipschitz property (21) of $\varphi$, there exists a constant $C_{2}$ such that

[TABLE]

provided

[TABLE]

Therefore, combining the requirements (27) and (28), we choose to define

[TABLE]

giving $L=O(|\log\varepsilon|)$ as $\varepsilon\rightarrow 0$.

Next, we need to choose the number of samples $N_{\ell}$ for each level. We aim to minimize the total expected computational cost, which is bounded by

[TABLE]

while at the same time ensuring that the total variance satisfies the bound

[TABLE]

Using a Lagrange multiplier, it is found that the optimal solution to the constrained optimization problem

[TABLE]

when the $N_{\ell}$ are treated as real variables is

[TABLE]

Rounding this up to an integer by defining

[TABLE]

we ensure that the required variance bound is satisfied. The resulting cost is then bounded by

[TABLE]

Since

[TABLE]

and

[TABLE]

we obtain the desired final result that there exists a constant $c_{3}$ such that

[TABLE]

∎

For Langevin SDEs, the computational cost can be reduced to $O(\varepsilon^{-2}).$

Theorem 5** (Langevin SDEs).**

If $\varphi$ satisfies the Lipschitz condition (21), and for the SDE, $m\!=\!d$, $g\equiv I_{m}$, $f$ satisfies Assumption 7, and the timestep function $h^{\delta}$ satisfies Assumption 4 with $\delta=M^{-\ell}$ for each level, then for each level $\ell,$ there exist constants $c_{1}$ and $c_{2}$ such that

[TABLE]

Furthermore, by choosing suitable $T_{\ell}$ and $N_{\ell}$ for each level $\ell$ in the MLMC estimator (20), one can achieve the MSE bound $\varepsilon^{2}$ at an expected computational cost bounded by

[TABLE]

for some constant $c_{3}\!>\!0$.

Proof.

Following a similar argument to the proof of Lemma 24, Theorem 3 implies $V_{\ell}\leq c_{1}\,M^{-2\ell}+c_{2}\,{\rm e}^{-2\lambda T_{\ell-1}},$ and by choosing $T_{\ell}$ to be

[TABLE]

we obtain $V_{\ell}\leq(c_{1}\!+\!c_{2})M^{-2\ell}.$ The computational cost of a single MLMC sample on level $\ell$ satisfies

[TABLE]

for any $0\!<\!\epsilon\!\ll\!1$ and some $C\!>\!0$. Therefore, the standard MLMC Theorem 1 in [3] is applicable with $\gamma\!<\!\beta$, giving an $O(\varepsilon^{-2})$ complexity. ∎

Note that the choice of $T_{\ell}$ (31) for the Langevin equation is different from (25) for SDEs with bounded volatility. In other words, the strong convergence result and the contractive convergence rate $\lambda$ together determine $T_{\ell}$. The difference in the variance convergence rate also affects the choice of $M$. Based on the analysis in [2], the optimal $M$ for SDEs with general $g$ is in the range $4-8$, while in the Langevin case the optimal $M$ is around $2$.

4 Numerical Experiment

In this section we present numerical results for the following scalar SDE:

[TABLE]

which satisfies the dissipativity condition (7) and the contractive condition (17). Our interest is to compute $\pi(\varphi)$ where $\varphi(x)=\|x\|$ satisfying Lipschitz condition.

Since the probability density function $\pi$ is

[TABLE]

we can use numerical integration to calculate an approximate value: $\varphi(\pi)\approx 0.44115$ with accuracy $10^{-5},$ and use this value as a benchmark for our numerical tests.

Following condition (9) we can set $x_{0}\!=\!0$, ${h_{\rm max}}\!=\!1$, $M\!=\!2$ and choose the adaptive function $h,\ h^{\delta}$ to be

[TABLE]

Next we need to determine $T_{\ell}$ for each level. Linear perturbations to the SDE satisfy the ODE:

[TABLE]

and therefore $\lambda\!\geq\!1$. Hence we choose to use

[TABLE]

to ensure that the truncation error is acceptably small.

The variance result (30) for the Langevin equation is illustrated in Figure 2. The exponential term dominates the variance initially, but as $T_{\ell}$ increase, the $M^{-2\ell}$ term eventually becomes the major part of the variance.

Figure 3 presents the MLMC results. The top right plot shows first order convergence for the weak error and the top left plot shows second order convergence for the multilevel correction variance. Hence the computational cost for RMS accuracy $\varepsilon$ is $O(\varepsilon^{-2})$ which is verified in the bottom right plot, while the bottom left plot shows the number of MLMC samples on each level as a function of the target accuracy.

5 Proofs

This section has the proofs of the three main theorems in this paper, one on stability, and two on the order of strong convergence.

5.1 Preliminaries

In this subsection, we introduce some inequalities and results we use frequently in the following sections.

5.1.1 Young inequality

For any $\zeta\!>\!0$ and $P,Q\!>\!1$ satisfying $\frac{1}{P}\!+\!\frac{1}{Q}=1,$ the following inequality holds for any $A,B\!>\!0,$

[TABLE]

In this paper, we use two particular cases. First, we take $P\!=\!Q\!=\!2$, and $\zeta^{2}\!=\!2\xi\!>\!0$ and get

[TABLE]

Second, for any $p\geq 2$ and $\xi>0,$ we take $P\!=\!\frac{p}{p-2},$ $Q\!=\!\frac{p}{2},$ $A\!=\!a^{p-2},$ $B\!=\!b^{2}$ and $\zeta\!=\!\xi^{(p-2)/p},$ and get

[TABLE]

When using these in proofs, we often keep $\xi$ arbitrary initially and choose it later to make one term sufficiently small, as needed.

5.1.2 Jensen inequality

One variant of the Jensen inequality is

[TABLE]

where $p\geq 1$ and $h_{k},\ u_{k}\geq 0.$ Its continuous version is

[TABLE]

5.1.3 Exponentially weighted supremum

For simplicity, for $\alpha\!>\!0$, we can define ${\widehat{M}}_{t}^{\alpha,p}\triangleq\sup_{0\leq s\leq t}{\rm e}^{\alpha ps}\|{\widehat{X}}_{s}\|^{p},$ and then by Young’s inequality (32), for any $\xi\!>\!0,$

[TABLE]

We can also define ${\overline{M}}_{t}^{\alpha,p}\triangleq\sup_{0\leq s\leq t}{\rm e}^{\alpha ps}\|{\overline{X}}_{s}\|^{p},$ which implies

[TABLE]

since ${\overline{X}}_{s}\!=\!{\widehat{X}}_{\underline{s}}$ and $|s\!-\!{\underline{s}}|\leq{h_{\rm max}},$ and

[TABLE]

provided $\gamma\!>\!\alpha\!>\!0$.

5.2 Theorem 1

Proof.

By theorem 1 in Part I [1], we know $T$ is almost surely attainable. Therefore we can directly analyse our discretization scheme without the $K$ truncation which was used in that paper. The proof proceeds in three steps. First, we derive an upper bound for ${\rm e}^{\alpha pt}\|{\widehat{X}}_{t}\|^{p}$. Second, we show that the moments ${\mathbb{E}}[{\widehat{M}}^{\alpha,p}_{t}]$ and ${\mathbb{E}}[{\overline{M}}^{\alpha,p}_{t}]$ are each bounded by $C_{p}{\rm e}^{\alpha pt}$ where $C_{p}$ is a constant which only depends on $p$, $x_{0}$, ${h_{\rm max}}$ and the constants $\alpha,\beta$ in Assumption 2. Finally, we get the uniform bound for ${\mathbb{E}}[\|{\widehat{X}}_{t}\|^{p}]$ and ${\mathbb{E}}[\|{\overline{X}}_{t}\|^{p}].$

The proof is given for $p\!\geq\!4$; the result for $0\!<\!p\!<\!4$ follows from Hölder’s inequality.

**Step 1: ** If we define $\phi(x)\triangleq x\!+\!h(x)f(x)$, then (4) gives

[TABLE]

Using condition (9) for $h$ then gives

[TABLE]

Since $1\!-\!2\alpha h_{n}\leq{\rm e}^{-2\alpha h_{n}}$ and $g$ and $h$ are both bounded, we multiply by ${\rm e}^{2\alpha{t_{n+1}}}$ on both sides to obtain

[TABLE]

Similarly, for the partial timestep from ${\underline{t}}$ to $t$, since $(t\!-\!{\underline{t}})\leq h_{n_{t}},$

[TABLE]

and therefore we obtain

[TABLE]

Summing (39) over multiple timesteps and then adding (41) gives

[TABLE]

Bounding the first summation using a Riemann integral, and re-writing the second as an Itô integral, raising both sides to the power $p/2$ and using Jensen’s inequality, we obtain

[TABLE]

**Step 2: ** For any $0\!\leq\!t\!\leq\!T,$ we take the supremum on both sides of inequality (5.2) and then take the expectation to obtain

[TABLE]

where

[TABLE]

[TABLE]

We now consider $I_{1},I_{2},I_{3},I_{4},I_{5}$ in turn.

[TABLE]

By the Burkholder-Davis-Gundy inequality, there exist constants $C_{p}^{1}$ such that

[TABLE]

Due to condition (9), for $u\!<\!{\underline{t}}$ we have

[TABLE]

and hence by Jensen’s inequality and the boundedness condition (8) of $g$ , we obtain

[TABLE]

Therefore, using Jensen’s inequality (35) with $\gamma\!=\!2\alpha$, followed by (38) with $\gamma\!=\!(1+4/p)\alpha$ and then (36) with $\xi\!=\!{\rm e}^{-\alpha pt/2}\zeta$, there exists a constant $C_{p}^{2}$ which is linearly dependent on $\zeta^{-1}$ such that

[TABLE]

Using Jensen inequality (34), we obtain

[TABLE]

where we define

[TABLE]

so that ${\mathbb{E}}[\|{\Delta W}_{k}\|^{p}]\leq c_{p}h_{k}^{p/2}$.

In considering $I_{4}$, we start by observing that for $t_{k}\!\leq\!s\!<\!t_{k+1}$

[TABLE]

In addition, using (40) and following the same argument as for $I_{2}$, we have

[TABLE]

Therefore, combining the estimation (43), (38) with $\gamma\!=\!2\alpha$ and (36) with $\xi={\rm e}^{-\alpha pt/2}\zeta$, there exists $C_{p}^{3}$ which is linearly dependent on $\zeta^{-1}$ such that

[TABLE]

Similarly, again using the same definition for $c_{p}$, we have

[TABLE]

Collecting together the bounds for $I_{1},I_{2},I_{3},I_{4},I_{5}$, we conclude that we can choose $\zeta\!>\!0$ sufficiently small so that there exist constants $C^{4}_{p}$ and $C_{p}^{5}$ such that

[TABLE]

and hence

[TABLE]

**Step 3: ** Due to the definition of ${\overline{M}}_{t}^{\alpha,p}$ and inequality (37), for any $t\geq 0,$

[TABLE]

and similarly

[TABLE]

∎

5.3 Theorem 2

Proof.

The approach which is followed is to bound the approximation error $e_{t}\triangleq{\widehat{X}}_{t}-X_{t}$ by terms which depend on ${\widehat{X}}_{s}\!-\!{\overline{X}}_{s}$, and then use local analysis within each timestep to bound these. The proof is for $2\leq\!p\!\leq\!p^{*}$; the result for $0\!<\!p\!<\!2$ follows from Hölder’s inequality.

We start by combining the original SDE (1) with (6) to obtain

[TABLE]

and then by Itô’s formula and Young’s inequality (32), together with $e_{0}=0,$ and $\lambda,\ \eta$ as defined in Assumption 5, we get

[TABLE]

Using the conditions in Assumption 5, (12) implies that

[TABLE]

(11) and (12) imply that

[TABLE]

(13) and Young inequality (32) implies that

[TABLE]

where $L(x,y)\triangleq\gamma(\|x\|^{q}+\|y\|^{q})+\mu$. Hence,

[TABLE]

where $\hat{L}(x,y)=\frac{2}{\lambda}L(x,y)^{2}\!+\!\frac{(p-1)\eta}{2}+\frac{\eta^{2}(p-1)^{2}}{\lambda}.$ Young inequality (33) implies

[TABLE]

Taking the expectation of each side yields

[TABLE]

By the Hölder inequality,

[TABLE]

and ${\mathbb{E}}\left[\hat{L}({\widehat{X}}_{s},{\overline{X}}_{s})^{p}\right]$ can be bounded by a constant $C_{p}^{1}$ due to the stability property in Theorem 1.

For any $s\!\in\![0,T]$, ${\widehat{X}}_{s}\!-\!{\overline{X}}_{s}=f({\widehat{X}}_{{\underline{s}}})(s\!-\!{\underline{s}})+g({\widehat{X}}_{{\underline{s}}})(W_{s}\!-\!W_{\underline{s}})$, and hence, by a combination of Jensen and Hölder inequalities, we get

[TABLE]

${\mathbb{E}}[\|f({\widehat{X}}_{{\underline{s}}})\|^{4p}]$ and ${\mathbb{E}}[\|g({\widehat{X}}_{{\underline{s}}})\|^{4p}]$ are both finite, due to stability and the polynomial bounds on the growth of $f(x)$ and $g(x)$. Furthermore, we have ${\mathbb{E}}[(s\!-\!{\underline{s}})^{4p}]\leq(\delta{h_{\rm max}})^{4p}\leq{h_{\rm max}}^{4p}\delta^{2p}$, and by standard results there is a constant $c_{p}$ such that ${\mathbb{E}}[\|W_{s}\!-\!W_{\underline{s}}\|^{4p}]={\mathbb{E}}[\ {\mathbb{E}}[\|W_{s}\!-\!W_{\underline{s}}\|^{4p}\,|\,{\cal F}_{\underline{s}}]\ ]\leq c_{p}(\delta{h_{\rm max}})^{2p}$. Hence, there exists a constant $C_{p}^{2}\!>\!0$ such that ${\mathbb{E}}[\,\|{\widehat{X}}_{s}\!-\!{\overline{X}}_{s}\|^{2p}]\leq C^{2}_{p}\,\delta^{p}$, and therefore equation (5.3) gives us

[TABLE]

which provides the final result:

[TABLE]

∎

5.4 Theorem 3

Proof.

The proof is given for $p\!\geq\!4$; the result for $0\!<\!p\!<\!4$ follows from Hölder’s inequality.

The error $e_{t}\triangleq{\widehat{X}}_{t}-X_{t}$ satisfies the SDE ${\rm d}e_{t}=\left(f({\overline{X}}_{t})\!-\!f(X_{t})\right){\rm d}t$ and hence by Itô’s formula,

[TABLE]

due to the one-sided Lipschitz condition (17), so therefore

[TABLE]

Within a single timestep, ${\widehat{X}}_{u}\!-\!{\overline{X}}_{u}=f({\overline{X}}_{u})(u\!-\!{\underline{u}})+(W_{u}\!-\!W_{\underline{u}})$, and therefore, following a similar approach to the proof of Theorem 4 in [1], Lemma 5 gives

[TABLE]

where $|R_{u}|\leq\left(\gamma\,(\|{\widehat{X}}_{u}\|^{q}\!+\!\|{\overline{X}}_{u}\|^{q})+\mu\right)\,\|e_{u}\|\,\|{\widehat{X}}_{u}\!-\!{\overline{X}}_{u}\|^{2}$, and hence

[TABLE]

where

[TABLE]

We now bound $I_{1},I_{2},I_{3},I_{4},I_{5}$ in turn. Noting that $u\!-\!{\underline{u}}\leq\delta{h_{\rm max}}$, by Young inequality (32) and Jensen inequality (35), we obtain

[TABLE]

The last integral is finite because of stability and the polynomial bounds on the growth of both $f$ and $\nabla f$, and hence there is a constant $C^{1}_{p}$ such that

[TABLE]

Similarly, using the Young inequality (32), Jensen inequality (35) and the Hölder inequality, we obtain

[TABLE]

where $L({\widehat{X}}_{u},{\overline{X}}_{u})\equiv\gamma\,(\|{\widehat{X}}_{u}\|^{q}\!+\!\|{\overline{X}}_{u}\|^{q})\!+\!\mu$. Hence, using stability and bounds on ${\mathbb{E}}\left[\|{\widehat{X}}_{u}\!-\!{\overline{X}}_{u}\|^{4p}\right]$ from the proof of Theorem 2, there is a constant $C^{2}_{p}$ such that

[TABLE]

The fact that ${\rm e}^{2\lambda u}\!-\!{\rm e}^{2\lambda{\underline{u}}}\leq 2\lambda\,{\rm e}^{2\lambda u}(u-{\underline{u}}),$ and Jensen inequality (35) gives

[TABLE]

and using the approach to estimating $I_{1},$ there is similarly a constant $C_{p}^{3}$ such that

[TABLE]

For the next term, $I_{4}$, we start by bounding $\|e_{s}\!-\!e_{\underline{s}}\|$. Since

[TABLE]

by Jensen’s inequality and Assumption 7 it follows that

[TABLE]

where $L({\overline{X}}_{u},X_{u})\equiv\gamma(\|{\overline{X}}_{u}\|^{q}\!+\!\|X_{u}\|^{q})+\mu$. We again have an $O(\delta^{p/2})$ bound for ${\mathbb{E}}[\|{\widehat{X}}_{s}\!-\!{\overline{X}}_{s}\|^{p}]$, while Theorem 2 proves that there is a constant $c_{p}$ such that

[TABLE]

Combining these, and using the Hölder inequality and the finite bound for ${\mathbb{E}}[L^{p}({\overline{X}}_{u},X_{u})]$ for all $p\!\geq\!2$, due to the usual stability results, we find that there is a different constant $c_{p}$ such that

[TABLE]

Now, by Jensen inequality (35), we obtain

[TABLE]

so using the Hölder inequality and the usual stability bounds, we conclude that there is a constant $C^{4}_{p}$ such that

[TABLE]

Lastly, considering that

[TABLE]

in the time interval $[t_{n},t_{n+1}]$ conditioned on $\mathcal{F}_{t_{n}},$ we have

[TABLE]

and therefore we can split $I_{5}$ into two parts:

[TABLE]

By the Burkholder-Davis-Gundy inequality, Young inequality (32) and Jensen inequality (35), there exists a constant $c_{p}^{1}$ such that

[TABLE]

Hence, there exists a constant $C^{51}_{p}$ such that

[TABLE]

Finally, for $I_{52}$, the Young inequality (32) and Jensen inequality (35) implies

[TABLE]

Thus, we find that there exists a constant $C^{52}_{p}>0$ such that

[TABLE]

Combining the five bounds, and choosing $\xi$ to be sufficiently small, we conclude that there is a constant $C_{p}$ such that

[TABLE]

and therefore

[TABLE]

∎

6 Conclusions and future work

In this paper, we have developed the adaptive Euler-Maruyama method for a class of ergodic SDEs with non-globally Lipschitz drifts and shown the moments and strong error of the numerical solutions are uniformly bounded in time $T.$ Moreover, we extend this adaptive scheme to Multi-level Monte Carlo for expectations with respect to the invariant measure by using different time intervals on different levels, and constructing an efficient coupling of the fine path and coarse paths with different simulation times. Numerical experiments support the theoretical analysis.

One direction for extension of the theory is to address SDEs which don’t satisfy the contractive property. For example, the SDEs arising from the double-well potential energy, in which case, the fine path and the coarse path may converge to different wells resulting in a poor coupling. Another possibility is to apply adaptive methods to SDEs with a discontinuous drift.