Strong convergence of multiscale truncated Euler-Maruyama method for super-linear slow-fast stochastic differential equations

Stochastic modelling plays an essential role in many branches of science and industry. Especially, super-linear stochastic differential equations (SDEs) are usually used to describe real-world systems in various applications, for examples, the stochastic Lotka-Volterra model in biology for the population growth [33], the elasticity of volatility model arising in finance for the asset price [26] and the stochastic Ginzburg-Landau equation stemming from statistical physics in the study of phase transitions [24]. In many fields, various factors change at different rates: some vary rapidly whereas others evolve slowly. As a result the separation of fast and slow time scales arises in chemistry, fluid dynamics, biology, physics, finance and other fields [5, 14, 17, 25]. Stochastic systems with this characteristic are studied extensively [10, 35, 36, 46] and are often modeled by the following slow-fast SDEs (SFSDEs),

with initial value $(x^{\varepsilon}(0),y^{\varepsilon}(0))=(x_{0},y_{0})\in\mathbb{R}^{n_{1}}\times\mathbb{R}^{n_{2}}$. Here, coefficients

are continuous, and $\{W^{1}(t)\}_{t\geq 0}$ and $\{W^{2}(t)\}_{t\geq 0}$ represent mutually independent $d_{1}$-dimensional and $d_{2}$-dimensional Brownian motions, respectively. The parameter $\varepsilon>0$ represents the ratio of nature time scales between $x^{\varepsilon}(t)$ and $y^{\varepsilon}(t)$. Especially, as $\varepsilon\ll 1$, $x^{\varepsilon}(t)$ and $y^{\varepsilon}(t)$ are called the slow component and fast component, respectively.

In various applications the time evolution of the slow component $x^{\varepsilon}(t)$ is under the spotlight. Due to the existence of super-linear coefficients and multiple time scales as well as the coupling of fast and slow components, it is almost impossible to anticipate the dynamics of slow component directly by solving the full system. Therefore, numerical methods or approximation techniques become the efficient tools. However, on account of the wide separation in time scale, standard computational schemes may be no longer applicable. Hence, our main aim is to construct the appropriate multiscale numerical scheme to approximate the slow component.

On analytical grounds, the averaging principle is essential to describe the asymptotic behavior of the slow component as $\varepsilon\rightarrow 0$. Precisely, assume that the frozen equation is described by

with initial value $y^{x}(0)=y_{0}$, where $x$ is regarded as a parameter. Let $y^{x,y_{0}}(s)$ denote the solution of the frozen equation (1.2) with initial value $y_{0}$. If the transition semigroup of $y^{x,y_{0}}(t)$ has a unique invariant probability measure and the following integral

exists. The averaging principle reveals that under suitable assumptions on the coefficients, the slow component $x^{\varepsilon}(t)$ converges to $\bar{x}(t)$, which is the solution of

Much of the study of the averaging principle can be traced back to the original work of Khasminskii [23] for a class of diffusion processes, in which an averaging principle was showed in weak sense. Subsequently, fruitful results on averaging principle are developed. The convergence between the solution of the averaged system and the exact slow component in probability was yielded in [43, 44]. The strong convergence in mean square was obtained in [16], and was further developed under an ergodicity-type assumption in [15]. Furthermore, order $1/4$ convergence rate was obtained in [13], order $1/2$ strong convergence rate and order $1$ in weak sense were yielded for the degenerate SFSDEs with the deterministic slow component in [8]. The order $1/2$ strong convergence rate was also yielded for SFSDEs whose coefficients depend on small parameter $\varepsilon$ in [30]. Recently, the averaging principle for SFSDEs with the super-linear growth coefficients was obtained in [32]. Furthermore, the strong averaging principles have been developed for various kinds of stochastic slow-fast systems, such as jump-diffusion processes [11, 47], stochastic partial differential equations [2, 9, 4], McKean-Vlasov SDEs [37], and so on.

The averaged equation derived from the averaging principle provides a substantial simplification of the original system. If the coefficient $\bar{b}(\cdot)$ is explicitly available, the averaged equation (1.4) can be immediately used to simulate the evolution of $x^{\varepsilon}(t)$. However, it is often impossible or impractical to do this because the dynamics of the frozen equation (1.2) are too complicated to acquire the invariant measure $\mu^{x}$ analytically. The heterogeneous multiscale method (HMM) [6, 7] was proposed to solve the averaged equation. Within the framework of the HMM, the basic idea is to solve the averaged equation (1.4) by a suitable macro solver, in which the averaged coefficient $\bar{b}$ is estimated on the fly by performing a series of constrained microscopic simulations.

Precisely, in 2003 Vanden-Eijnden [41] proposed numerical schemes for Multiscale dynamical systems with stochastic effects. To fully justify the scheme, in 2005, E et al. [8] provided a thorough analysis of convergence and efficiency of the scheme for degenerate SFSDEs of which the slow dynamics is deterministic. In 2006, Givon et al. [13] developed the Projective integration schemes for SFSDEs with the slow diffusion coefficient only depending on the slow component. And in 2008, Givon et al. [12] extended the projective integration schemes to jump-diffusion systems. Furthermore, in 2010, Liu [29] generalized the analysis of E et al. for the fully coupled SFSDEs, whose slow diffusion coefficient depends on both slow and fast components. Bréhier [2, 3] developed the HMM scheme for the slow-fast parabolic stochastic partial differential equations, which promotes the theory in [8] to the case of the infinite dimension.

Even though the averaged equation is not available explicitly, the selection of macro solver relies on its structure. Provided that the coefficients of averaged equation are global Lipschitz continuous, the Euler-Maruyama (EM) scheme is very popular used as the macro solver owing to the simple algebraic structure and the cheap computational cost [41, 13, 8, 29, 12]. However, super-linear growth coefficients are common phenomena in multiscale systems, which may derive a super-linear averaged equation by averaging principle. Hutzenthaler et al. [21] pointed out that the EM approximation errors for a large class of super-linear SDEs diverge to infinity in $p$th moment for any $p\geq 1$. Hence, the classical EM scheme is no longer suitable as macro solver for SFSDEs with super-linear coefficients. Besides, chosen as the macro solver for the sup-linear averaged equation, implicit scheme tends to make the algorithm and implementation more involved [8]. As a consequence, there is an urgent need to construct an explicit multiscale numerical scheme for super-linear SFSDEs.

Fortunately, great achievements have been made in the research of explicit numerical methods for super-linear SDEs, for examples, the tamed EM scheme [20, 22, 38, 39], the tamed Milstein scheme [45], the stopped EM scheme [31], the truncated EM scheme [27, 28, 34] and therein. So far the ability of these modified EM methods to approximate the solutions of super-linear diffusion systems have been shown extensively. Inspired by the above works, using modified EM method we are devoted to constructing an explicit multiscale numerical method suitable for super-linear SFSDEs.

In fact, HMM relies heavily on the structure of the averaged equation for the slow variable. Thus we have to overcome two major obstacles: the unexplicit form and super-linear structure of $\bar{b}$ . Inspired by the idea from [34], we design a truncation device to modify the super-linear coefficient $b$ of original system in advance, so as to achieve the modification of $\bar{b}$. This modification can avoid possible large excursion due to the super-linearity of $\bar{b}$. Then fitting into the framework of HMM, we construct an explicit multiscale numerical scheme involving three subroutines as follows.

• 1.

  The truncated EM (TEM) scheme is selected as the macro solver to evolute the macro dynamics $\bar{x}(t)\approx x^{\varepsilon}(t)$ in which the modified averaged coefficient is required to be estimated at each macro time step.

• 2.

  An appropriate numerical scheme is chosen as the micro solver to solve the frozen equation to produce the data used for approximating the modified coefficient.

• 3.

  An estimator is established to obtain the desired approximation of the modified averaged coefficient.

Following this line, we construct an easily implementable explicit multiscale numerical scheme for a class of super-linear SFSDEs and obtain its strong convergence.

The rest of this paper is organized as follows. Section 2 gives some notations, hypotheses and preliminaries. Section 3 proposes an explicit multiscale numerical method. Section 4 provides some important pre-estimates. Section 5 yields the strong convergence of MTEM scheme. Section 6 focuses on the error analysis of the explicit MTEM scheme and presents an important example. Section 7 shows two numerical examples and carries out some numerical experiments to verify our theoretical results. Section 8 concludes this paper.

Throughout this paper, we use the following notations. Let $(\Omega,\cal{F},\mathbb{P})$ be a complete probability space with a natural filtration $\{\mathcal{F}_{t}\}_{t\geq 0}$ satisfying the usual conditions (i.e. it is right continuous and increasing while $\mathcal{F}_{0}$ contains all $\mathbb{P}$-null sets), and $\mathbb{E}$ be the expectation corresponding to $\mathbb{P}$. Let $|\cdot|$ denote the Euclidean norm in $\mathbb{R}^{n}$ and the trace norm in $\mathbb{R}^{n\times d}.$ If $A$ is a vector or matrix, we denote its transpose by $A^{T}$. For a set $\mathbb{D}$, let $I_{\mathbb{D}}(x)=1$ if $x\in\mathbb{D}$ and $0$ otherwise. We set $\inf\emptyset=\infty$, where $\emptyset$ is empty set. Moreover, for any $a,b\in\mathbb{R}$, we define $a\vee b=\max\{a,b\}$ and $a\wedge b=\min\{a,b\}$. We use $C$ and $C_{p}$ to denote the generic positive constants, which may take different values at different appearances, where $C_{p}$ is used to emphasize that the constant depends on the parameter $p$. In addition, $C,C_{p}$ are independent of constants $\Delta_{1}$, $\Delta_{2}$, $n$ and $M$ that occur in the next section. $C_{R}$ usually denotes some positive function increasing with respect to $R$.

Let $\mathcal{P}_{p}(\mathbb{R}^{n_{2}})$ be the set of all probability measures on $\mathbb{R}^{n_{2}}$ denoted by $\mathcal{P}(\mathbb{R}^{n_{2}})$, with finite $p$-th moment, i.e.,

which is a Polish space under the Wasserstein distance

where $\mathcal{C}(\mu_{1},\mu_{2})$ stands for the set of all probability measures on $\mathbb{R}^{n_{2}}\times\mathbb{R}^{n_{2}}$ with marginals $\mu_{1}$ and $\mu_{2}$, respectively.

To state the main results, we impose some hypotheses on the coefficients $b$, $\sigma$ of slow system and $f$ and $g$ of fast system.

• (S1)

  There exists a constant $\theta_{1}\geq 1$ such that for any $R>0$, $x_{1},x_{2}\in\mathbb{R}^{n_{1}}$ with $|x_{1}|\vee|x_{2}|\leq R$,

  $\displaystyle|b(x_{1},y)-b(x_{2},y)|+|\sigma(x_{1})-\sigma(x_{2})|\leq L_{R}|x_{1}-x_{2}|(1+|y|^{\theta_{1}}),\$

  here $L_{R}$ is a positive constant dependent on $R$.

• (S2 )

  There exist constants $\theta_{2}\geq 1$ and $K_{1}>0$ such that for any $x\in\mathbb{R}^{n_{1}}$ and $y_{1},y_{2}\in\mathbb{R}^{n_{2}}$.

  $\displaystyle|b(x,y_{1})-b(x,y_{2})|\leq K_{1}|y_{1}-y_{2}|\big{(}1+|x|^{\theta_{2}}+|y_{1}|^{\theta_{2}}+|y_{2}|^{\theta_{2}}\big{)}.$

• (S3)

  There exists a constant $K_{2}>0$ such that for any $x\in\mathbb{R}^{n_{1}}$,

  $\displaystyle|\sigma(x)|\leq K_{2}(1+|x|).$

• (S4)

  There exist constants $\theta_{3},\theta_{4}\geq 1$ and $K_{3}>0$ such that for any $x\in\mathbb{R}^{n_{1}},y\in\mathbb{R}^{n_{2}}$,

  $\displaystyle|b(x,y)|\leq K_{3}(1+|x|^{\theta_{3}}+|y|^{\theta_{4}}).$

• (S5)

  There exist constants $K_{4}>0$ and $\lambda>0$ such that for any $x\in\mathbb{R}^{n_{1}},y\in\mathbb{R}^{n_{2}}$,

  $\displaystyle x^{T}b(x,y)\leq K_{4}(1+|x|^{2})+\lambda|y|^{2}.\$

• (F1)

  The functions $f$ and $g$ are globally Lipschitz continuous, namely, for any $x_{1},x_{2}\in\mathbb{R}^{n_{1}}$ and $y_{1},y_{2}\in\mathbb{R}^{n_{2}}$, there exists a positive constant $L$ such that

  $\displaystyle|f(x_{1},y_{1})-f(x_{2},y_{2})|\vee|g(x_{1},y_{1})-g(x_{2},y_{2})|\leq L(|x_{1}-x_{2}|+|y_{1}-y_{2}|).$

• (F2)

  There exists a constant $\beta>0$ such that for any $x\in\mathbb{R}^{n_{1}}$ and $y_{1},y_{2}\in\mathbb{R}^{n_{2}}$,

  $\displaystyle 2(y_{1}-y_{2})^{T}\big{(}f(x,y_{1})-f(x,y_{2})\big{)}+|g(x,y_{1})-g(x,y_{2})|^{2}\leq-\beta|y_{1}-y_{2}|^{2}.$

• $({\bf F3})$

  For some fixed $k\geq 2$, there exist constants $\alpha_{k}>0$ and $L_{k}>0$ such that for any $x\in\mathbb{R}^{n_{1}}$, $y\in\mathbb{R}^{n_{2}}$,

  $\displaystyle y^{T}f(x,y)+\frac{k-1}{2}|g(x,y)|^{2}\leq-\alpha_{k}|y|^{2}+L_{k}(1+|x|^{2}).$

An averaged equation (1.4) derived from the averaging principle, which provides a substantial simplification of original system (1.1). Then the numerical scheme for system (1.1) can be proposed by solving the averaged equation (1.4) numerically. For this purpose, we cite some known results on the averaging principle firstly.

With the help of the strong averaging principle, this section is devoted to constructing an easily implementable multiscale numerical scheme for the slow component of original SFSDE (1.1). One notices from $({\bf S4})$ that for any $y\in\mathbb{R}^{n_{2}}$,

Then for any $u\geq 1$ and $x\in\mathbb{R}^{n_{1}}$ with $|x|\leq u$

where $\varphi(u)=1+u^{(\theta_{3}\vee\theta_{4}-1)}$, and $\theta_{3},\theta_{4}$ are given in $({\bf S4})$. Then for any step size $\Delta_{1}\in(0,1]$, define

where $x/|x|={\bf 0}\in\mathbb{R}^{n_{1}}$ if $x={\bf 0}$, and $K$ is a constant satisfying $K\geq~{}1+|x_{0}|^{(\theta_{3}\vee\theta_{4}-1)}$. Clearly, for any $x\in\mathbb{R}^{n_{1}}$,

Moreover, under (S4), (F1)-(F3) with $k\geq\theta_{4}$ by the definition (1.3) we derive from the above inequality and (2.1) that

where $\mu^{x^{*}}$ is the unique invariant probability measure of the frozen equation (1.2) with the fixed parameter $x^{*}$, and the last step used the increasing of $\varphi$.

Because the analytical form of $\bar{b}(x^{*})$ is unobtainable, using the ergodicity of the frozen equation (1.2), we approximate $\bar{b}(x^{*})$ by the time average of $b(x^{*},\cdot)$ with respect to the numerical solution of the frozen equation (1.2) with fixed parameter $x^{*}$. For convenience, for an integer $M>0$, we introduce an average function

where $h=\{h_{m}\}_{m=1}^{\infty}$ is a $\mathbb{R}^{n_{2}}$-valued sequence. Within the framework of HMM, we design an easily implementable multiscale numerical scheme involving a macro solver and a micro solver as well as an estimator. For clarity, we illustrate it as follows. Let $\Delta_{1}$ and $\Delta_{2}$ denote macro time step size and micro time step size, respectively.

• (1)

  Macro solver: For the known $\!X_{n}$, since the drift coefficient $\bar{b}$ of the averaged equation may be sup-linear, the truncated EM scheme is selected as macro solver to make a macro step and get $X_{n+1}$. Then we have

  $$X_{n+1}=X_{n}+B_{n}\Delta_{1}+\sigma(X_{n})\Delta W^{1}_{n},$$

  where $B_{n}$ ia an approximation of $\bar{b}(X^{*}_{n})$ that we obtain in third step, and $n_{1}(u):=\lfloor u/\Delta_{1}\rfloor$ for any $u\geq 0$ with $\lfloor t/\delta\rfloor$ the integer part of $t/\delta$, and $\Delta W^{1}_{n}=W^{1}((n+1)\Delta_{1})-W^{1}(n\Delta_{1})$.

• (2)

  Micro solver: To obtain $B_{n}$ at each macro time step, for the known $X_{n}\in\mathbb{R}^{n_{1}}$, use the EM method to solve the frozen equation (1.2) with parameter $x=X^{*}_{n}$ fixed. Therefore, the micro solver is given by

  $$\!\!\!\begin{cases}Y^{X^{*}_{n}}_{0}=y_{0},\\ Y^{X^{*}_{n}}_{m+1}=Y^{X^{*}_{n}}_{m}\!+f(X^{*}_{n},Y^{X^{*}_{n}}_{m})\Delta_{2}\!+g(X^{*}_{n},Y^{X^{*}_{n}}_{m})\Delta W^{2}_{n,m},~{}~{}m=0,1,\cdot\cdot\cdot,\end{cases}$$

  where $\{W^{2}_{n}(\cdot)\}_{n\geq 0}$ is a mutually independent Brownian motion sequence and also independent of $W^{1}(t)$, and $\Delta W^{2}_{n,m}=W^{2}_{n}((m+1)\Delta_{2})-W^{2}_{n}(m\Delta_{2})$.

• (3)

  Estimator: For the known $X_{n}$ and $Y^{X^{*}_{n}}:=\{Y^{X_{n}^{*}}_{m}\}_{m\geq 1}$, let

  $$B_{n}=B_{M}(X^{*}_{n},Y^{X^{*}_{n}})$$

  as an approximation of $\bar{b}(X^{*}_{n})$, where $B_{M}(\cdot,\cdot)$ is defined by (3.4) and $M$ denotes the number of micro time steps used for this approximation.

Overall, for any given $\Delta_{1},\Delta_{2}\in(0,1]$ and integer $M\geq 1$ define the multiscale TEM scheme (MTEM) as follows: for any $n\geq 0$,

By this scheme we define the continuous approximation processes

Note that $\bar{X}(n\Delta_{1})={X}(n\Delta_{1})=X_{n}$, that is, $\bar{X}(t)$ and $X(t)$ coincide with the discrete solution at the grid points, respectively.

In order to better study the strong convergence of the MTEM scheme, we need to study some important properties of the averaged coefficient $\bar{b}(x)$ and its estimator $B_{M}(x,Y^{x}_{n})$ (defined in later) in advance. In this section, we mainly provide some pre-estimates for $\bar{b}(x)$ and $B_{M}(x,Y^{x}_{n})$.

By virtue of Lemma 2.1, we show that the drift term $\bar{b}$ of the averaged equation (1.4) inherits the local Lipschitz continuity. ∎