Semi-implicit Euler-Maruyama method for non-linear time-changed stochastic differential equations
Chang-Song Deng Wei Liu
math.NA cs.NA math.PR 60H10, 65C30, 60J60










F
G
P



theoremTheorem[section]
lemma[theorem]Lemma
coro[theorem]Corollary
defn[theorem]Definition
assp[theorem]Assumption
cond[theorem]Condition
expl[theorem]Example
prop[theorem]Proposition
rmk[theorem]Remark
conj[theorem]Conjecture











α
ε ζ η øθ
ϑ κ̨ łλ μ ν
ξ  ρ̊ σ
ϕ φ   ω
√ ı ∀ ȷ∅

β δ υ ≡
∨ ∧


Δ ØΘ ŁΛ
Ξ Σ Ω
∂ ∇ ∃ ×
⋁ ⋀

⊘ ⊕ ⊗ ⊖
∘ ⊙ \ ∞
⋂ ⋃ ± ­∓
⟨ ⟩

Proof. 

trace
diag
  for 
=0.3in=1










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


School of Mathematics and Statistics 
Wuhan University, Wuhan 430072, China.


Wei Liu[Corresponding author. Email: weiliu@shnu.edu.cn, lwbvb@hotmail.com]


Department of Mathematics 
Shanghai Normal University, Shanghai, 200234, China.

    
===================================================================================================================================================================================================================================================================




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.

Key words: time-changed stochastic differential equations, semi-implicit Euler-Maruyama method, strong convergence, mean square polynomial stability, subordinator.
2010 MSC2010: 60H10, 65C30, 60J60.




§ INTRODUCTION


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

    X(t) = f(E(t),X(t))  E(t) + g(E(t),X(t))  B(E(t)).

Here the coefficients f and g satisfy some regularity
conditions (to be specified in Section <ref>),
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 <ref>.
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, <cit.> is the first paper to study the finite time strong convergence of numerical methods for time-changed SDEs by directly discretizing the equations. In <cit.>, the authors used the duality principle established in <cit.> to construct the Euler-Maruyama (EM) method. In a very recent work <cit.>, 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

    X(t) = (X(t) - X^3(t))   E(t)
    + X(t)   B(E(t)),

where some cubic term appears in the drift coefficient. Moreover, the EM is proved to be divergent to SDEs with super-linear growing coefficients <cit.>.
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 <cit.> and the references therein.
Stabilities in different senses for SDEs driven
by time-changed Brownian motion have been discussed in <cit.>. See <cit.> 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 <ref> 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 <ref>, and the mean square stabilities of both underlying and numerical solutions are shown in Subsection <ref>. In Section <ref>, we present numerical simulations to demonstrate the
theoretical results derived
in Section <ref>.



§ PRELIMINARIES
 

Throughout this paper, unless otherwise specified, we will use the following notation.
Let |·| be the Euclidean norm in ^d and x,y be the inner product of vectors x,y∈^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|=√(trace(A^TA)). For two real numbers u and v, we
use u v=min(u,v) and u∨ v=max(u,v).
Moreover, let (Ω, , ) be a complete probability space with a filtration {_t}_t ≥ 0 satisfying the usual conditions (that is, it is right continuous and increasing while _0 contains all -null sets). Let B(t)= (B_1(t), B_2(t), ..., B_m(t))^T be an m-dimensional
_t-adapted standard Brownian motion. Let  denote the expectation under the probability measure .


Let D(t) be an _t-adapted subordinator (without killing), i.e. a nondecreasing Lévy process on [0,∞) starting
at D(0)=0. The Laplace transform of D(t) is of
the form

    ^-rD(t) = ^-t ϕ(r),   r>0, t≥0,

where the characteristic (Laplace)
exponent ϕ:(0,∞)→(0,∞) is a Bernstein
function with ϕ(0+):=lim_r↓0ϕ(r)=0,
i.e. a C^∞-function such
that (-1)^n-1ϕ^(n)≥0 for all n∈.
Every such ϕ has a unique Lévy–Khintchine representation

    ϕ(r)
        =ϑ r+∫_(0,∞)(1-^-rx)  ν( x),   r>0,

where ϑ≥0 is the drift parameter and
ν is a Lévy measure on (0,∞)
satisfying ∫_(0,∞)(1∧ x)
 ν( x)<∞. We will focus on the case
that t↦ D(t) is a.s. strictly increasing,
i.e. ϑ>0 or ν(0,∞)=∞;
obviously, this is also equivalent
to ϕ(∞):=lim_r→∞ϕ(r)=∞.


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

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

We call E(t) an inverse subordinator associated with the Bernstein function ϕ. Note that t↦ 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↦ E(t) is constant. We remark that the jumps of t↦ D(t) correspond to flat pieces
of t↦ 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. <cit.>).



Consider the following time-changed SDE

    X(t) = f(E(t),X(t))  E(t)
    + g(E(t),X(t))  B(E(t)),    t∈ [0,T],

with |X(0)|^γ < 0 for any γ∈ (0,∞), where f:[0,∞) ×^d →^d and g:[0,∞) ×^d→^d× m are measurable coefficients. We will need the following assumptions
on the drift and diffusion coefficients.

There exists a constant K_1 > 0 such that,
for all t≥0 and x,y∈^d,

    x - y,   f(t,x) - f(t,y) ≤ K_1
    |x - y|^2.




There exist constants K_2 >0, a ≥ 2 and γ∈ (0,2] such that, for all t,s≥0
and x∈^d,

    | f(t,x) - f(s,x) |^2 ≤ K_2 ( 1 + |x|^a ) |t - s|^γ

and

    | g(t,x) - g(s,x) |^2 ≤ K_2 ( 1 + |x|^2 ) |t - s|^γ.




Assume that there exist constants K_3 > 0 and b≥ 0
such that, for all t≥0 and x,y∈^d,

    | f(t,x) - f(t,y) |^2 ≤ K_3 ( 1 + |x|^b + |y|^b ) |x - y|^2

and

    | g(t,x) - g(t,y) |^2 ≤ K_3  |x - y|^2.




Assume that there exist constant p ≥ 2 and K_4 > 0 such that, for all t≥0 and x∈^d,

    x, f(t,x)  + p-1/2 |g(t,x)|^2 ≤ K_4 (1 + |x|^2).



By Assumption <ref>, we can see that there exists a constant K_5 > 0 such that

    |f(t,x)|^2 ≤ K_5 (1 + |x|^a)

and

    |g(t,x)|^2 ≤ K_5 (1 + |x|^2)

for all t≥0 and x∈^d.



According to the duality principle in <cit.>, the time-changed SDE (<ref>) and the classical SDE of Itô type

    Y(t) = f(t,Y(t))  t + g(t,Y(t))  B(t),   Y(0) = X(0),

have a deep connection. The next lemma states such a relation more precisely, which is borrowed from Theorem 4.2 in <cit.>.

Suppose Assumptions <ref> to <ref> hold. If Y(t) is the unique solution to the SDE (<ref>), then the time-changed process Y(E(t)), which is an _E(t)-semimartingale, is the unique solution to the time-changed SDE (<ref>). On the other hand, if X(t) is the unique solution to the time-changed SDE (<ref>), then the process X(D(t)), which is an _t-semimartingale, is the unique solution to the SDE (<ref>).

The plan to numerically approximate the time-changed SDE (<ref>) in this paper is as follows. Firstly, we construct the numerical method for the SDE (<ref>). Secondly, we discretize the inverse subordinator E(t). Then the combination of the numerical solution of the SDE (<ref>) and the discretized inverse subordinator is used to approximate the solution to the time-changed SDE (<ref>).
The semi-implicit EM method for (<ref>) is defined as

    y_i+1 = y_i + f(t_i+1,y_i+1)h + g(t_i,y_i)
    Δ B_i,   i∈,

with y_0 = Y(0), where Δ 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 <ref>, the semi-implicit EM method (<ref>)
is well defined for any h ∈ (0,1/K_1) (see for example <cit.>). 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 ∈ (0,1/K_1).
We also define the piecewise continuous numerical solution by y(t):= y_i for t ∈ [t_i,t_i+1), i∈.
We follow the idea in <cit.> 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 )+ Δ_i with D_h(0) = 0, where Δ_i is independently identically sequence with Δ_i = D(h) in distribution. The procedure is stopped when

    T ∈ [ D_h(t_n), D_h(t_n+1)),

for some n. Then the approximate E_h(t) to E(t) is generated by

    E_h(t) = (min{n; D_h(t_n) > t} - 1)h,

for t ∈ [0,T]. It is easy to see

    E_h(t) = ih,  when t ∈[ D_h(t_i), D_h(t_i+1)).


The next lemma will be used as the approximation error of E_h(t) to E(t), whose proof can be found in <cit.>.

For any t>0,

    E(t) - h ≤ E_h(t) ≤ E(t)   a.s.



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

For any δ>0, there
exists C=C(δ)>0 such that

    ^δ E(t)≤^Ct  for all t≥1.



By the definition of E(t), it is clear that

    ( E(t) ≤ s ) = ( D(s) ≥ t ),   t,s ≥ 0.

Note that

    ^δ E(t)   =
    ∫_0^∞(^δ E(t) > r )  r
       =1+∫_1^∞( E(t) >1/δlog r
    )  r
       =1+∫_1^∞( D(1/δlog r)<t
    )  r
       =1+δ∫_0^∞( D(r) < t )
    ^δ r  r.

Denote by ϕ^-1 the inverse function
of ϕ. By the Chebyshev inequality,

    ( D(r) < t )
       =(^-ϕ^-1(2δ) D(r) > ^-ϕ^-1(2δ) t)
       ≤^ϕ^-1(2δ) t  ^-ϕ^-1(2δ) D(r)
       =^ϕ^-1(2δ) t^-rϕ(ϕ^-1(2δ))
       =^ϕ^-1(2δ) t-2δ r.

Thus, for all t > 0,

    ^δ E(t)≤1+δ^ϕ^-1(2δ) t∫_0^∞^-δ r  r
        =1+^ϕ^-1(2δ) t,

which immediately implies the assertion. 


The following result is taken from <cit.>.


Suppose that Assumptions <ref> to <ref> hold. Then the solution to (<ref>) satisfies

    |Y(t)|^p ≤ 2^p-2/2(1 +   |Y(0)|^p)^pK_4t  for all t≥0.



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


Suppose that Assumptions <ref> to <ref> hold. Then for any q∈(1,2p/a] and
t,s≥0 with |t - s|≤1,

    |Y(t) - Y(s)|^q ≤ C |t - s|^q/2^Ct,

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



For any 0≤ s < t, we derive from (<ref>) that

    Y(t) - Y(s) = ∫_s^t f(r,Y(r))  r + ∫_s^t g(r,Y(r))  B(r).

By the elementary inequality

    |∑_i=1^nu_i|^q
        ≤ n^q-1∑_i=1^n|u_i|^q,
           u_i∈^d

with n=2, the Hölder inequality
and <cit.>, we get

    |Y(t) - Y(s)|^q ≤    |2(t-s)|^q-1∫_s^t | f(r,Y(r)) |^q  r
       + 2^q/2 - 1 |q(q-1)|^q/2 |t - s|^q/2-1∫_s^t | g(r,Y(r)) |^q  r

Combining this with (<ref>),
(<ref>) and Lemma <ref>,
we obtain

    |Y(t) - Y(s)|^q    ≤ C( |t-s|^q + |t - s|^q ^Ct + |t - s|^q/2 +  |t - s|^q/2^Ct)
       ≤ C |t - s|^q/2^Ct,

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





§ MAIN RESULTS
 



 §.§ Strong convergence
 
Briefly speaking, the following theorem states the strong convergence with the rate of (γ∧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 <ref>, we need to track the temporal variable t as we will replace it by E(t) in Theorem <ref>. In addition, it seems that no such a result exists on the semi-implicit EM method for non-autonomous SDEs.

Suppose that Assumptions <ref> to <ref> hold with p ≥ 2 (a∨ b) and the step size satisfies h <1/(2(K_1+1)). Then the semi-implicit EM method (<ref>) is convergent
to (<ref>) with

    | Y(t) - y(t) |^2 ≤ C h^γ∧1 ^C t,   t≥0,

where C is a constant independent of t and h.


From (<ref>) and (<ref>), it holds that
for i=1,2,…,

    Y(t_i+1) - y_i+1 =   ( Y(t_i) - y_i ) + ∫_t_i^t_i+1( f(s,Y(s)) - f(t_i+1,y_i+1) )  s
       + ∫_t_i^t_i+1( g(s,Y(s)) - g(t_i,y_i) )  B(s).

Taking square on both sides yields

    | Y(t_i+1) - y_i+1|^2 = I_1 + I_2,

where

    I_1   :=  Y(t_i+1) - y_i+1,  ∫_t_i^t_i+1( f(s,Y(s)) - f(t_i+1,y_i+1) )
      s
       =∫_t_i^t_i+1 Y(t_i+1) - y_i+1,    f(s,Y(s)) - f(t_i+1,y_i+1)   s

and

    I_2 :=  Y(t_i+1) - y_i+1,  ( Y(t_i) - y_i ) +  ∫_t_i^t_i+1( g(s,Y(s)) - g(t_i,y_i) )  B(s) .


To estimate I_1, we rewrite the integrand of I_1 into three parts

    Y(t_i+1) - y_i+1,    f(s,Y(s)) - f(t_i+1,y_i+1) 
       = Y(t_i+1) - y_i+1,    f(t_i+1,Y(t_i+1)) - f(t_i+1,y_i+1) 
         +  Y(t_i+1) - y_i+1,    f(s,Y(t_i+1)) - f(t_i+1,Y(t_i+1)) 
         +  Y(t_i+1) - y_i+1,    f(s,Y(s)) - f(s,Y(t_i+1)) 
       =:I_11 + I_12 + I_13.

Using Assumption <ref>, we obtain

    I_11≤ K_1 | Y(t_i+1) - y_i+1|^2.

Applying the elementary inequality

    u,v≤|u|^2+|v|^2/2,
           u,v∈^d,

we have

    I_12≤1/2| Y(t_i+1) - y_i+1|^2 + 1/2| f(s,Y(t_i+1)) - f(t_i+1,Y(t_i+1)) |^2.

By Assumption <ref>, we can see

    | f(s,Y(t_i+1)) - f(t_i+1,Y(t_i+1)) |^2 ≤ K_2 ( 1 + |Y(t_i+1)|^a ) |s - t_i+1|^γ.

Thus,

    I_12≤1/2| Y(t_i+1) - y_i+1|^2 + K_2/2( 1 + |Y(t_i+1)|^a ) |s - t_i+1|^γ.

Applying the elementary inequality (<ref>) and Assumption <ref> gives

    I_13   ≤1/2| Y(t_i+1) - y_i+1|^2 + 1/2| f(s,Y(s)) - f(s,Y(t_i+1)) |^2 
       ≤1/2| Y(t_i+1) - y_i+1|^2 + K_3/2( 1 + |Y(s)|^b + |Y(t_i+1)|^b ) |Y(s) - Y(t_i+1) |^2.

Combining the upper bound estimates of I_11, I_12 and I_13, we conclude that

    I_1    ≤∫_t_i^t_i+1(  (K_1 + 1) | Y(t_i+1) - y_i+1|^2 + K_2/2( 1 + |Y(t_i+1)|^a ) |s - t_i+1|^γ
         + K_3/2( 1 + |Y(s)|^b + |Y(t_i+1)|^b ) |Y(s) - Y(t_i+1) |^2 )  s.

By the Hölder inequality, we find

    ( ( 1 + |Y(s)|^b + |Y(t_i+1)|^b  )|Y(s) - Y(t_i+1) |^2 )
               ≤(( 1 + |Y(s)|^b + |Y(t_i+1)|^b  )^2)^1/2(|Y(s) - Y(t_i+1) |^4 )^1/2.

Taking expectations on both sides of (<ref>) and applying Lemmas <ref> and <ref>, we obtain

    I_1    ≤ (K_1 + 1)h | Y(t_i+1) - y_i+1|^2 + C h^γ + 1 + C h^γ + 1 ^C t_i+1 + C h^2 ^Ct_i+1
       ≤ (K_1 + 1) h | Y(t_i+1) - y_i+1|^2 + C h^γ + 1 ^Ct_i+1,

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 (<ref>) again, we have

    I_2    ≤1/2| Y(t_i+1) - y_i+1|^2 + 1/2|  ( Y(t_i) - y_i ) +  ∫_t_i^t_i+1( g(s,Y(s)) - g(t_i,y_i) )  B(s)  |^2 
       =:1/2| Y(t_i+1) - y_i+1|^2 + 1/2 I_21.

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

    I_21≤ |Y(t_i) - y_i|^2 +  ∫_t_i^t_i+1| g(s,Y(s)) - g(t_i,y_i) |^2  s.

Rewriting the integrand of the second term on the right hand side, and using the elementary inequality (<ref>)
with n=3 and q=2 and Assumptions <ref> and <ref>, we can see

    | g(s,Y(s)) - g(t_i,y_i) |^2    ≤ 3 ( | g(s,Y(s)) - g(s,Y(t_i)) |^2 + | g(s,Y(t_i)) - g(t_i,Y(t_i)) |^2
         + | g(t_i,Y(t_i)) - g(t_i,y_i) |^2 ) 
       ≤ 3 ( K_3 |Y(s) - Y(t_i)|^2 + K_2 (1 + |Y(t_i)|^2)|s - t_i|^γ + K_3 |Y(t_i) - y_i|^2 ).

Now applying Lemmas <ref> and <ref> gives

    I_2 ≤1/2 | Y(t_i+1) - y_i+1|^2 + 1+3K_3 h/2 | Y(t_i) - y_i |^2 + C h^(1+γ)∧2
     ^C t_i+1.

Combining (<ref>) and (<ref>) yields

    | Y(t_i+1) - y_i+1|^2 ≤   (1/2 + h(K_1 + 1)) | Y(t_i+1) - y_i+1|^2 
       + 1+3K_3 h/2 | Y(t_i) - y_i |^2 + C h^(1+γ)∧2
     ^C t_i+1,

which implies that

    | Y(t_i+1) - y_i+1|^2 ≤1+3K_3h/1 - 2h(K_1 + 1)(| Y(t_i) - y_i |^2 +
    Ch^(1+γ)∧2 ^C t_i).

Now summing both sides from 0 to i-1 yields

    ∑_l=1^i | Y(t_l) - y_l|^2 ≤1+3K_3h/1 - 2h(K_1 + 1)(∑_l=0^i-1| Y(t_l) - y_l|^2 + iC h^(1+γ)∧2 ^C t_i).

Due to the fact that ih = t_i ≤^C t_i, from combining same terms together on both sides we can derive

    | Y(t_i) - y_i |^2 ≤h(3K_3 + 2K_1 +2)/1 - 2h(K_1 + 1)∑_l=0^i-1| Y(t_l) - y_l|^2 + C h^γ∧1
     ^C t_i.

By the discrete version of the Gronwall inequality, we have

    | Y(t_i) - y_i |^2 ≤ C h^γ∧1 ^C t_i.

Moveover, when t ∈ [t_i,t_i+1)
for some i=1,2,…, Lemma <ref> and (<ref>) yield

    | Y(t) - y(t) |^2     =| Y(t) - y_i |^2
       ≤2 | Y(t) - Y(t_i) |^2 + 2 | Y(t_i) - y_i |^2 
       ≤ Ch ^Ct+C h^γ∧1 ^Ct_i
       ≤ C h^γ∧1 ^C t.

Therefore, the proof is completed. 




Suppose that Assumptions <ref> to <ref> hold with p > 2 (a∨ 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 (<ref>) with

    | X(T) - y(E_h(T))|^2 ≤ C h^γ∧1^C T,

where C is a constant independent of T and h.


By Lemma <ref> and
(<ref>) with n=2 and q=2,

    | X(T) - y(E_h(T))|^2    = | Y(E(T)) - y(E_h(T))|^2 
       ≤  2 | Y(E(T)) - Y(E_h(T))|^2 + 2 | Y(E_h(T) - y(E_h(T))|^2.

By Lemmas <ref>, <ref> and <ref>,  we can see

    | Y(E(T)) - Y(E_h(T))|^2 ≤ C h  ^C E(T)≤ C h^C(T∨1).

On the other hand, it holds from Lemmas <ref>
and <ref> and Theorem <ref> that

    | Y(E_h(T) - y(E_h(T))|^2 ≤ C h^γ∧1  ^C E_h(T)≤ C h^γ∧1  ^C E(T)≤ C h^γ∧1^C(T∨1).

Combining (<ref>) and (<ref>), we obtain the required assertion.  



 §.§ Stability
 

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




A function F:(0,∞)→(0,∞) is said
to be regularly varying at zero with
index α∈ if for any c>0,

    lim_s↓0F(cs)/F(s)
        =c^α.

Denote by RV_α the class of
all regularly varying functions at 0.
A function F∈RV_0
is said to be slowly varying at 0. It is clear
that every F∈RV_α
can be rewritten as

    F(s)=s^αℓ(s),

where ℓ is a slowly varying function at 0.

In the following, we will assume that the Bernstein
function ϕ∈RV_α with
α∈(0,1). Typical examples are

  
  * Let ϕ(r)=r^αlog^β(1+r)
     with 0<α<1 and 0≤β<1-α.
     Then ϕ∈RV_α+β;


  
  * Let ϕ(r)=r^αlog^-β(1+r)
     with 0<β<α<1.
     Then ϕ∈RV_α-β;

  
  * Let ϕ(r)=log(1+r^α)
     with 0<α<1.
     Then ϕ∈RV_α;

  
  * Let ϕ(r)=r^α(1+r)^-α
     with 0<α<1.
     Then ϕ∈RV_α.

We refer the reader to <cit.> for
more examples of such Bernstein functions.




    If the Bernstein function
    ϕ∈RV_α with α∈(0,1), then for any λ>0
    
    lim_t→∞log ^-λ E(t)/log t
            =-α.




    Denote by ℒ_t[F(t)] the Laplace
    transform of a function F(t).
    It follows from <cit.> that
    for any s>0 and λ>0,
    
    ℒ_t[ ^-λ E(t)](s)
            =ϕ(s)/s[ϕ(s)+λ].

    Since ϕ∈RV_α, we get
    
    sℒ_t[ ^-λ E(t)](s)
            =ϕ(s)/ϕ(s)+λ  ∼  ϕ(s)/λ
            =1/λ s^αℓ(s),
               s↓0,

    where ℓ is a slowly varying
    function at 0. Combining this with Karamata's Tauberian
    theorem (cf. <cit.>), it holds that
    
    ^-λ E(t)  ∼  1/λΓ(1-α) t^-αℓ(1/t),
               t→∞.

    Noting that t↦ℓ(1/t) is slowly
    varying at ∞,
    one has (see <cit.>)
    
    lim_t→∞ℓ(1/t)/log t
            =0,

    which, together with (<ref>), implies
    the desired limit.  



Assume that the Bernstein function ϕ∈RV_α with α∈(0,1), and that there exists a
constant L_1 > 0 such that

    x, f(t,x)  + 1/2 |g(t,x)|^2 ≤
    -L_1 |x|^2,   (x,t) ∈^d × [0,∞).

Then

    lim sup_t →∞log |X(t)|^2/log t≤ -α.

In other words, the solution to (<ref>) is mean square polynomially stable.


Given (<ref>),
by <cit.> we know
that the solution to (<ref>) is mean
square exponentially stable

    |Y(t)|^2 ≤^-L_1 t  |Y(0)|^2.

Using Lemma <ref>, we obtain

    |X(t)|^2 =  |Y(E(t))|^2 ≤ ^-L_1 E(t) · |Y(0)|^2.

It remains to apply Lemma <ref> to
complete the proof. 



It is interesting to observe that the time-changed SDEs (<ref>) is polynomially stable while the dual SDEs (<ref>) 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.

Assume that the Bernstein function ϕ∈RV_α with α∈(0,1), and that there exist positive constants L_2 and L_3 with 2L_2 > L_3 such that

    x, f(t,x) ≤ - L_2 |x|^2   and   |g(t,x)|^2 ≤ L_3 |x|^2,    (x,t) ∈^d × [0,∞).

Then

    lim sup_t →∞log |Y(E_h(t))|^2 /log t≤ -α.

That is to say, the numerical solution to (<ref>) is mean square polynomially stable.


Assume that (<ref>) holds with 2L_2 > L_3, the standard approach (see for example <cit.>) gives

    |Y(t_i)|^2 ≤^-L_4 t_i |Y(0)|^2,

where L_4 = (2L_2 - L_3)/(1 + 2 L_2). Now,
replacing t_i by E_h(t) and using
Lemma <ref>, we have

    |Y(E_h(t))|^2 ≤|Y(0)|^2 · ^-L_4 E_h(t)≤|Y(0)|^2 ·^L_4h · ^-L_4 E(t).

Now, the application of Lemma <ref> yields the desired assertion. 

It is not hard to see that (<ref>) together with 2L_2 > L_3 indicates (<ref>) in Theorem <ref>. Hence, it can be seen from Theorem <ref> that the semi-implicit EM method can preserve the mean square polynomial stability of the underlying time-changed SDE.




§ 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 ϕ(r)=r^0.9.

A one-dimensional nonlinear autonomous time-changed SDE

    X(t) = (X(t) - X^3(t))  E(t) + X(t)   B(E(t)),   with X(0)=1,

is considered.

It is not hard to check that Assumptions <ref> to <ref> hold for (<ref>). Therefore, by Theorem <ref> 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 (<ref>) 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 ∈ [ih, (i+1)h) for i=1,2,3,..., to the duel SDE

    Y(t) = (Y(t) - Y^3(t))  t + Y(t)  B(t),   with Y(0)=1.

Step 2. One path of the subordinator D(t) is simulated with the same step size h. (see for example <cit.>).
Step 3. The E_h(t) is found by using (<ref>).
Step 4. The combination, y(E_h(t)), is used to approximate (<ref>).
One path of the 0.9-stable subordinator D(t) is plotted using h^-6 in Figure <ref>. The corresponding inverse subordinator E(t) is drawn in Figure <ref>. One path of the numerical solution to Example <ref> is displayed in Figure <ref>.



Now we illustrate the strong convergence and the convergence rate. Since the explicit form of the true solution to (<ref>) 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

    1/N∑_i=1^N| y_i(E_h_0(T)) - y_i(E_h(T)) |.

Two hundreds (N=200) sample paths are used to draw Loglog plot of the L^1 error against the step sizes in Figure <ref>. 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.



A one-dimensional nonlinear time-changed SDE

    X(t) = (-X(t) - X^3(t))  E(t) + X(t)   B(E(t)),   with X(0)=5,

is considered.

It is not hard to check that (<ref>) is satisfied, thus the underlying time-changed SDE is stable in the mean square sense. In addition, (<ref>) holds for (<ref>) 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 <ref> that the second moments of the solution tends to 0 as the time t advances, which indicates the numerical solution is mean square stable. In addition, five sample paths are displayed in Figure <ref>.




§ 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.


99

BGT87
N. H. Bingham, C. M. Goldie, and J. L. Teugels.
Regular variation, volume 27 of Encyclopedia of
  Mathematics and its Applications.
Cambridge University Press, Cambridge, 1987.

GM2010
J. Gajda and M. Magdziarz.
Fractional Fokker-Planck equation with tempered α-stable
  waiting times: Langevin picture and computer simulation.
Phys. Rev. E (3), 82(1):011117, 6, 2010.

HMY2007a
D. J. Higham, X. Mao, and C. Yuan.
Almost sure and moment exponential stability in the numerical
  simulation of stochastic differential equations.
SIAM J. Numer. Anal., 45(2):592–609 (electronic), 2007.

Hu1996a
Y. Hu.
Semi-implicit Euler-Maruyama scheme for stiff stochastic
  equations.
In Stochastic analysis and related topics, V (Silivri,
  1994), volume 38 of Progr. Probab., pages 183–202. Birkhäuser
  Boston, Boston, MA, 1996.

HJK2011a
M. Hutzenthaler, A. Jentzen, and P. E. Kloeden.
Strong and weak divergence in finite time of Euler's method for
  stochastic differential equations with non-globally Lipschitz continuous
  coefficients.
Proc. R. Soc. Lond. Ser. A Math. Phys. Eng. Sci.,
  467(2130):1563–1576, 2011.

JK2019
S. Jin and K. Kobayashi.
Strong approximation of stochastic differential equations driven by a
  time-changed Brownian motion with time-space-dependent coefficients.
J. Math. Anal. Appl., 476(2):619–636, 2019.

JUMthesis2015
E. Jum.
Numerical Approximation of Stochastic Differential Equations
  Driven by Lévy Motion with Infinitely Many Jumps.
PhD thesis, University of Tennessee - Knoxville, 2015.

JK2016
E. Jum and K. Kobayashi.
A strong and weak approximation scheme for stochastic differential
  equations driven by a time-changed Brownian motion.
Probab. Math. Statist., 36(2):201–220, 2016.

Kobayashi2011
K. Kobayashi.
Stochastic calculus for a time-changed semimartingale and the
  associated stochastic differential equations.
J. Theoret. Probab., 24(3):789–820, 2011.

LMYY2018
X. Li, Q. Ma, H. Yang, and C. Yuan.
The numerical invariant measure of stochastic differential equations
  with Markovian switching.
SIAM J. Numer. Anal., 56(3):1435–1455, 2018.

LM2015
W. Liu and X. Mao.
Numerical stationary distribution and its convergence for nonlinear
  stochastic differential equations.
J. Comput. Appl. Math., 276:16–29, 2015.

Mag2009
M. Magdziarz.
Stochastic representation of subdiffusion processes with
  time-dependent drift.
Stochastic Process. Appl., 119(10):3238–3252, 2009.

MOW2011
M. Magdziarz, S. Orzel, and A. Weron.
Option pricing in subdiffusive Bachelier model.
J. Stat. Phys., 145(1):187–203, 2011.

M2008a
X. Mao.
Stochastic differential equations and applications.
Horwood Publishing Limited, Chichester, second edition, 2008.

MS2013a
X. Mao and L. Szpruch.
Strong convergence rates for backward Euler-Maruyama method for
  non-linear dissipative-type stochastic differential equations with
  super-linear diffusion coefficients.
Stochastics, 85(1):144–171, 2013.

MS2004
M. M. Meerschaert and H.-P. Scheffler.
Limit theorems for continuous-time random walks with infinite mean
  waiting times.
J. Appl. Probab., 41(3):623–638, 2004.

NN2017
E. Nane and Y. Ni.
Stability of the solution of stochastic differential equation driven
  by time-changed Lévy noise.
Proc. Amer. Math. Soc., 145(7):3085–3104, 2017.

NN2018
E. Nane and Y. Ni.
Path stability of stochastic differential equations driven by
  time-changed Lévy noises.
ALEA Lat. Am. J. Probab. Math. Stat., 15(1):479–507, 2018.

NS2014
A. Neuenkirch and L. Szpruch.
First order strong approximations of scalar SDEs defined in a
  domain.
Numer. Math., 128(1):103–136, 2014.

SSV12
R. L. Schilling, R. Song, and Z. Vondraček.
Bernstein functions. Theory and applications, volume 37 of De Gruyter Studies in Mathematics.
Walter de Gruyter & Co., Berlin, second edition, 2012.

Sch1997a
H. Schurz.
Stability, stationarity, and boundedness of some implicit
  numerical methods for stochastic differential equations and applications.
Logos Verlag Berlin, Berlin, 1997.

UHMK2018
S. Umarov, M. Hahn, and K. Kobayashi.
Beyond the triangle: Brownian motion, Ito calculus, and
  Fokker-Planck equation—fractional generalizations.
World Scientific Publishing Co. Pte. Ltd., Hackensack, NJ, 2018.

WW2010
L. Wang and X. Wang.
Convergence of the semi-implicit Euler method for stochastic
  age-dependent population equations with Poisson jumps.
Appl. Math. Model., 34(8):2034–2043, 2010.

Wu2016arXiv1
Q. Wu.
Stability analysis for a class of nonlinear time-changed systems.
Cogent Math., 3:Art. ID 1228273, 10, 2016.



