Semi-implicit Euler-Maruyama method for non-linear time-changed stochastic differential equations

TL;DR

This paper introduces and analyzes a semi-implicit Euler-Maruyama method for approximating non-linear time-changed stochastic differential equations, proving convergence, stability, and demonstrating effectiveness through simulations.

Contribution

The paper develops a semi-implicit EM scheme for complex SDEs with super-linear drift and time-change, establishing convergence, stability, and preservation of asymptotic properties.

Findings

Proved strong convergence of the method.

Established mean square polynomial stability.

Numerical simulations confirm theoretical results.

Abstract

The semi-implicit Euler-Maruyama (EM) method is investigated to approximate a class of time-changed stochastic differential equations, whose drift coefficient can grow super-linearly and diffusion coefficient obeys the global Lipschitz condition. The strong convergence of the semi-implicit EM is proved and the convergence rate is discussed. When the Bernstein function of the inverse subordinator (time-change) is regularly varying at zero, we establish the mean square polynomial stability of the underlying equations. In addition, the numerical method is proved to be able to preserve such an asymptotic property. Numerical simulations are presented to demonstrate the theoretical results.

Equations

$$\mathrm{d}X(t)=f(E(t),X(t))\mathrm{d}E(t)+g(E(t),X(t))\mathrm{d}B(E(t)).$$

$$\mathrm{d}X(t)=\left(X(t)-X^3(t)\right)\mathrm{d}E(t)+X(t)\mathrm{d}B(E(t)),$$

$$\mathbb{E}{\mathrm{e}}^{-rD(t)}={\mathrm{e}}^{-t\varphi (r)},r>0,t\ge 0,$$

$$\varphi (r)=\vartheta r+{\int }_{(0,\infty )}\left(1-{\mathrm{e}}^{-rx}\right)\nu (\mathrm{d}x),r>0,$$

$$E(t):=\inf \{s\ge 0;D(s)>t\},t\ge 0.$$

$$\mathrm{d}X(t)=f(E(t),X(t))\mathrm{d}E(t)+g(E(t),X(t))\mathrm{d}B(E(t)),t\in [0,T],$$

$$⟨x-y,f(t,x)-f(t,y)⟩\le K_1|x-y{|}^2.$$

$${\left|f(t,x)-f(s,x)\right|}^2\le K_2\left(1+|x{|}^a\right)|t-s{|}^{\gamma }$$

$${\left|g(t,x)-g(s,x)\right|}^2\le K_2\left(1+|x{|}^2\right)|t-s{|}^{\gamma }.$$

$${\left|f(t,x)-f(t,y)\right|}^2\le K_3\left(1+|x{|}^b+|y{|}^b\right)|x-y{|}^2$$

$${\left|g(t,x)-g(t,y)\right|}^2\le K_3|x-y{|}^2.$$

$$⟨x,f(t,x)⟩+\frac{p-1}{2}|g(t,x){|}^2\le K_4(1+|x{|}^2).$$

$$|f(t,x){|}^2\le K_5(1+|x{|}^a)$$

$$|g(t,x){|}^2\le K_5(1+|x{|}^2)$$

$$\mathrm{d}Y(t)=f(t,Y(t))\mathrm{d}t+g(t,Y(t))\mathrm{d}B(t),Y(0)=X(0),$$

$$y_{i+1}=y_i+f(t_{i+1},y_{i+1})h+g(t_i,y_i)\Delta B_i,i\in \mathbb{N},$$

$$T\in [D_h(t_n),D_h(t_{n+1})),$$

$$E_{h}(t)=\big{(}\min\{n;D_{h}(t_{n})>t\}-1\big{)}h,$$

$$E_h(t)=ih,\text{when t∈[Dh(ti),Dh(ti+1))}.$$

$$E(t)-h\le E_h(t)\le E(t)\text{a.s.}$$

$$\mathbb{E}{\mathrm{e}}^{\delta E(t)}\le {\mathrm{e}}^{Ct}\text{for all t≥1}.$$

$$\mathbb{P}\left(E(t)\le s\right)=\mathbb{P}\left(D(s)\ge t\right),t,s\ge 0.$$

$$\mathbb{E}{\mathrm{e}}^{\delta E(t)}$$

$$=1+{\int }_1^{\infty }\mathbb{P}\left(E(t)>\frac{1}{\delta }\log r\right)\mathrm{d}r$$

$$=1+{\int }_1^{\infty }\mathbb{P}\left(D\left(\frac{1}{\delta }\log r\right)<t\right)\mathrm{d}r$$

$$=1+\delta {\int }_0^{\infty }\mathbb{P}\left(D(r)<t\right){\mathrm{e}}^{\delta r}\mathrm{d}r.$$

$$\mathbb{P}\left(D(r)<t\right)$$

$$\le {\mathrm{e}}^{{\varphi }^{-1}(2\delta )t}\mathbb{E}{\mathrm{e}}^{-{\varphi }^{-1}(2\delta )D(r)}$$

$$={\mathrm{e}}^{{\varphi }^{-1}(2\delta )t}{\mathrm{e}}^{-r\varphi ({\varphi }^{-1}(2\delta ))}$$

$$={\mathrm{e}}^{{\varphi }^{-1}(2\delta )t-2\delta r}.$$

$$\mathbb{E}{\mathrm{e}}^{\delta E(t)}\le 1+\delta {\mathrm{e}}^{{\varphi }^{-1}(2\delta )t}{\int }_0^{\infty }{\mathrm{e}}^{-\delta r}\mathrm{d}r=1+{\mathrm{e}}^{{\varphi }^{-1}(2\delta )t},$$

$$\mathbb{E}|Y(t){|}^p\le 2^{\frac{p-2}{2}}\left(1+\mathbb{E}|Y(0){|}^p\right){\mathrm{e}}^{p{K}_{4}t}\text{for all t≥0}.$$

$$\mathbb{E}|Y(t)-Y(s){|}^q\le C|t-s{|}^{q/2}{\mathrm{e}}^{Ct},$$

Full text

Key words: time-changed stochastic differential equations, semi-implicit Euler-Maruyama method, strong convergence, mean square polynomial stability, subordinator.

2010 MSC2010: 60H10, 65C30, 60J60.

1 Introduction

In this paper, we study the numerical approximations to a class of time-changed stochastic differential equations (SDEs) which are of the form

[TABLE]

Here the coefficients $f$ and $g$ satisfy some regularity conditions (to be specified in Section 2), $B(t)$ represents a standard Brownian motion, and $E(t)$ is an independent time-change given by an inverse subordinator. The rigorous mathematical definitions are postponed to Section 2.

Since it is in general impossible to derive the explicit solution to such SDEs, numerical approximations become extremely important when one applies them to model uncertain phenomenon in real life. This paper aims to construct a numerical method for these time-changed SDEs. The strong convergence with the convergence rate and the mean square stability of the numerical method are investigated.

To our best knowledge, [8] is the first paper to study the finite time strong convergence of numerical methods for time-changed SDEs by directly discretizing the equations. In [8], the authors used the duality principle established in [9] to construct the Euler-Maruyama (EM) method. In a very recent work [6], the authors studied the EM method for a larger class of time-changed SDEs without the duality principle. However, both of these two works required the coefficients of the time-changed SDEs to satisfy the global Lipschitz condition. This requirement rules out many interesting SDEs like

[TABLE]

where some cubic term appears in the drift coefficient. Moreover, the EM is proved to be divergent to SDEs with super-linear growing coefficients [5].

To cope with such super-linearity, we propose the semi-implicit EM method to approximate the SDEs driven by time-changed Brownian motions in this paper. It should be noted that the semi-implicit EM (also called the backward Euler method) have been studied for approximating different types of SDEs driven by Brownian motions, see [3, 4, 10, 11, 15, 19, 21, 23] and the references therein.

Stabilities in different senses for SDEs driven by time-changed Brownian motion have been discussed in [24]. See [17, 18] for related results when the driven process is a time-changed Lévy process. As far as we know, however, there is no result concerning the stability analysis for numerical methods for time-changed SDEs.

In the three papers mentioned above, the global Lipschitz condition was required for the coefficients of the equations. In this paper, we study the the mean square stability of the underlying time-changed SDEs, where the global Lipschitz condition on the drift coefficients is not required. Then, we investigate the capability of the semi-implicit EM method to reproduce such a property under the similar condition.

The main contributions of this paper are as follows.

• •

The semi-implicit EM method is proved to be convergent to a class of time-changed SDEs and the convergence rate is explicitly given.

• •

We establish the mean square stability of the underlying time-changed SDEs. In addition, the numerical solution is proved to be able to preserve such a property.

The rest of this paper is organized as follows. Section 2 is devoted to some mathematical preliminaries for the time-changed SDEs to be considered in this paper, and some necessary lemmas. The strong convergence of the numerical method is proved in Subsection 3.1, and the mean square stabilities of both underlying and numerical solutions are shown in Subsection 3.2. In Section 4, we present numerical simulations to demonstrate the theoretical results derived in Section 3.

2 Preliminaries

Throughout this paper, unless otherwise specified, we will use the following notation. Let $|\cdot|$ be the Euclidean norm in $\mathbb{R}^{d}$ and $\langle x,y\rangle$ be the inner product of vectors $x,y\in\mathbb{R}^{d}$. If A is a vector or matrix, its transpose is denoted by $A^{T}$. If A is a matrix, its trace norm is denoted by $|A|=\sqrt{trace(A^{T}A)}$. For two real numbers $u$ and $v$, we use $u\wedge v=\min(u,v)$ and $u\vee v=\max(u,v)$.

Moreover, let $(\Omega,{\cal F},\mathbb{P})$ be a complete probability space with a filtration $\left\{{\cal F}_{t}\right\}_{t\geq 0}$ satisfying the usual conditions (that is, it is right continuous and increasing while ${\cal F}_{0}$ contains all $\mathbb{P}$-null sets). Let $B(t)=(B_{1}(t),B_{2}(t),...,B_{m}(t))^{T}$ be an $m$-dimensional ${\cal F}_{t}$-adapted standard Brownian motion. Let $\mathbb{E}$ denote the expectation under the probability measure $\mathbb{P}$.

Let $D(t)$ be an ${\cal F}_{t}$-adapted subordinator (without killing), i.e. a nondecreasing Lévy process on $[0,\infty)$ starting at $D(0)=0$. The Laplace transform of $D(t)$ is of the form

[TABLE]

where the characteristic (Laplace) exponent $\phi:(0,\infty)\rightarrow(0,\infty)$ is a Bernstein function with $\phi(0+):=\lim_{r\downarrow 0}\phi(r)=0$, i.e. a $C^{\infty}$-function such that $(-1)^{n-1}\phi^{(n)}\geq 0$ for all $n\in\mathbb{N}$. Every such $\phi$ has a unique Lévy–Khintchine representation

[TABLE]

where $\vartheta\geq 0$ is the drift parameter and $\nu$ is a Lévy measure on $(0,\infty)$ satisfying $\int_{(0,\infty)}(1\wedge x)\,\nu(\mathrm{d}x)<\infty$. We will focus on the case that $t\mapsto D(t)$ is a.s. strictly increasing, i.e. $\vartheta>0$ or $\nu(0,\infty)=\infty$; obviously, this is also equivalent to $\phi(\infty):=\lim_{r\rightarrow\infty}\phi(r)=\infty$.

Let $E(t)$ be the (generalized, right-continuous) inverse of $D(t)$, i.e.

[TABLE]

We call $E(t)$ an inverse subordinator associated with the Bernstein function $\phi$. Note that $t\mapsto E(t)$ is a.s. continuous and nondecreasing.

We always assume that $B(t)$ and $D(t)$ are independent. The process $B(E(t))$ is called a time-changed Brownian motion, which is trapped whenever $t\mapsto E(t)$ is constant. We remark that the jumps of $t\mapsto D(t)$ correspond to flat pieces of $t\mapsto E(t)$. Due to these traps, the time-change slows down the original Brownian motion $B(t)$, and $B(E(t))$ is understood as a subdiffusion in the literature (cf. [16, 22]).

Consider the following time-changed SDE

[TABLE]

with $\mathbb{E}|X(0)|^{\gamma}<0$ for any $\gamma\in(0,\infty)$, where $f:[0,\infty)\times\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$ and $g:[0,\infty)\times\mathbb{R}^{d}\rightarrow\mathbb{R}^{d\times m}$ are measurable coefficients. We will need the following assumptions on the drift and diffusion coefficients.

Assumption 2.1

There exists a constant $K_{1}>0$ such that, for all $t\geq 0$ and $x,y\in\mathbb{R}^{d}$,

[TABLE]

Assumption 2.2

There exist constants $K_{2}>0$, $a\geq 2$ and $\gamma\in(0,2]$ such that, for all $t,s\geq 0$ and $x\in\mathbb{R}^{d}$,

[TABLE]

and

[TABLE]

Assumption 2.3

Assume that there exist constants $K_{3}>0$ and $b\geq 0$ such that, for all $t\geq 0$ and $x,y\in\mathbb{R}^{d}$,

[TABLE]

and

[TABLE]

Assumption 2.4

Assume that there exist constant $p\geq 2$ and $K_{4}>0$ such that, for all $t\geq 0$ and $x\in\mathbb{R}^{d}$,

[TABLE]

By Assumption 2.3, we can see that there exists a constant $K_{5}>0$ such that

[TABLE]

and

[TABLE]

for all $t\geq 0$ and $x\in\mathbb{R}^{d}$.

According to the duality principle in [9], the time-changed SDE (2.1) and the classical SDE of Itô type

[TABLE]

have a deep connection. The next lemma states such a relation more precisely, which is borrowed from Theorem 4.2 in [9].

Lemma 2.5

Suppose Assumptions 2.1 to 2.3 hold. If $Y(t)$ is the unique solution to the SDE (2.4), then the time-changed process $Y(E(t))$, which is an ${\cal F}_{E(t)}$-semimartingale, is the unique solution to the time-changed SDE (2.1). On the other hand, if $X(t)$ is the unique solution to the time-changed SDE (2.1), then the process $X(D(t))$, which is an ${\cal F}_{t}$-semimartingale, is the unique solution to the SDE (2.4).

The plan to numerically approximate the time-changed SDE (2.1) in this paper is as follows. Firstly, we construct the numerical method for the SDE (2.4). Secondly, we discretize the inverse subordinator $E(t)$. Then the combination of the numerical solution of the SDE (2.4) and the discretized inverse subordinator is used to approximate the solution to the time-changed SDE (2.1).

The semi-implicit EM method for (2.4) is defined as

[TABLE]

with $y_{0}=Y(0)$, where $\Delta B_{i}$ is the Brownian increment following the normal distribution with the mean 0 and the variance $h>0$ and $t_{i}=ih$.

Note that under Assumption 2.1, the semi-implicit EM method (2.5) is well defined for any $h\in(0,1/K_{1})$ (see for example [15]). To be more precisely, this means that given $y_{i}$ is known a unique $y_{i+1}$ can be found. Throughout the paper, we always assume $h\in(0,1/K_{1})$.

We also define the piecewise continuous numerical solution by $y(t):=y_{i}$ for $t\in[t_{i},t_{i+1})$, $i\in\mathbb{N}$.

We follow the idea in [2] to approximate the inverse subordinator $E(t)$ in a time interval $[0,T]$ for any given $T>0$. Firstly, we simulate the path of $D(t)$ by $D_{h}(t_{i})=D_{h}(t_{i-1})+\Delta_{i}$ with $D_{h}(0)=0$, where $\Delta_{i}$ is independently identically sequence with $\Delta_{i}=D(h)$ in distribution. The procedure is stopped when

[TABLE]

for some $n$. Then the approximate $E_{h}(t)$ to $E(t)$ is generated by

[TABLE]

for $t\in[0,T]$. It is easy to see

[TABLE]

The next lemma will be used as the approximation error of $E_{h}(t)$ to $E(t)$, whose proof can be found in [8, 12].

Lemma 2.6

For any $t>0$,

[TABLE]

The following lemma states that the inverse subordinator $E(t)$ is known to have the finite exponential moment, which was proved in [8, 13]. Here, we give an alternative proof, which can, furthermore, provide an explicit upper bound.

Lemma 2.7

For any $\delta>0$, there exists $C=C(\delta)>0$ such that

[TABLE]

*Proof. * By the definition of $E(t)$, it is clear that

[TABLE]

Note that

[TABLE]

Denote by $\phi^{-1}$ the inverse function of $\phi$. By the Chebyshev inequality,

[TABLE]

Thus, for all $t>0$,

[TABLE]

which immediately implies the assertion.

The following result is taken from [14, Theorem 4.1, p. 59].

Lemma 2.8

Suppose that Assumptions 2.1 to 2.4 hold. Then the solution to (2.4) satisfies

[TABLE]

The next lemma is easy; for the sake of completeness and our readers’ convenience, we give a brief proof.

Lemma 2.9

Suppose that Assumptions 2.1 to 2.4 hold. Then for any $q\in(1,2p/a]$ and $t,s\geq 0$ with $|t-s|\leq 1$,

[TABLE]

where $C>0$ is a constant independent of $t$ and $s$.

*Proof. * For any $0\leq s<t$, we derive from (2.4) that

[TABLE]

By the elementary inequality

[TABLE]

with $n=2$, the Hölder inequality and [14, Theorem 7.1, p. 39], we get

[TABLE]

Combining this with (2.2), (2.3) and Lemma 2.8, we obtain

[TABLE]

where $C>0$ is a generic constant independent of $t$ and $s$ that may change from line to line. This completes the proof.

3 Main results

3.1 Strong convergence

Briefly speaking, the following theorem states the strong convergence with the rate of $(\gamma\wedge 1)/2$ of the semi-implicit EM method, which is not surprising. But to our best knowledge, it seems that no existing result fulfills our needs in this paper. In Theorem 3.1, we need to track the temporal variable $t$ as we will replace it by $E(t)$ in Theorem 3.2. In addition, it seems that no such a result exists on the semi-implicit EM method for non-autonomous SDEs.

Theorem 3.1

Suppose that Assumptions 2.1 to 2.4 hold with $p\geq 2(a\vee b)$ and the step size satisfies $h<1/(2(K_{1}+1))$. Then the semi-implicit EM method (2.5) is convergent to (2.4) with

[TABLE]

where $C$ is a constant independent of $t$ and $h$.

*Proof. * From (2.4) and (2.5), it holds that for $i=1,2,\ldots$,

[TABLE]

Taking square on both sides yields

[TABLE]

where

[TABLE]

and

[TABLE]

To estimate $I_{1}$, we rewrite the integrand of $I_{1}$ into three parts

[TABLE]

Using Assumption 2.1, we obtain

[TABLE]

Applying the elementary inequality

[TABLE]

we have

[TABLE]

By Assumption 2.2, we can see

[TABLE]

Thus,

[TABLE]

Applying the elementary inequality (3.1) and Assumption 2.3 gives

[TABLE]

Combining the upper bound estimates of $I_{11}$, $I_{12}$ and $I_{13}$, we conclude that

[TABLE]

By the Hölder inequality, we find

[TABLE]

Taking expectations on both sides of (3.2) and applying Lemmas 2.8 and 2.9, we obtain

[TABLE]

where (and in what follows) $C$ is a generic constant independent of $i$ and the step size $h$ that may change from line to line.

Next, we bound $I_{2}$. Applying the elementary inequality (3.1) again, we have

[TABLE]

Taking expectation on both sides and using the Itô isometry, it follows that

[TABLE]

Rewriting the integrand of the second term on the right hand side, and using the elementary inequality (2.7) with $n=3$ and $q=2$ and Assumptions 2.2 and 2.3, we can see

[TABLE]

Now applying Lemmas 2.8 and 2.9 gives

[TABLE]

Combining (3.3) and (3.4) yields

[TABLE]

which implies that

[TABLE]

Now summing both sides from [math] to $i-1$ yields

[TABLE]

Due to the fact that $ih=t_{i}\leq\mathrm{e}^{Ct_{i}}$, from combining same terms together on both sides we can derive

[TABLE]

By the discrete version of the Gronwall inequality, we have

[TABLE]

Moveover, when $t\in[t_{i},t_{i+1})$ for some $i=1,2,\ldots$, Lemma 2.9 and (3.5) yield

[TABLE]

Therefore, the proof is completed.

Theorem 3.2

Suppose that Assumptions 2.1 to 2.4 hold with $p>2(a\vee b)$ and the step size satisfies $h<1/(2(K_{1}+1))$. Then the combination of the semi-implicit EM solution, $y(t)$, and the discretized inverse subordinator, $E_{h}(t)$, i.e. $y(E_{h}(t))$, converges strongly to the solution of (2.1) with

[TABLE]

where $C$ is a constant independent of $T$ and $h$.

*Proof. * By Lemma 2.5 and (2.7) with $n=2$ and $q=2$,

[TABLE]

By Lemmas 2.6, 2.7 and 2.9, we can see

[TABLE]

On the other hand, it holds from Lemmas 2.6 and 2.7 and Theorem 3.1 that

[TABLE]

Combining (3.6) and (3.7), we obtain the required assertion.

3.2 Stability

In the section, we always assume the existence and uniqueness of the solutions to (2.1) and (2.4). In fact, Assumptions 2.1 to 2.3 are sufficient to guarantee it, but we do not use them explicitly.

A function $F:(0,\infty)\rightarrow(0,\infty)$ is said to be regularly varying at zero with index $\alpha\in\mathbb{R}$ if for any $c>0$,

[TABLE]

Denote by $\operatorname{RV}_{\alpha}$ the class of all regularly varying functions at [math]. A function $F\in\operatorname{RV}_{0}$ is said to be slowly varying at [math]. It is clear that every $F\in\operatorname{RV}_{\alpha}$ can be rewritten as

[TABLE]

where $\ell$ is a slowly varying function at [math].

In the following, we will assume that the Bernstein function $\phi\in\operatorname{RV}_{\alpha}$ with $\alpha\in(0,1)$. Typical examples are

• •

Let $\phi(r)=r^{\alpha}\log^{\beta}(1+r)$ with $0<\alpha<1$ and $0\leq\beta<1-\alpha$. Then $\phi\in\operatorname{RV}_{\alpha+\beta}$;

• •

Let $\phi(r)=r^{\alpha}\log^{-\beta}(1+r)$ with $0<\beta<\alpha<1$. Then $\phi\in\operatorname{RV}_{\alpha-\beta}$;

• •

Let $\phi(r)=\log\left(1+r^{\alpha}\right)$ with $0<\alpha<1$. Then $\phi\in\operatorname{RV}_{\alpha}$;

• •

Let $\phi(r)=r^{\alpha}(1+r)^{-\alpha}$ with $0<\alpha<1$. Then $\phi\in\operatorname{RV}_{\alpha}$.

We refer the reader to [20, Chapter 16] for more examples of such Bernstein functions.

Lemma 3.3

If the Bernstein function $\phi\in\operatorname{RV}_{\alpha}$ with $\alpha\in(0,1)$, then for any $\lambda>0$

[TABLE]

*Proof. * Denote by $\mathcal{L}_{t}[F(t)]$ the Laplace transform of a function $F(t)$. It follows from [8, (3.10)] that for any $s>0$ and $\lambda>0$,

[TABLE]

Since $\phi\in\operatorname{RV}_{\alpha}$, we get

[TABLE]

where $\ell$ is a slowly varying function at [math]. Combining this with Karamata’s Tauberian theorem (cf. [1, Theorem 1.7.6]), it holds that

[TABLE]

Noting that $t\mapsto\ell(1/t)$ is slowly varying at $\infty$, one has (see [1, Proposition 1.3.6 (i)])

[TABLE]

which, together with (3.8), implies the desired limit.

Theorem 3.4

Assume that the Bernstein function $\phi\in\operatorname{RV}_{\alpha}$ with $\alpha\in(0,1)$, and that there exists a constant $L_{1}>0$ such that

[TABLE]

Then

[TABLE]

In other words, the solution to (2.1) is mean square polynomially stable.

*Proof. * Given (3.9), by [14, Theorem 4.4, p. 130] we know that the solution to (2.4) is mean square exponentially stable

[TABLE]

Using Lemma 2.5, we obtain

[TABLE]

It remains to apply Lemma 3.3 to complete the proof.

Remark 3.5

It is interesting to observe that the time-changed SDEs (2.1) is polynomially stable while the dual SDEs (2.4) is stable in the exponential rate. This may be due to the effect of the time-changed Brownian, which slows down the diffusion.

Now, we present our result about the stability of the semi-implicit EM method.

Theorem 3.6

Assume that the Bernstein function $\phi\in\operatorname{RV}_{\alpha}$ with $\alpha\in(0,1)$, and that there exist positive constants $L_{2}$ and $L_{3}$ with $2L_{2}>L_{3}$ such that

[TABLE]

Then

[TABLE]

That is to say, the numerical solution to (2.1) is mean square polynomially stable.

*Proof. * Assume that (3.10) holds with $2L_{2}>L_{3}$, the standard approach (see for example [15]) gives

[TABLE]

where $L_{4}=(2L_{2}-L_{3})/(1+2L_{2})$. Now, replacing $t_{i}$ by $E_{h}(t)$ and using Lemma 2.6, we have

[TABLE]

Now, the application of Lemma 3.3 yields the desired assertion.

Remark 3.7

It is not hard to see that (3.10) together with $2L_{2}>L_{3}$ indicates (3.9) in Theorem 3.4. Hence, it can be seen from Theorem 3.6 that the semi-implicit EM method can preserve the mean square polynomial stability of the underlying time-changed SDE.

4 Numerical simulations

In this section, we will present two numerical examples. The first example is used to illustrate the strong convergence as well as the convergence rate. The second example demonstrates the mean square stability of the numerical stability. Throughout this section, we focus on the case that $E(t)$ is an inverse $0.9$-stable subordinator with Bernstein function $\phi(r)=r^{0.9}$.

Example 4.1

A one-dimensional nonlinear autonomous time-changed SDE

[TABLE]

is considered.

It is not hard to check that Assumptions 2.1 to 2.3 hold for (4.1). Therefore, by Theorem 3.2 the numerical solution proposed in this paper is strongly convergent to the underlying solution with the rate of $1/2$.

For a given step size $h$, one path of the numerical solution to (4.1) is simulated in the following way.

Step 1. The semi-implicit EM method with the step size $h$ is used to simulate the numerical solution, $y(t)=y_{i}$, when $t\in[ih,(i+1)h)$ for $i=1,2,3,...$, to the duel SDE

[TABLE]

Step 2. One path of the subordinator $D(t)$ is simulated with the same step size $h$. (see for example [7]).

Step 3. The $E_{h}(t)$ is found by using (2.6).

Step 4. The combination, $y(E_{h}(t))$, is used to approximate (4.1).

One path of the $0.9$-stable subordinator $D(t)$ is plotted using $h^{-6}$ in Figure 1(a). The corresponding inverse subordinator $E(t)$ is drawn in Figure 1(b). One path of the numerical solution to Example 4.1 is displayed in Figure 1(c).

Now we illustrate the strong convergence and the convergence rate. Since the explicit form of the true solution to (4.1) is hard to obtain. The numerical solution with a small step size, $h_{0}=10^{-6}$, is regarded as the true solution. The step sizes of $h=10^{-2}$, $10^{-3}$ and $10^{-4}$ are used to calculated the numerical solutions. For the given step size $h$, the $L^{1}$ strong error is calculated by

[TABLE]

Two hundreds ($N=200$) sample paths are used to draw Loglog plot of the $L^{1}$ error against the step sizes in Figure 2. The red solid line is the reference line with the slope of 1/2. It can be seen that the strong convergence rate is approximately 1/2. A simple regression also shows that the rate is 0.4996, which is in line with the theoretical one.

Example 4.2

A one-dimensional nonlinear time-changed SDE

[TABLE]

is considered.

It is not hard to check that (3.9) is satisfied, thus the underlying time-changed SDE is stable in the mean square sense. In addition, (3.10) holds for (4.2) indicates the numerical solution is also mean square stable.

One hundred paths are used to draw the mean square of the numerical solutions from $t=0$ to $t=10$. It is clear in Figure 3(a) that the second moments of the solution tends to [math] as the time $t$ advances, which indicates the numerical solution is mean square stable. In addition, five sample paths are displayed in Figure 3(b).

Acknowledgement

Chang-Song Deng would like to thank the National Natural Science Foundation of China (11401442, 11831015).

Wei Liu would like to thank the National Natural Science Foundation of China (11701378, 11871343), Chenguang Program supported by both Shanghai Education Development Foundation and Shanghai Municipal Education Commission (16CG50), and Shanghai Gaofeng & Gaoyuan Project for University Academic Program Development, for their financial support.