In this paper, we study the numerical approximations to a class of time-changed stochastic differential equations (SDEs) which are of the form
Here the coefficients and satisfy some regularity conditions (to be specified in Section 2), represents a standard Brownian motion, and 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,  is the first paper to study the finite time strong convergence of numerical methods for time-changed SDEs by directly discretizing the equations. In , the authors used the duality principle established in  to construct the Euler-Maruyama (EM) method. In a very recent work , 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
where some cubic term appears in the drift coefficient. Moreover, the EM is proved to be divergent to SDEs with super-linear growing coefficients .
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 . 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.
Throughout this paper, unless otherwise specified, we will use the following notation. Let be the Euclidean norm in and
be the inner product of vectors. If A is a vector or matrix, its transpose is denoted by . If A is a matrix, its trace norm is denoted by . For two real numbers and , we use and .
be a complete probability space with a filtrationsatisfying the usual conditions (that is, it is right continuous and increasing while contains all -null sets). Let be an -dimensional -adapted standard Brownian motion. Let denote the expectation under the probability measure .
Let be an -adapted subordinator (without killing), i.e. a nondecreasing Lévy process on starting at . The Laplace transform of is of the form
where the characteristic (Laplace) exponent is a Bernstein function with , i.e. a -function such that for all . Every such has a unique Lévy–Khintchine representation
where is the drift parameter and is a Lévy measure on satisfying . We will focus on the case that is a.s. strictly increasing, i.e. or ; obviously, this is also equivalent to .
Let be the (generalized, right-continuous) inverse of , i.e.
We call an inverse subordinator associated with the Bernstein function . Note that is a.s. continuous and nondecreasing.
We always assume that and are independent. The process is called a time-changed Brownian motion, which is trapped whenever is constant. We remark that the jumps of correspond to flat pieces of . Due to these traps, the time-change slows down the original Brownian motion , and is understood as a subdiffusion in the literature (cf. [16, 22]).
Consider the following time-changed SDE
with for any , where and are measurable coefficients. We will need the following assumptions on the drift and diffusion coefficients.
There exists a constant such that, for all and ,
There exist constants , and such that, for all and ,
Assume that there exist constants and such that, for all and ,
Assume that there exist constant and such that, for all and ,
By Assumption 2.3, we can see that there exists a constant such that
for all and .
have a deep connection. The next lemma states such a relation more precisely, which is borrowed from Theorem 4.2 in .
Suppose Assumptions 2.1 to 2.3 hold. If is the unique solution to the SDE (2.4), then the time-changed process , which is an -semimartingale, is the unique solution to the time-changed SDE (2.1). On the other hand, if is the unique solution to the time-changed SDE (2.1), then the process , which is an -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 . 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
with , whereand .
Note that under Assumption 2.1, the semi-implicit EM method (2.5) is well defined for any (see for example ). To be more precisely, this means that given is known a unique can be found. Throughout the paper, we always assume .
We also define the piecewise continuous numerical solution by for , .
We follow the idea in  to approximate the inverse subordinator in a time interval for any given . Firstly, we simulate the path of by with , where is independently identically sequence with in distribution. The procedure is stopped when
for some . Then the approximate to is generated by
for . It is easy to see
For any ,
The following lemma states that the inverse subordinator
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.
For any , there exists such that
Proof. By the definition of , it is clear that
Denote by the inverse function of . By the Chebyshev inequality,
Thus, for all ,
which immediately implies the assertion.
The following result is taken from [14, Theorem 4.1, p. 59].
The next lemma is easy; for the sake of completeness and our readers’ convenience, we give a brief proof.
Proof. For any , we derive from (2.4) that
By the elementary inequality
with , the Hölder inequality and [14, Theorem 7.1, p. 39], we get
where is a generic constant independent of and 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 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 as we will replace it by in Theorem 3.2. In addition, it seems that no such a result exists on the semi-implicit EM method for non-autonomous SDEs.
To estimate, we rewrite the integrand of into three parts
Using Assumption 2.1, we obtain
Applying the elementary inequality
By Assumption 2.2, we can see
Combining the upper bound estimates of , and , we conclude that
By the Hölder inequality, we find
where (and in what follows) is a generic constant independent of and the step size that may change from line to line.
Next, we bound . Applying the elementary inequality (3.1) again, we have
Taking expectation on both sides and using the Itô isometry, it follows that
which implies that
Now summing both sides from to yields
Due to the fact that , from combining same terms together on both sides we can derive
By the discrete version of the Gronwall inequality, we have
Therefore, the proof is completed.
A function is said to be regularly varying at zero with index if for any ,
Denote by the class of all regularly varying functions at . A function is said to be slowly varying at . It is clear that every can be rewritten as
where is a slowly varying function at .
In the following, we will assume that the Bernstein function with . Typical examples are
Let with and . Then ;
Let with . Then ;
Let with . Then ;
Let with . Then .
We refer the reader to [20, Chapter 16] for more examples of such Bernstein functions.
If the Bernstein function with , then for any
Proof. Denote by the Laplace transform of a function . It follows from [8, (3.10)] that for any and ,
Since , we get
where is a slowly varying function at . Combining this with Karamata’s Tauberian theorem (cf. [1, Theorem 1.7.6]), it holds that
Noting that is slowly varying at , one has (see [1, Proposition 1.3.6 (i)])
which, together with (3.8), implies the desired limit.
Assume that the Bernstein function with , and that there exists a constant such that
In other words, the solution to (2.1) is mean square polynomially stable.