- Email: [email protected]
Strong predictor–corrector Euler–Maruyama methods for stochastic differential equations with Markovian switching
Journal of Computational and Applied Mathematics 237 (2013) 5–17
Contents lists available at SciVerse ScienceDirect
Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam
Strong predictor–corrector Euler–Maruyama methods for stochastic differential equations with Markovian switching Haibo Li, Lili Xiao, Jun Ye ∗ Department of Mathematical Sciences, Tsinghua University, Beijing, 100084, China
article
abstract
info
Article history: Received 11 April 2010 Received in revised form 30 June 2012
In this paper numerical methods for solving stochastic differential equations with Markovian switching (SDEwMSs) are developed by pathwise approximation. The proposed family of strong predictor–corrector Euler–Maruyama methods is designed to overcome the propagation of errors during the simulation of an approximate path. This paper not only shows the strong convergence of the numerical solution to the exact solution but also reveals the order of the error under some conditions on the coefficient functions. A natural analogue of the p-stability criterion is studied. Numerical examples are given to illustrate the computational efficiency of the new predictor–corrector Euler–Maruyama approximation. © 2012 Elsevier B.V. All rights reserved.
MSC: 65C05 60H10 Keywords: Strong predictor–corrector Euler–Maruyama methods Markovian switching Numerical solutions
1. Introduction Stochastic differential equations with Markovian switching (SDEwMSs) arise in mathematical models of hybrid systems that possess frequent unpredictable structural changes. One of the distinct features of such systems is that the underlying dynamics are subject to change with respect to certain configurations. Such models have been used with great success in a variety of application areas, including flexible manufacturing systems, electric power networks, risk theory, financial engineering and insurance modeling; we refer the readers to Cheng et al. [1], Ghosh et al. [2], Jobert and Rogers [3], Mao and Yuan [4], Rolski et al. [5], Smith [6], Wu et al. [7], Yang and Yin [8], Zhao et al. [9] and references therein. Generally, although the fundamental theories such as the existence and uniqueness of the solution as well as stability of SDEwMSs have been well studied, most of SDEwMSs cannot be solved analytically. Thus, appropriate numerical approximation methods such as the Euler (or Euler–Maruyama) method are needed to apply SDEwMSs in practice or to study their properties. Yuan and Mao [10] first considered the numerical solutions of the following stochastic differential equation with Markovian switching dy(t ) = f (y(t ), r (t ))dt + g (y(t ), r (t ))dW (t ),
t ≥0
(1.1)
with initial conditions y(0) = y0 ∈ R and r (0) = r0 ∈ S = {1, 2, . . . , N }, f and g are sufficiently smooth so that Eq. (1.1) has a unique solution. Here y(t ) is referred to the state while r (t ) is regarded as the mode. The system will switch from one mode to another in a random way, and the switching between the modes is governed by a Markov chain. They proved the mean-square convergence of the Euler–Maruyama (EM) approximation for this hybrid stochastic systems, and the order of error was also estimated. Yin et al. [11] extended (1.1) to a family of more general jump-diffusions with Markovian switching, d
∗
Corresponding author. Tel.: +86 10 62788974. E-mail address: [email protected] (J. Ye).
0377-0427/$ – see front matter © 2012 Elsevier B.V. All rights reserved. doi:10.1016/j.cam.2012.07.001
6
H. Li et al. / Journal of Computational and Applied Mathematics 237 (2013) 5–17
and proved the numerical solutions based on the finite-difference procedure converge weakly to the desired limit by means of a martingale problem formulation. During recent years, there also exist extensive studies which prove the convergence of the Euler–Maruyama method applied to some stochastic differential equation with some additional features, like including predictor–corrector or linearimplicit methods and including some sort of delay, jumps, Markovian switching or combinations thereof; see for example, Bruti-Liberati and Platen [12,13], Hou et al. [14], Li and Hou [15], Mao and Yuan [4], Rathinasamy and Balachandran [16], Schurz [17–19] among others. The corresponding proof is basically the same each time, the only novelty coming from changing it a bit to deal with the additional feature. It is well known that the Euler–Maruyama method and most other explicit schemes for solving stochastic differential equations (SDEs) work unreliably and sometimes generate large errors; see for instance Klauder and Petersen [20], Petersen [21], Milstein et al. [22]. Implicit and predictor–corrector schemes are designed to achieve improved numerical stability and turn out to be better suited to simulation task. Generally, implicit schemes usually cost significant computational time and are sometimes not reliably accomplished; however, this phenomenon can be avoided when using some appropriate discrete time schemes, including predictor–corrector methods. In Kloeden and Platen [23], predictor–corrector methods have been proposed as weak discrete time approximations for solving SDEs, which can be used in Monte Carlo simulation. For the strong discrete time approximation of solutions of SDEs, a family of predictor–corrector Euler methods has been developed by Bruti-Liberati and Platen [13]. However, there are no strong predictor–corrector methods available for SDEwMSs yet. In this paper, we develop a new family of strong predictor–corrector Euler–Maruyama (PCEM) methods for SDEwMS (1.1), which are shown to converge with strong order 0.5, and demonstrate their performance by considering some examples. The rest of the paper is arranged as follows. In Section 2 we introduce some necessary notations and define a family of strong predictor–corrector Euler–Maruyama approximate solutions to SDEwMSs. In Section 3 we show that the PCEM solutions converge to the exact solution in L2 under the global Lipschitz condition and reveal that the order of convergence is 0.5. In Section 4 we extend the PCEM convergence results to multi-dimensional case under certain conditions. In Section 5 the numerical stability of SDEwMSs will be introduced and discussed. In Section 6 some numerical examples are given and compared for simulated paths with different degrees of implicitness to illustrate the computational efficiency of the predictor–corrector Euler–Maruyama approximation. Finally, some concluding remarks and future works are provided in Section 7. 2. Preliminary and algorithm Let (Ω , F , {Ft }t ≥0 , P) be a complete probability space with a filtration {Ft }t ≥0 satisfying the usual conditions. Suppose that there is a finite set S = {1, 2, . . . , N }, representing the possible regimes of the environment. We work with a finite-time horizon [0, T ] for some T > 0. Consider the dynamic system given by (1.1) with initial values y(0) = y0 ∈ Rd and r (0) = i0 ∈ S, where f (·, ·) : Rd × S → Rd , g (·, ·) : Rd × S → Rd×m , W = {W (t ) = (W 1 (t ), . . . , W m (t ))T , t ≥ 0} is an m-dimensional Ft -adapted Wiener process, and r = {r (t ), t ≥ 0} is a continuous-time Markov chain taking value in a finite state space S with the generator Q = (qij )N ×N given by P {r (t + δ) = j|r (t ) = i} =
qij δ + o(δ), 1 + qii δ + o(δ),
if i ̸= j, if i = j,
provided δ ↓ 0, and
−qii =
qij < +∞.
i̸=j
We assume that W and r are independent. Throughout this paper, we denote by | · | the Euclidean norm for vectors and
∥ · ∥ the trace norm for matrices. 2.1. Existence and uniqueness Under certain conditions we can establish the existence of a pathwise unique solution of (1.1). Here we make the following global Lipschitz (GL) and linear growth (LG) assumptions.
(H 1) GL. There exists a constant L1 > 0, for all (x, i), (y, i) ∈ Rd × S, such that |f (x, i) − f (y, i)|2 + ∥ g (x, i) − g (y, i) ∥2 ≤ L1 |x − y|2 . (H 2) LG. There exists a constant L2 > 0, for all (x, i) ∈ Rd × S, such that |f (x, i)|2 ∨ ∥ g (x, i) ∥2 ≤ L2 (1 + |x|2 ). Remarks 2.1. It is easy to show that if f (·, ·), g (·, ·) satisfy the GL condition, then they also satisfy the LG condition, but, for the convenience of our later statements, we explicitly require it here.
H. Li et al. / Journal of Computational and Applied Mathematics 237 (2013) 5–17
7
The following theorem guarantee the existence and uniqueness of the solution to Eq. (1.1), which are of use in studying some numerical schemes. Theorem 2.1. If f (x, i), g (x, i) satisfy the conditions (H 1), (H 2), and suppose W , r be independent. Then there exists a unique d-dimensional Ft -adapted right-continuous process y = {y(t ), t ≥ 0} with left-hand limits which satisfies Eq. (1.1) such that y(0) = y0 ∈ Rd and r (0) = i0 ∈ S a.s. Proof. See Theorem 3.13 in [4].
2.2. Algorithm Now we turn our attention to the numerical algorithm. For convenience, we first consider one-dimensional SDEwMSs. Given 0 < ∆ < 1 as a step size, we denote by {ti }i≥1 the usual equidistant time discretization of a bounded interval [0, T ], i.e. t0 = 0, ti − ti−1 = ∆, if tn−1 < T ≤ tn then set tn = T . For simplicity, in this paper we denote by nu the largest integer n such that tn ≤ u, for all u ∈ [0, T ]. Denote ∆Wtk = W (tk+1 ) − W (tk ). For a given partition {tk }k≥1 , {r (tk )}k≥1 is a discrete Markov chain with the transition probability matrix (P (i, j))N ×N ; here P (i, j) = P (r (tk+1 ) = j|r (tk ) = i) is the ijth entry of the matrix e(tk+1 −tk )Q . Thus we could use the following recursion procedure to simulate the discrete Markov chain {r (tk )}k≥1 ; suppose r (tk ) = i1 and generate a random number ξ which is uniformly distributed in [0, 1], then we define
r (tk+1 ) =
i 2 ,
if i2 ∈ S − {N } and
N ,
if
i2 −1
P (i1 , j) ≤ ξ <
j =1 N −1
i2
P (i1 , j),
j =1
P (i1 , j) ≤ ξ .
j =1
Repeating this procedure a trajectory of {r (tk )}k≥1 can be simulated. Now we can introduce a new family of strong predictor–corrector Euler–Maruyama (PCEM) methods for SDEwMSs. Given initial values Yt0 = y0 ∈ R and rt0 = i0 ∈ S, the proposed family of strong PCEM is given by the predictor
Yt
k+1
= Ytk + f (Ytk , rtk )∆ + g (Ytk , rtk )∆Wtk ,
(2.1)
and by the corrector Ytk+1 = Ytk + {θ f¯η ( Ytk+1 , rtk ) + (1 − θ )f¯η (Ytk , rtk )}∆ + {ηg ( Ytk+1 , rtk ) + (1 − η)g (Ytk , rtk )}∆Wtk .
(2.2)
Here parameters θ , η ∈ [0, 1] denote the degree of implicitness in the drift and the diffusion coefficients, respectively, and f¯η (x, i) is defined as f¯η (x, i) = f (x, i) − ηg (x, i)
∂ g (x, i) , ∂x
η ∈ [0, 1],
(2.3)
which is called the corrected drift function. This scheme can be written in the form Ytk+1 = Ytk + f (Ytk , rtk )∆ + g (Ytk , rtk )∆Wtk +
4
Rl,k ,
(2.4)
l=1
where R1,k = θ{f ( Ytk+1 , rtk ) − f (Ytk , rtk )}∆,
R2,k = −θ η
(2.5)
g ( Ytk+1 , rtk )
∂ g ( Ytk+1 , rtk ) ∂ g (Ytk , rtk ) − g (Ytk , rtk ) ∆, ∂x ∂x
(2.6)
∂ g (Ytk , rtk ) ∆, ∂x = η g ( Ytk+1 , rtk ) − g (Ytk , rtk ) ∆Wtk .
R3,k = −ηg (Ytk , rtk )
(2.7)
R4,k
(2.8)
For each t ∈ [tk , tk+1 ), let Y¯ (t ) = Ytk ,
r¯ (t ) = rtk ,
R¯ l (t ) =
R l ,k
∆
,
l = 1, 2, 3,
R¯ 4 (t ) =
R4,k
∆Wtk
.
(2.9)
Therefore, we can define the continuous approximation solution Y (t ) on the entire interval [0, T ] by Y (t ) = y0 +
t
f (Y¯ (s), r¯ (s))ds + 0
t
g (Y¯ (s), r¯ (s))dW (s) + 0
3 l =1
t
R¯ l (s)ds + 0
t
R¯ 4 (s)dW (s). 0
(2.10)
8
H. Li et al. / Journal of Computational and Applied Mathematics 237 (2013) 5–17
Note that Y (tk ) = Y¯ (tk ) = Ytk , which means that Y (t ) and Y¯ (t ) coincide with the discrete approximate solution at the gridpoints. Remarks 2.2. The major advantage of the above PCEM approximate schemes is that there are flexible degrees of implicitness parameters θ and η to choose for simulating paths properly. For the case θ = η = 0 one recovers the Euler–Maruyama scheme which is well discussed in [10]. 3. Convergence with the global Lipschitz (GL) and linear growth (LG) conditions In this section, we will prove that the numerical solution Y converges to the exact solution y in L2 as step size ∆ ↓ 0, and the order of convergence is one-half. To begin with, we need the following GL condition and LG condition for f¯η (·, ·).
(H 3)GL: There exists a constant L3 > 0, for all (x, i), (y, i) ∈ R × S, such that |f¯η (x, i) − f¯η (y, i)|2 ≤ L3 |x − y|2 . (H 4)LG: There exists a constant L4 > 0, for all (x, i) ∈ R × S, such that |f¯η (x, i)|2 ≤ L4 (1 + |x|2 ). Remarks 3.1. For the notation convenience, we denote L = max{L1 , L2 , L3 , L4 }. We are now ready to present the key results of this section which are stated as following. Theorem 3.1. Assume the SDEwMS (1.1) defined on (Ω , F , {Ft }t ≥0 , P) satisfying (a) W , r are independent, (b) f (·, ·), g (·, ·), f¯η (·, ·) satisfy conditions (H 1), (H 2), (H 3) and (H 4), with initial conditions y(0) = y0 ∈ R and r (0) = i0 ∈ S; then the unique strong solution y and the numerical solution Y obtained in Section 2.2 satisfying:
E
sup |Y (t ) − y(t )|2
≤ Ce40L(4+T )T ∆ + o(∆),
(3.1)
0≤t ≤T
where C = 240ML2 (4 + T )T + 10ML(4L + T + 4LT )T + max(−qii ) · 40ML(4 + T ), i∈S
is a positive constant, dependent on y0 , Q , L and T , but independent of ∆ with M = (1 + 7|y0 |2 )e14L(2+8L+T +4LT )T .
(3.2)
In order to give the proof of this theorem, we first provide a number of useful lemmas. The first two lemmas show that the continuous approximation has bounded moments in a strong sense. The latter lemmas play an important role in proving the strong convergence result, which mainly refer to Bruti-Liberati and Platen [13]. Lemma 3.1. Under conditions (H 1), (H 2), (H 3) and (H 4), we have
E
sup |Y (t )|
2
≤ M,
(3.3)
0≤t ≤T
where M is given by (3.2). Proof. By Eq. (2.10), the fundamental inequality and Hölder inequality, for any T ′ ∈ [0, T ], we have
E
sup |Y (t )|
2
2
≤ 7|y0 | + 7TE
sup
0≤t ≤T ′
0≤t ≤T ′
+ 7T
3 l=1
E
t
|R¯ l (s)|2 ds + 7E
sup 0≤t ≤T ′
0≤t ≤T
0 t
2 ¯ |f (Y (s), r¯ (s))| ds + 7E sup
0
t 2 g (Y¯ (s), r¯ (s))dW (s) ′ 0
t 2 sup R¯ 4 (s)dW (s) . ′
0 ≤t ≤T
0
(3.4)
H. Li et al. / Journal of Computational and Applied Mathematics 237 (2013) 5–17
9
For the second part of Eq. (3.4), by the LG condition, Eq. (2.9) and Y (tk ) = Y¯ (tk ) = Ytk , we have
E
T′
|f (Y¯ (s), r¯ (s))| ds ≤ L 2
sup 0 ≤t ≤T ′
t
sup |Y¯ (s)|2
1+E
0
du
0≤s≤u
0 T′
≤L
sup |Y (s)|2
1+E
du.
(3.5)
0≤s≤u
0
Similarly, for the third part of Eq. (3.4), by Doob’s martingale inequality,
E
t 2 t 2 sup g (Y¯ (s), r¯ (s))dW (s) ≤ 4E sup |g (Y¯ (s), r¯ (s))| ds ′ ′
0 ≤t ≤T
0 ≤t ≤T
0
0
T′ ≤ 4L
1+E
sup |Y (s)|
2
du.
(3.6)
0≤s≤u
0
Now, we focus our attention on the remaining parts of (3.4). By Eqs. (2.5), (2.9), the GL condition and the LG condition, we have
E
t
sup
0≤t ≤T ′
2 2 ¯ |R1 (s)| ds ≤ θ LE
0
≤ 2θ LE 2
T′
E
0≤tk ≤tnu
0
T′
E
sup
≤ 2θ L (∆ + ∆)E 2
|f (Ytk , rtk )∆| |Ftk + E 2
0≤tk ≤tnu
0
2 2
2 sup Ytk+1 − Ytk du
T′
|g (Ytk , rtk )∆Wtk | |Ftk du
sup 0≤tk ≤tnu
1 + sup |Ytk | Ftnu du . 0≤tk ≤tnu 2
E 0
2
(3.7)
So by the assumption 0 < ∆ < 1, θ ∈ [0, 1] and Y (tk ) = Y¯ (tk ) = Ytk , we have
E
|R¯ 1 (s)| ds ≤ 4L2 2
sup 0≤t ≤T ′
t
T′
sup |Y (s)|2
1+E
0
du.
(3.8)
0≤s≤u
0
With the same kind of the steps above, from Eq. (2.6), we can easily show that
E
|R¯ 2 (s)| ds ≤ 4L2 2
sup 0≤t ≤T ′
t
T′
sup |Y (s)|2
1+E
0
du.
(3.9)
0≤s≤u
0
Second, by the LG condition, η ∈ [0, 1], Eqs. (2.7), (2.9) and Y (tk ) = Y¯ (tk ) = Ytk , we obtain
E
|R¯ 3 (s)| ds ≤ η2 L 2
sup 0≤t ≤T ′
t
0
T′
1+E T′
du
0≤s≤u
0
≤L
sup |Y¯ (s)|2
1+E
sup |Y (s)|2
du.
(3.10)
0≤s≤u
0
Moreover, by Doob’s martingale inequality, Eq. (2.8) and the GL condition, we have
E
t 2 t 2 ¯ ¯ sup R4 (s)dW (s) ≤ 4E sup |R4 (s)| ds ′ ′
0≤t ≤T
0≤t ≤T
0
≤ 4η LE 2
0
T′
E 0
2 sup Ytk+1 − Ytk du .
(3.11)
0≤tk ≤tnu
Thus, with the similar step in Eqs. (3.7), (3.8) and η ∈ [0, 1], we obtain
E
t 2 sup R¯ 4 (s)dW (s) ≤ 16L2 ′
0≤t ≤T
0
T′
1+E
0
sup |Y (s)|2 0≤s≤u
du.
(3.12)
10
H. Li et al. / Journal of Computational and Applied Mathematics 237 (2013) 5–17
Substituting (3.5), (3.6), (3.8), (3.9), (3.10), (3.12) into (3.4), we get
sup |Y (t )|
1+E
2
0≤t ≤T ′
≤ 1 + 7|y0 |2 + 14L(2 + 8L + T + 4LT ) ×
T′
1+E
0
sup |Y (s)|2
du.
(3.13)
0≤s≤u
Therefore, from the Gronwall inequality we obtain
E
sup |Y (t )|2
≤ (1 + 7|y0 |2 )e14L(2+8L+T +4LT )T .
(3.14)
0≤t ≤T
The proof is complete.
Lemma 3.2. Under conditions (H 1), (H 2), (H 3) and (H 4), we have
E
max |Ytk |2
0≤k≤nT
≤ M,
(3.15)
where M is given by (3.2). This is an immediate result of Lemma 3.1, since Y (tk ) = Y¯ (tk ) = Ytk . Lemma 3.3. Under conditions (H 1), (H 2), (H 3) and (H 4), there exists a constant M which is dependent on y0 , L and T , but independent of ∆, such that 4 l=1
2 Rl,k ≤ ML(4L + T + 4LT )T ∆ + o(∆), E max 1≤n≤nT 0≤k≤n−1
(3.16)
where M is given by (3.2). Proof. By the Cauchy–Schwarz inequality and the GL condition, we can obtain from Eq. (2.5) that
2 R1,k E max 1≤n≤nT 0≤k≤n−1 2 = E max θ {f ( Ytk+1 , rtk ) − f (Ytk , rtk )}∆ 1≤n≤nT 0≤k≤n−1 2 2 |θ ∆| ≤ E max f ( Ytk+1 , rtk ) − f (Ytk , rtk ) 1≤n≤nT
0≤k≤n−1
0≤k≤n−1
≤ θ ∆E 2
∆
0≤k≤nT −1
f ( Yt
2 , rtk ) − f (Ytk , rtk ) k+1
0≤k≤nT −1
≤ θ LT ∆E 2
Yt
k+1
2 − Ytk
0≤k≤nT −1
≤ 2θ LT ∆E 2
2 2 E |f (Ytk , rtk )∆| |Ftk + E |g (Ytk , rtk )∆Wtk | |Ftk .
(3.17)
0≤k≤nT −1
Then by using the Cauchy–Schwarz inequality, the Itô’s isometry, the LG condition, Lemma 3.2 and θ ∈ [0, 1], we get
2 tn T E max R1,k ≤ 2θ 2 L2 T (∆ + 1)∆E E 1 + |Ytnz |2 |Ftnz dz 1≤n≤nT 0 0≤k≤n−1 tn T ≤ 2θ 2 L2 T (∆ + 1)∆ E 1 + max |Ytn |2 dz 0
≤ 2ML T ∆ + o(∆). 2
2
0≤n≤nT
(3.18)
H. Li et al. / Journal of Computational and Applied Mathematics 237 (2013) 5–17
11
With the similar steps as in (3.17) and (3.18), by the GL condition and the LG condition, we have
2 E max R2,k ≤ 2ML2 T 2 ∆ + o(∆). 1≤n≤nT 0≤k≤n−1
(3.19)
It is also easy to show by the LG condition and Lemma 3.2 that
2 E max R3,k ≤ MLT 2 ∆ + o(∆). 1≤n≤nT 0≤k≤n−1
(3.20)
By using Doob’s martingale inequality and Eq. (2.8), we have
2 2 tk + 1 E max R4,k = E max g ( Ytk+1 , rtk ) − g (Ytk , rtk ) dW (z ) η 1≤n≤nT 1≤n≤nT t k 0≤k≤n−1 0≤k≤n−1 2 T G(z )dW (z ) , ≤ 4η2 E 0
where
G(z ) =
g ( Ytk+1 , rtk ) − g (Ytk , rtk ) I[tk ,tk+1 ) (z ),
0≤k≤nT −1
and IΛ be the indicator function for set Λ. Obviously, G(z ) ≤
g ( Ytk+1 , rtk ) − g (Ytk , rtk ) I[0,T ] (z ).
max 0≤k≤nT −1
By Itô’s isometry, the GL condition, with similar steps in (3.17), (3.18) and η ∈ [0, 1], we have
2 T E |G(z )|2 dz R4,k ≤ 4η2 E max 1≤n≤nT 0 0≤k≤n−1 T 2 ≤ 4η2 E max g ( Ytk+1 , rtk ) − g (Ytk , rtk ) I[0,T ] (z )dz 0≤k≤nT −1 0 2 ≤ 4η2 LTE max Ytk+1 − Ytk
0≤k≤nT −1
≤ 4ML2 T ∆ + o(∆). Thus, substituting (3.18), (3.19), (3.20), (3.21) into (3.17), the required assertion follows.
(3.21)
Lemma 3.4. Under conditions (H 1), (H 2), (H 3) and (H 4), for any t ∈ [tk , tk+1 ), we have 3 l =1
E sup t ≤s≤t k
s
2 s 2 R4,k dz + E sup dW (z ) ≤ o(∆). ∆ tk ≤s≤t tk ∆Wtk
R l ,k
tk
(3.22)
The proof of this lemma is similar to that in Lemma 3.3. Now we are in a position to prove Theorem 3.1. Proof of Theorem 3.1. From Eq. (2.10), we have Z (t ) = E
sup |Y (s) − y(s)|2
0≤s≤t
s s = E sup (f (Y¯ (z ), r¯ (z )) − f (y(z ), r (z )))dz + (g (Y¯ (z ), r¯ (z )) − g (y(z ), r (z )))dW (z ) 0≤s≤t 0 0 +
3 l =1
s
R¯ l (z )dz + 0
0
s
2 R¯ 4 (z )dW (z ) .
(3.23)
12
H. Li et al. / Journal of Computational and Applied Mathematics 237 (2013) 5–17
Let n = [t /∆], the integer part of t /∆. Then, by the fundamental inequality, the Hölder inequality, Doob’s martingale inequality and Eq. (2.9), we have Z (t ) = E
sup |Y (s) − y(s)|2
0≤s≤t
t |g (Y¯ (z ), r¯ (z )) − g (y(z ), r (z ))|2 dz |f (Y¯ (z ), r¯ (z )) − f (y(z ), r (z ))|2 dz + 40E 0 0 2 s 2 4 3 R l ,k + 10 E max Rl,k + 10 E sup dz 1≤m≤n tk ≤s≤t tk ∆ l =1 0≤k≤m−1 l =1
t
≤ 10TE
+ 10E sup t ≤s≤t k
s tk
R4,k
∆Wtk
2 dW (z ) .
(3.24)
We focus on the last three parts of the right side. From Lemmas 3.3 and 3.4, we have 10
4 l =1
2 s 2 3 R l ,k + 10 E sup R l ,k dz E max 1≤m≤n ∆ tk ≤s≤t t k l =1 0≤k≤m−1
+ 10E sup t ≤s≤t k
s tk
2 dW (z ) ≤ 10ML(4L + T + 4LT )T ∆ + o(∆). ∆Wtk R4,k
(3.25)
Let t ∈ [tk , tk+1 ). From (2.1)–(2.10), it is easy to know that Y¯ (tk ) and I{r (t )̸=r (tk )} are conditionally independent with respect to the σ -algebra generated by r (tk ). So by the same procedures as in Theorem 3.1 in [10] or Theorem 4.1 in [4], we can obtain t
|f (Y¯ (s), r¯ (s)) − f (y(s), r (s))|2 ds ≤ 2L
E 0
E |Y¯ (s) − y(s)|2 ds + max(−qii ) · 4ML∆ + o(∆),
(3.26)
E |Y¯ (s) − y(s)|2 ds + max(−qii ) · 4ML∆ + o(∆).
(3.27)
i∈S
0 t
|b(Y¯ (s), r¯ (s)) − b(y(s), r (s))|2 ds ≤ 2L
E
t
0
t
i∈S
0
Substituting (3.25), (3.26), (3.27) into (3.24) shows that
E
sup |Y (s) − y(s)|2
≤ 20L(4 + T )
0≤s≤t
t
E |Y¯ (s) − y(s)|2 ds 0
+ 10ML(4L + T + 4LT )T + max(−qii ) · 40ML(4 + T ) ∆ + o(∆).
(3.28)
i∈S
Note that E |Y¯ (s) − y(s)|2 ≤ 2E |Y (s) − y(s)|2 + 2E |Y¯ (s) − Y (s)|2 .
(3.29)
Suppose tk ≤ s < tk+1 , by (2.4), (2.10), Lemmas 3.3 and 3.4, we can then show in the same way as above that E |Y (s) − Y¯ (s)|2 ≤ E |Ytk+1 − Ytk |2 ≤ 6ML∆ + o(∆).
(3.30)
Substituting (3.29), (3.30), into (3.28) immediately shows that
E
sup |Y (s) − y(s)|2
0≤s≤t
t ≤ 40L(4 + T ) (E |Y (s) − y(s)|2 + E |Y¯ (s) − Y (s)|2 )ds 0 + 10ML(4L + T + 4LT )T + max(−qii ) · 40ML(4 + T ) ∆ + o(∆) i∈S t 2 ≤ 40L(4 + T ) E sup |Y (s) − y(s)| du 0
0≤s≤u
+ 240ML2 (4 + T )T + 10ML(4L + T + 4LT )T + max(−qii ) · 40ML(4 + T ) ∆ + o(∆). i∈S
(3.31)
H. Li et al. / Journal of Computational and Applied Mathematics 237 (2013) 5–17
13
Therefore, from the Gronwall inequality we obtain that
E
sup |Y (t ) − y(t )|
2
≤ Ce40L(4+T )T ∆ + o(∆),
0 ≤t ≤T
where C = 240ML2 (4 + T )T + 10ML(4L + T + 4LT )T + max(−qii ) · 40ML(4 + T ). i∈S
The proof is complete.
4. The general multi-dimensional case The results derived in Section 3 can be easily generalized to the multi-dimensional case, we just summarize the related numerical schemes and the convergence results in this section. Consider the solution y = {(y1 (t ), . . . , yd (t ))T , t ≥ 0} of the d-dimensional SDEwMS y(t ) = y(0) +
t
f (y(s), r (s))ds +
0
m j=1
t
g j (y(s), r (s))dW j (s), 0
for t ≥ 0. Here y(0) ∈ Rd and r (0) = r0 ∈ S denotes the deterministic initial values, W j = {W j (t ), t ≥ 0}, j ∈ {1, 2, . . . , m} is a Wiener process. r = {r (t ), t ≥ 0} is a Markov chain. The function f : Rd × S → Rd has the kth component f k (·, ·). The function g j : Rd × S → Rd , j ∈ {1, 2, . . . , m} has the kth component g k,j (·, ·). We define the function f¯η for η = {η1 ∈ [0, 1], . . . , ηd ∈ [0, 1]} with the kth component f¯ηk (x, i) = f k (x, i) − ηk
m d
g i ,j 1 ( x , i )
j1 ,j2 =1 i=1
∂ g i,j2 (x, i) ∂ xi
for (x, i) ∈ R × S, which satisfies the following GL condition and the LG condition. (H 3′ )GL. There exists a constant L5 > 0, for all (x, i), (y, i) ∈ Rd × S, such that |f¯η (x, i) − f¯η (y, i)|2 ≤ L5 |x − y|2 . d
(H 4′ )LG. There exists a constant L6 > 0, for all (x, i) ∈ Rd × S, such that |f¯η (x, i)|2 ≤ L6 (1 + |x|2 ). The kth component of the proposed family of strong PCEM schemes is given by the predictor
Ytk
n+1
= Ytkn + f k (Ytn , rtn )∆ +
m
g k,j (Ytn , rtn )∆Wtn , j
j =1
and by the corrector Ytkn+1 = Ytkn + {θk f¯ηk ( Ytn+1 , rtn ) + (1 − θk )f¯ηk (Ytn , rtn )}∆ +
m {ηk g k,j ( Ytn+1 , rtn ) + (1 − ηk )g k,j (Ytn , rtn )}∆Wtn , j =1
for θk , ηk ∈ [0, 1], k ∈ {1, 2, . . . , d}. Hence we can define the approximation solution Y similarly to Eq. (2.10); then we can derive the following theorem analogically. Theorem 4.1. Assume the SDEwMS (1.1) defined on (Ω , F , {Ft }t ≥0 , P ) satisfying (a) W , r are independent, (b) f (·, ·), g (·, ·), f¯η (·, ·) satisfy conditions (H 1), (H 2), (H 3′ ) and (H 4′ ), with initial conditions y(0) = y0 ∈ Rd and r (0) = i0 ∈ S; then the unique strong solution y and the numerical solution Y satisfy:
E
sup |Y (t ) − y(t )|
0 ≤t ≤T
2
≤ C ∆ + o(∆),
where C is a positive constant dependent on y0 , Q , L and T , but independent of ∆. 5. Numerical stability In this section we consider numerical stability issues, extending the analysis by Platen and Shi [24] to the Markovian switching case. When simulating discrete time approximations of solutions of SDEwMSs, numerical stability is clearly as
14
H. Li et al. / Journal of Computational and Applied Mathematics 237 (2013) 5–17
important as numerical efficiency. There have been various efforts made in the literature trying to study numerical stability for a given scheme approximating solutions of SDEs; see, for instance, Bruti-Liberati and Platen [13], Higham [25], Hofmann and Platen [26], Saito and Mitsui [27], Platen and Shi [24], and Schurz [17–19] among others. Generally, for analyzing numerical stability, some specifically designed test equations are necessary to be introduced to allow insight into the stability behavior of the numerical methods. In deterministic theory the test equation is generally scalar and linear, and indicates that the stability theory has relevance to the behavior near equilibrium on a nonlinear problem. In stochastic theory, the situation is more complex, it is not very clear to which extend the study of a scalar linear test problem is relevant to the nonlinear equations. In this section, we only specialize to a linear test problem and show the connections of stability properties between SDEs and SDEwMSs. For the purpose of comparison, we first consider that the test SDE which has been used by Platen and Shi [24] is a linear and scalar Itô-SDE with multiplicative noise defined as dX (t ) =
1−
3 2
α λX (t )dt + α|λ|X (t )dW (t ),
t ≥ 0,
(5.1)
where X (0) = x0 ≥ 0, λ < 0 and α ∈ [0, 1). Eq. (5.1) can be derived by linearizing about the fixed point X (t ) ≡ 1 of the following SDE
3
dX (t ) = − 1 −
2
α λX (t )(1 − X (t ))dt −
α|λ|X (t )(1 − X (t ))dW (t ),
t ≥ 0,
(5.2)
which is a normalized version of a stochastic population model (see Higham [25]). It is worth noting that the test Eq. (5.1) is just a special case of more general bilinear systems discussed in [17] when investigating the mean square stability and A-stability of numerical methods for SDEs. Obviously, Eq. (5.1) has an explicit solution of the form
X (t ) = x0 exp (1 − α)λt +
α|λ|W (t ) ,
t ≥ 0.
(5.3)
As a unified criterion, Platen and Shi [24] proposed the concept of p-stability criterion, which means that a process is p-stable if in the long run its pth moment vanishes. Hence, for p > 0, λ < 0, p-stability in Eq. (5.1) may be characterized by lim E (|X (t )|p ) = 0
t →∞
0≤α<
iff
1 1+
p 2
.
Since for different combinations of values of λ∆, α and p with given time step size ∆, a discrete time approximation Y and the original continuous process X have different stability properties. To explore these differences, the concept of stability region is introduced. The stability region, denoted by Γ , is determined by those triplets (λ∆, α, p) ∈ (−∞, 0) × [0, 1) × (0, ∞) for which the discrete time approximation Y is p-stable with time step size ∆, when applied to the test Eq. (5.1). By defining the random variable
Y (tn+1 ) , Gn+1 (λ∆, α) = Y (tn )
which is called the transfer function of the approximation Y at time tn , Platen and Shi [24] derived the following useful result which can determine the stability regions for given schemes by the following theorem. Theorem 5.1. A discrete time approximation is for given λ < 0, α ∈ [0, 1) and p > 0, p-stable if and only if E ((Gn+1 (λ∆, α))p ) < 1. In the spirit of Platen and Shi [24], the stability regions for a range of schemes of SDEwMSs are discussed in this paper. Here we consider that the test process Xr = {X (t , r (t )), t ≥ 0} satisfies the linear and scalar Itô-SDEwMS with multiplicative noise dX (t ) =
1−
3 2
α(r (t )) λ(r (t ))X (t )dt + α(r (t ))|λ(r (t ))|X (t )dW (t ),
(5.4)
for t ≥ 0, where r = {r (t ), t ≥ 0} is a Markov chain taking values in a finite state space S = {1, 2, . . . , N }, r (0) = i0 ∈ S , X (0) = x0 ≥ 0 and (α(r (t )), λ(r (t ))) ∈ {(α(i), λ(i)), α(i) ∈ [0, 1), λ(i) < 0, i = 1, 2, . . . , N }. It is well known that the explicit solution of the test equation (5.4) is (see Mao and Yuan [4]) X (t ) = x0 exp
t
(1 − α(r (t )))λ(r (t ))dt + 0
t α(r (t ))|λ(r (t ))|dW (t ) , 0
for t ≥ 0. Now, we introduce the following stability criterion.
(5.5)
H. Li et al. / Journal of Computational and Applied Mathematics 237 (2013) 5–17
15
Definition 5.1. Given p > 0, for all y0 ∈ R and r0 ∈ S, a process Yr = {Y (t , r (t )), t ≥ 0} is called state-p-stable if for each i = 1, 2, . . . , N , Yi = {Y (t , i), t ≥ 0} satisfies lim E (|Y (t , i)|p ) = 0.
t →∞
Remarks 5.1. Actually, if the process Yr is state-p-stable, then it must be p-stable, that is lim E (|Y (t , r (t ))|p ) = 0,
t →∞
∀y0 ∈ R, r0 ∈ S ,
but the reverse is not true; see Pang et al. [28]. For the convenience of numerical comparison, we use this state-p-stability concept in this paper. For more details about p-stability of SDEwMSs, we refer the reader to Pang et al. [28] and references therein. Definition 5.2. The state-stability region Γs is determined by those triplets (λ∆, α, p) ∈ (−∞, 0) × [0, 1) × (0, ∞) for which the discrete time approximation Yr is state-p-stable, when applied to the test Eq. (5.4), for each i = 1, 2, . . . , N, (α(i), λ(i)∆, p) ∈ Γs with time step size ∆. Then we can obtain the following conclusion. Theorem 5.2. The state-stability region Γs , which is generated by one algorithm applied to the test equation (5.4), is the same as the stability region Γ , which is generated by the same algorithm applied to the test equation (5.1). Proof. Γ ⊆ Γs is obvious. To prove Γs ⊆ Γ , suppose for all i = 1, 2, . . . , N, (α(i), λ(i)∆, p) ∈ Γs . By Definition 5.2, we can see that Yr is state-p-stable. And by Definition 5.1, we know for each i = 1, 2, . . . , N , Yi = {Y (t , i), t > 0} is p-stable. So (α(i), λ(i)∆, p) ∈ Γ , for all i ∈ {1, 2, . . . , N }. Immediately, we have Γs ⊆ Γ . Remarks 5.2. By the conclusions derived by Platen and Shi [24], we can also see that the PCEM methods are more efficient than the EM method under these conditions in SDEwMSs case when the parameters selection are reasonable. 6. Numerical examples In this section, we discuss two numerical examples to illustrate our theory established in the previous sections. Let us now consider several combinations of parameters θ and η in Eq. (2.2), the names of the methods listed below are similar to those used by Bruti-Liberati and Platen [13] and Platen and Shi [24]. (1) θ = 0, η = 0 (called the EM scheme), (2) θ = 12 , η = 21 (called the symmetric PCEM scheme), (3) (4) (5) (6)
θ θ θ θ
= 12 , η = 0 (called the semi-drift-implicit PCEM scheme), = 1, η = 0 (called the drift-implicit PCEM method), = 0, η = 21 (called the semi-diffusion-implicit PCEM scheme), = 1, η = 1 (called the fully implicit PCEM scheme).
For a given problem, we will compare the simulated paths for these different degrees of implicitness. If these paths differ significantly from each other, then some numerical stability problem is likely to be present and one needs to make an effort in providing extra numerical stability for further research. Example 6.1. Let W be a Wiener process. Let r be a right-continuous Markov chain taking values in S = {1, 2} with generator
Q =
−0.5 0.5
0.5 . −0.5
And W and r are assumed to be independent. Consider a one-dimensional linear SDEwMS dy(t ) = y(t )a(r (t ))ds + y(t )b(r (t ))dW (t ),
t ≥ 0,
(6.1)
where a(i), b(i), for all i ∈ S, are constants, and the initial value y0 ∈ R, r (0) = i0 ∈ S. It is well known that Eq. (6.1) has an explicit solution y(t ) = y0 exp
t
a(r (s))ds + 0
t
b(r (s))dW (s) − 0
1 2
t
b (r (s))ds , 2
t ≥ 0.
(6.2)
0
For simulation reason, it is convenient to transform (6.2) into the following recursion form with y(t0 ) = y0 ,
√
y(tk+1 ) ≈ y(tk ) exp (tk+1 − tk )a(rtk ) + b(rtk ) tk+1 − tk ξk+1 −
1 2
(tk+1 − tk )b2 (rtk ) .
(6.3)
16
H. Li et al. / Journal of Computational and Applied Mathematics 237 (2013) 5–17
Table 1 Estimation of the errors between the numerical and exact solutions.
E suptk ∈[0,10] |Y(θi ,ηi ) (tk ) − y(tk )|2 ∆ \ θi , ηi
0, 0
1 2
,
0.1 0.05 0.01 0.005 0.001
4.0e−03 1.0e−03 4.7e−05 1.4e−05 1.3e−06
5.8e−06 5.0e−07 4.2e−09 7.8e−10 2.4e−11
1 2
1 2
,0
2.1e−04 7.5e−05 1.0e−05 4.9e−06 9.5e−07
1, 0
0,
4.7e−03 1.1e−03 4.7e−05 1.4e−05 1.3e−06
3.6e−03 8.8e−04 3.4e−05 8.6e−06 3.4e−07
1, 0
0,
5.4e+07 1.3e+07 5.2e+05 1.3e+05 5.3e+03
2.9e+07 9.6e+06 4.9e+05 1.3e+05 5.2e+03
1 2
1, 1 5.2e−03 1.2e−03 5.0e−05 1.5e−05 1.3e−06
Table 2 Estimation of the errors between the numerical and exact solutions.
E suptk ∈[0,10] |Y(θi ,ηi ) (tk ) − y(tk )|2 ∆ \ θi , ηi
0, 0
1 2
,
0.1 0.05 0.01 0.005 0.001
3.0e+07 9.9e+06 5.1e+05 1.3e+05 5.5e+03
6.1e+04 5.6e+03 3.8e+01 6.3e+00 1.7e−01
1 2
1 2
,0
1.5e+05 3.0e+04 2.3e+03 9.5e+02 1.6e+02
1 2
1, 1 5.4e+07 1.3e+07 5.2e+05 1.3e+05 5.3e+03
Here each ξk+1 is an independent standard √ normal random variable, so that the increment ∆Wtk = W (tk+1 ) − W (tk ) of the Wiener process W can be generated by tk+1 − tk ξk+1 . Notice that y(tk+1 ) in (6.3) is not the exact value of y(t ) at the division points tk+1 , because r (s) is not necessarily constant on s ∈ [tk , tk+1 ]. However, since P {r (tk+1 ) = i|r (tk ) = i} = 1 + qii (tk+1 − tk ) + o(tk+1 − tk ) → 1 as ∆ ↓ 0, for sufficiently small ∆, it is reasonable to use (6.3) as an approximation of the exact solution of y. Choose initial values y0 = 10, i0 = 1, T is fixed at 10, a(1) = −0.5, a(2) = −1,
b(1) = 0.1 b(2) = 0.1.
By applying the previously described procedure in Section 2, the trajectory of the approximate solution Y with given step size ∆ can be constructed. To carry out the numerical simulation we successively choose the step size ∆ as in the following Table 1, and for each ∆, we repeatedly simulate and compute suptk ∈[0,10] |Y(θi ,ηi ) (tk ) − y(tk )|2 , (i = 1, 2, 3, 4, 5, 6) for 106 times, then calculate the
sample mean E [suptk ∈[0,10] |Y(θi ,ηi ) (tk ) − y(tk )|2 ] (i = 1, 2, 3, 4, 5, 6). The results are listed in the following Table 1. Clearly, we can easily show that Example 6.1 satisfies the conditions of Eq. (5.4), so it will be state-p-stable, and the simulation results listed in Table 1 just illustrate the stability properties in a certain extent. On the one hand, the numerical method reveals that the numerical solution Y defined by the strong PCEM methods converge to the exact solution y in L2 as step size ∆ ↓ 0. On the other hand, when the parameters’ selection is reasonable, the PCEM methods are much more efficient than the EM method, which strongly demonstrate our results. Example 6.2. Consider the same 1-dimensional linear SDEwMS dy(t ) = y(t )a(r (t ))ds + y(t )b(r (t ))dW (t ),
t ≥ 0.
(6.4)
This time, choose initial values y0 = 10, r0 = 1, T = 10, and
a(1) = 0.5, a(2) = 1,
b(1) = 0.1 b(2) = 0.1
By applying the previously described procedure, the trajectory of the approximate solution Y with given step size ∆ can be constructed. To carry out the numerical simulation we successively choose the step size ∆ as in the following Table 2, and for each ∆, we repeatedly simulate and compute suptk ∈[0,10] |Y(θi ,ηi ) (tk ) − y(tk )|2 , (i = 1, 2, 3, 4, 5, 6) for 106 times; then calculate the
sample mean E [suptk ∈[0,10] |Y(θi ,ηi ) (tk ) − y(tk )|2 ] (i = 1, 2, 3, 4, 5, 6). The results are listed in the following Table 2. Since we can easily show that Example 6.2 does not satisfy the conditions of Eq. (5.4), it will not be state-p-stable, even not p-stable, the computer simulation results illustrate this point to some extent. Nevertheless, when the parameters θ and η are selected reasonably, for example, θ = η = 12 in this example, some PCEM methods are much more efficient than the EM method to a certain extent.
H. Li et al. / Journal of Computational and Applied Mathematics 237 (2013) 5–17
17
7. Concluding remarks In this paper, we investigated the PCEM methods based on jump-adapted time discretizations to handle scenario simulation of the exact solutions of SDEwMSs under global Lipschitz conditions. We showed that these numerical schemes give strong L2 -convergence with order 0.5. Furthermore, we established the numerical state p-stability criterion of the PCEM methods for a specifically designed linear test equation. It has been shown that such schemes exhibit good computational efficiency when the parameters are selected. There are several interesting problems that deserve further investigation. For example, it would clearly be of interest to extend the strong convergence theory to the case where coefficients are not globally Lipschitz continuous, it remained an open question whether the approximation using the PCEM methods also converges in the strong sense if the coefficients of the equations are not globally Lipschitz continuous. Another important problem is how to develop and analyze the various numerical stability criteria for such schemes when making decisions about the selection of an appropriate combination of parameters for strong PCEM methods for a given simulation task. Since for each possible combination of these parameters, the schemes are of strong order 0.5. Nevertheless, the numerical results show that the symmetric PCEM scheme exhibits good numerical properties. The most common method for clarifying this issue is likely to study the shapes of stability regions for related approximating schemes. It is worthwhile to compare the observed stability regions with these numerical methods, which are known to show good numerical stability properties. We hope to address at least some of the specific issues in subsequent work. Acknowledgments The author gratefully acknowledges many useful comments and suggestions from the referees and editors. References [1] C.T. Cheng, C.P. Ou, K.W. Chau, Combining a fuzzy optimal model with a genetic algorithm to solve multi-objective rainfall-runoff model calibration, J. Hydrol. 268 (1–4) (2002) 72–86. [2] M.K. Ghosh, A. Arapostathis, S.I. Marcus, Ergodic control of switching diffusions, SIAM J. Control Optim. 35 (6) (1997) 1952–1988. [3] A. Jobert, L.C.G. Rogers, Option pricing with Markov-modulated dynamics, SIAM J. Control Optim. 44 (6) (2006) 2063–2078. [4] X. Mao, C. Yuan, Stochastic Differential Equations with Markovian Switching, Imperial College Press, 2006. [5] T. Rolski, H. Schmidili, V. Schmidt, J. Teugels, Stochastic Processes for Insurance and Finance, Wiley Online Library, 2000. [6] D.R. Smith, Markov-switching and stochastic volatility diffusion models of short-term interest rates, J. Bus. Econom. Statist. 20 (2) (2002) 183–197. [7] C.L. Wu, K.W. Chau, J.S. Huang, Modelling coupled water and heat transport in a soil-mulch-plant-atmosphere continuum (SMPAC) system, Appl. Math. Model. 31 (2) (2007) 152–169. [8] H. Yang, G. Yin, Ruin probability for a model under Markovian switching regime, in: Probability, Finance and Insurance: Proceedings of a Workshop at the University of Hong Kong, Hong Kong, 15–17 July 2002, World Scientific, 2004, pp. 206–217. [9] M.Y. Zhao, C.T. Cheng, K.W. Chau, G. Li, Multiple criteria data envelopment analysis for full ranking units associated to environment impact assessment, Int. J. Environ. Pollut. 28 (3) (2006) 448–464. [10] C. Yuan, X. Mao, Convergence of the Euler–Maruyama method for stochastic differential equations with Markovian switching, Math. Comput. Simulation 64 (2) (2004) 223–235. [11] G. Yin, Q.S. Song, Z. Zhang, Numerical solutions for jump-diffusions with regime switching, Stochastics 77 (1) (2005) 61–79. [12] N. Bruti-Liberati, E. Platen, Strong approximations of stochastic differential equations with jumps, J. Comput. Appl. Math. 205 (2) (2007) 982–1001. [13] N. Bruti-Liberati, E. Platen, Srong predictor–corrector Euler methods for stochastic differetial equations, Stoch. Dyn. 8 (3) (2008) 561–581. [14] Z. Hou, J. Tong, Z. Zhang, Convergence of jump-diffusion non-linear differential equation with phase semi-Markovian switching, Appl. Math. Model. 33 (9) (2009) 3650–3660. [15] R. Li, Y. Huo, Convergence and stability of numerical solutions to SDDEs with Markovian switching, Appl. Math. Comput. 175 (2) (2006) 1080–1091. [16] A. Rathinasamy, K. Balachandran, Mean square stability of semi-implicit Euler method for linear stochastic differential equations with multiple delays and Markovian switching, Appl. Math. Comput. 206 (2) (2008) 968–979. [17] H. Schurz, Asymptotical mean square stability of an equilibrium point of some linear numerical solutions with multiplicative noise, Stoch. Anal. Appl. 14 (3) (1996) 313–353. [18] H. Schurz, Stability, Stationarity, and Boundedness of Some Implicit Numerical Methods for Stochastic Differential Equations and Applications, LogosVerlag, Berlin, 1997. [19] H. Schurz, The invariance of asymptotic laws of stochastic systems under discretization, Z. Angew. Math. Mech. 79 (6) (1999) 375–382. [20] J.R. Klauder, W.P. Petersen, Numerical integration of multiplicative-noise stochastic differential equations, SIAM J. Numer. Anal. 22 (1985) 1153–1166. [21] W.P. Petersen, Some experiments on numerical simulations of stochastic differential equations and a new algorithm, J. Comput. Phys. 113 (1) (1994) 75–81. [22] G.N. Milstein, E. Platen, H. Schurz, Balanced implicit methods for stiff stochastic systems, SIAM J. Numer. Anal. 35 (3) (1998) 1010–1019. [23] P.E. Kloeden, E. Platen, Numerical Solution of Stochastic Differential Equations, Springer-Verlag, New York, 1999. [24] E. Platen, L. Shi, On the numerical stability of simulation methods for SDEs, Research Paper Series, Quantitative Finance Research Centre, University of Technology, Sydney, 2008, 234, pp. 1–25. [25] D.J. Higham, Mean-square and asymptotic stability of numerical methods for stochastic ordinary differential equations, SIAM J. Numer. Anal. 38 (3) (2000) 753–769. [26] N. Hofmann, E. Platen, Stability of weak numerical schemes for stochastic differential equations, Comput. Math. Appl. 28 (10–12) (1994) 45–57. [27] Y. Saito, T. Mitsui, Stability analysis of numerical schemes for stochastic differential equations, SIAM J. Numer. Anal. 33 (6) (1996) 2254–2267. [28] S. Pang, F. Deng, X. Mao, Almost sure and moment exponential stability of Euler–Maruyama discretizations for hybrid stochastic differential equations, J. Comput. Appl. Math. 213 (1) (2008) 127–141.
Contents lists available at SciVerse ScienceDirect
Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam
Strong predictor–corrector Euler–Maruyama methods for stochastic differential equations with Markovian switching Haibo Li, Lili Xiao, Jun Ye ∗ Department of Mathematical Sciences, Tsinghua University, Beijing, 100084, China
article
abstract
info
Article history: Received 11 April 2010 Received in revised form 30 June 2012
In this paper numerical methods for solving stochastic differential equations with Markovian switching (SDEwMSs) are developed by pathwise approximation. The proposed family of strong predictor–corrector Euler–Maruyama methods is designed to overcome the propagation of errors during the simulation of an approximate path. This paper not only shows the strong convergence of the numerical solution to the exact solution but also reveals the order of the error under some conditions on the coefficient functions. A natural analogue of the p-stability criterion is studied. Numerical examples are given to illustrate the computational efficiency of the new predictor–corrector Euler–Maruyama approximation. © 2012 Elsevier B.V. All rights reserved.
MSC: 65C05 60H10 Keywords: Strong predictor–corrector Euler–Maruyama methods Markovian switching Numerical solutions
1. Introduction Stochastic differential equations with Markovian switching (SDEwMSs) arise in mathematical models of hybrid systems that possess frequent unpredictable structural changes. One of the distinct features of such systems is that the underlying dynamics are subject to change with respect to certain configurations. Such models have been used with great success in a variety of application areas, including flexible manufacturing systems, electric power networks, risk theory, financial engineering and insurance modeling; we refer the readers to Cheng et al. [1], Ghosh et al. [2], Jobert and Rogers [3], Mao and Yuan [4], Rolski et al. [5], Smith [6], Wu et al. [7], Yang and Yin [8], Zhao et al. [9] and references therein. Generally, although the fundamental theories such as the existence and uniqueness of the solution as well as stability of SDEwMSs have been well studied, most of SDEwMSs cannot be solved analytically. Thus, appropriate numerical approximation methods such as the Euler (or Euler–Maruyama) method are needed to apply SDEwMSs in practice or to study their properties. Yuan and Mao [10] first considered the numerical solutions of the following stochastic differential equation with Markovian switching dy(t ) = f (y(t ), r (t ))dt + g (y(t ), r (t ))dW (t ),
t ≥0
(1.1)
with initial conditions y(0) = y0 ∈ R and r (0) = r0 ∈ S = {1, 2, . . . , N }, f and g are sufficiently smooth so that Eq. (1.1) has a unique solution. Here y(t ) is referred to the state while r (t ) is regarded as the mode. The system will switch from one mode to another in a random way, and the switching between the modes is governed by a Markov chain. They proved the mean-square convergence of the Euler–Maruyama (EM) approximation for this hybrid stochastic systems, and the order of error was also estimated. Yin et al. [11] extended (1.1) to a family of more general jump-diffusions with Markovian switching, d
∗
Corresponding author. Tel.: +86 10 62788974. E-mail address: [email protected] (J. Ye).
0377-0427/$ – see front matter © 2012 Elsevier B.V. All rights reserved. doi:10.1016/j.cam.2012.07.001
6
H. Li et al. / Journal of Computational and Applied Mathematics 237 (2013) 5–17
and proved the numerical solutions based on the finite-difference procedure converge weakly to the desired limit by means of a martingale problem formulation. During recent years, there also exist extensive studies which prove the convergence of the Euler–Maruyama method applied to some stochastic differential equation with some additional features, like including predictor–corrector or linearimplicit methods and including some sort of delay, jumps, Markovian switching or combinations thereof; see for example, Bruti-Liberati and Platen [12,13], Hou et al. [14], Li and Hou [15], Mao and Yuan [4], Rathinasamy and Balachandran [16], Schurz [17–19] among others. The corresponding proof is basically the same each time, the only novelty coming from changing it a bit to deal with the additional feature. It is well known that the Euler–Maruyama method and most other explicit schemes for solving stochastic differential equations (SDEs) work unreliably and sometimes generate large errors; see for instance Klauder and Petersen [20], Petersen [21], Milstein et al. [22]. Implicit and predictor–corrector schemes are designed to achieve improved numerical stability and turn out to be better suited to simulation task. Generally, implicit schemes usually cost significant computational time and are sometimes not reliably accomplished; however, this phenomenon can be avoided when using some appropriate discrete time schemes, including predictor–corrector methods. In Kloeden and Platen [23], predictor–corrector methods have been proposed as weak discrete time approximations for solving SDEs, which can be used in Monte Carlo simulation. For the strong discrete time approximation of solutions of SDEs, a family of predictor–corrector Euler methods has been developed by Bruti-Liberati and Platen [13]. However, there are no strong predictor–corrector methods available for SDEwMSs yet. In this paper, we develop a new family of strong predictor–corrector Euler–Maruyama (PCEM) methods for SDEwMS (1.1), which are shown to converge with strong order 0.5, and demonstrate their performance by considering some examples. The rest of the paper is arranged as follows. In Section 2 we introduce some necessary notations and define a family of strong predictor–corrector Euler–Maruyama approximate solutions to SDEwMSs. In Section 3 we show that the PCEM solutions converge to the exact solution in L2 under the global Lipschitz condition and reveal that the order of convergence is 0.5. In Section 4 we extend the PCEM convergence results to multi-dimensional case under certain conditions. In Section 5 the numerical stability of SDEwMSs will be introduced and discussed. In Section 6 some numerical examples are given and compared for simulated paths with different degrees of implicitness to illustrate the computational efficiency of the predictor–corrector Euler–Maruyama approximation. Finally, some concluding remarks and future works are provided in Section 7. 2. Preliminary and algorithm Let (Ω , F , {Ft }t ≥0 , P) be a complete probability space with a filtration {Ft }t ≥0 satisfying the usual conditions. Suppose that there is a finite set S = {1, 2, . . . , N }, representing the possible regimes of the environment. We work with a finite-time horizon [0, T ] for some T > 0. Consider the dynamic system given by (1.1) with initial values y(0) = y0 ∈ Rd and r (0) = i0 ∈ S, where f (·, ·) : Rd × S → Rd , g (·, ·) : Rd × S → Rd×m , W = {W (t ) = (W 1 (t ), . . . , W m (t ))T , t ≥ 0} is an m-dimensional Ft -adapted Wiener process, and r = {r (t ), t ≥ 0} is a continuous-time Markov chain taking value in a finite state space S with the generator Q = (qij )N ×N given by P {r (t + δ) = j|r (t ) = i} =
qij δ + o(δ), 1 + qii δ + o(δ),
if i ̸= j, if i = j,
provided δ ↓ 0, and
−qii =
qij < +∞.
i̸=j
We assume that W and r are independent. Throughout this paper, we denote by | · | the Euclidean norm for vectors and
∥ · ∥ the trace norm for matrices. 2.1. Existence and uniqueness Under certain conditions we can establish the existence of a pathwise unique solution of (1.1). Here we make the following global Lipschitz (GL) and linear growth (LG) assumptions.
(H 1) GL. There exists a constant L1 > 0, for all (x, i), (y, i) ∈ Rd × S, such that |f (x, i) − f (y, i)|2 + ∥ g (x, i) − g (y, i) ∥2 ≤ L1 |x − y|2 . (H 2) LG. There exists a constant L2 > 0, for all (x, i) ∈ Rd × S, such that |f (x, i)|2 ∨ ∥ g (x, i) ∥2 ≤ L2 (1 + |x|2 ). Remarks 2.1. It is easy to show that if f (·, ·), g (·, ·) satisfy the GL condition, then they also satisfy the LG condition, but, for the convenience of our later statements, we explicitly require it here.
H. Li et al. / Journal of Computational and Applied Mathematics 237 (2013) 5–17
7
The following theorem guarantee the existence and uniqueness of the solution to Eq. (1.1), which are of use in studying some numerical schemes. Theorem 2.1. If f (x, i), g (x, i) satisfy the conditions (H 1), (H 2), and suppose W , r be independent. Then there exists a unique d-dimensional Ft -adapted right-continuous process y = {y(t ), t ≥ 0} with left-hand limits which satisfies Eq. (1.1) such that y(0) = y0 ∈ Rd and r (0) = i0 ∈ S a.s. Proof. See Theorem 3.13 in [4].
2.2. Algorithm Now we turn our attention to the numerical algorithm. For convenience, we first consider one-dimensional SDEwMSs. Given 0 < ∆ < 1 as a step size, we denote by {ti }i≥1 the usual equidistant time discretization of a bounded interval [0, T ], i.e. t0 = 0, ti − ti−1 = ∆, if tn−1 < T ≤ tn then set tn = T . For simplicity, in this paper we denote by nu the largest integer n such that tn ≤ u, for all u ∈ [0, T ]. Denote ∆Wtk = W (tk+1 ) − W (tk ). For a given partition {tk }k≥1 , {r (tk )}k≥1 is a discrete Markov chain with the transition probability matrix (P (i, j))N ×N ; here P (i, j) = P (r (tk+1 ) = j|r (tk ) = i) is the ijth entry of the matrix e(tk+1 −tk )Q . Thus we could use the following recursion procedure to simulate the discrete Markov chain {r (tk )}k≥1 ; suppose r (tk ) = i1 and generate a random number ξ which is uniformly distributed in [0, 1], then we define
r (tk+1 ) =
i 2 ,
if i2 ∈ S − {N } and
N ,
if
i2 −1
P (i1 , j) ≤ ξ <
j =1 N −1
i2
P (i1 , j),
j =1
P (i1 , j) ≤ ξ .
j =1
Repeating this procedure a trajectory of {r (tk )}k≥1 can be simulated. Now we can introduce a new family of strong predictor–corrector Euler–Maruyama (PCEM) methods for SDEwMSs. Given initial values Yt0 = y0 ∈ R and rt0 = i0 ∈ S, the proposed family of strong PCEM is given by the predictor
Yt
k+1
= Ytk + f (Ytk , rtk )∆ + g (Ytk , rtk )∆Wtk ,
(2.1)
and by the corrector Ytk+1 = Ytk + {θ f¯η ( Ytk+1 , rtk ) + (1 − θ )f¯η (Ytk , rtk )}∆ + {ηg ( Ytk+1 , rtk ) + (1 − η)g (Ytk , rtk )}∆Wtk .
(2.2)
Here parameters θ , η ∈ [0, 1] denote the degree of implicitness in the drift and the diffusion coefficients, respectively, and f¯η (x, i) is defined as f¯η (x, i) = f (x, i) − ηg (x, i)
∂ g (x, i) , ∂x
η ∈ [0, 1],
(2.3)
which is called the corrected drift function. This scheme can be written in the form Ytk+1 = Ytk + f (Ytk , rtk )∆ + g (Ytk , rtk )∆Wtk +
4
Rl,k ,
(2.4)
l=1
where R1,k = θ{f ( Ytk+1 , rtk ) − f (Ytk , rtk )}∆,
R2,k = −θ η
(2.5)
g ( Ytk+1 , rtk )
∂ g ( Ytk+1 , rtk ) ∂ g (Ytk , rtk ) − g (Ytk , rtk ) ∆, ∂x ∂x
(2.6)
∂ g (Ytk , rtk ) ∆, ∂x = η g ( Ytk+1 , rtk ) − g (Ytk , rtk ) ∆Wtk .
R3,k = −ηg (Ytk , rtk )
(2.7)
R4,k
(2.8)
For each t ∈ [tk , tk+1 ), let Y¯ (t ) = Ytk ,
r¯ (t ) = rtk ,
R¯ l (t ) =
R l ,k
∆
,
l = 1, 2, 3,
R¯ 4 (t ) =
R4,k
∆Wtk
.
(2.9)
Therefore, we can define the continuous approximation solution Y (t ) on the entire interval [0, T ] by Y (t ) = y0 +
t
f (Y¯ (s), r¯ (s))ds + 0
t
g (Y¯ (s), r¯ (s))dW (s) + 0
3 l =1
t
R¯ l (s)ds + 0
t
R¯ 4 (s)dW (s). 0
(2.10)
8
H. Li et al. / Journal of Computational and Applied Mathematics 237 (2013) 5–17
Note that Y (tk ) = Y¯ (tk ) = Ytk , which means that Y (t ) and Y¯ (t ) coincide with the discrete approximate solution at the gridpoints. Remarks 2.2. The major advantage of the above PCEM approximate schemes is that there are flexible degrees of implicitness parameters θ and η to choose for simulating paths properly. For the case θ = η = 0 one recovers the Euler–Maruyama scheme which is well discussed in [10]. 3. Convergence with the global Lipschitz (GL) and linear growth (LG) conditions In this section, we will prove that the numerical solution Y converges to the exact solution y in L2 as step size ∆ ↓ 0, and the order of convergence is one-half. To begin with, we need the following GL condition and LG condition for f¯η (·, ·).
(H 3)GL: There exists a constant L3 > 0, for all (x, i), (y, i) ∈ R × S, such that |f¯η (x, i) − f¯η (y, i)|2 ≤ L3 |x − y|2 . (H 4)LG: There exists a constant L4 > 0, for all (x, i) ∈ R × S, such that |f¯η (x, i)|2 ≤ L4 (1 + |x|2 ). Remarks 3.1. For the notation convenience, we denote L = max{L1 , L2 , L3 , L4 }. We are now ready to present the key results of this section which are stated as following. Theorem 3.1. Assume the SDEwMS (1.1) defined on (Ω , F , {Ft }t ≥0 , P) satisfying (a) W , r are independent, (b) f (·, ·), g (·, ·), f¯η (·, ·) satisfy conditions (H 1), (H 2), (H 3) and (H 4), with initial conditions y(0) = y0 ∈ R and r (0) = i0 ∈ S; then the unique strong solution y and the numerical solution Y obtained in Section 2.2 satisfying:
E
sup |Y (t ) − y(t )|2
≤ Ce40L(4+T )T ∆ + o(∆),
(3.1)
0≤t ≤T
where C = 240ML2 (4 + T )T + 10ML(4L + T + 4LT )T + max(−qii ) · 40ML(4 + T ), i∈S
is a positive constant, dependent on y0 , Q , L and T , but independent of ∆ with M = (1 + 7|y0 |2 )e14L(2+8L+T +4LT )T .
(3.2)
In order to give the proof of this theorem, we first provide a number of useful lemmas. The first two lemmas show that the continuous approximation has bounded moments in a strong sense. The latter lemmas play an important role in proving the strong convergence result, which mainly refer to Bruti-Liberati and Platen [13]. Lemma 3.1. Under conditions (H 1), (H 2), (H 3) and (H 4), we have
E
sup |Y (t )|
2
≤ M,
(3.3)
0≤t ≤T
where M is given by (3.2). Proof. By Eq. (2.10), the fundamental inequality and Hölder inequality, for any T ′ ∈ [0, T ], we have
E
sup |Y (t )|
2
2
≤ 7|y0 | + 7TE
sup
0≤t ≤T ′
0≤t ≤T ′
+ 7T
3 l=1
E
t
|R¯ l (s)|2 ds + 7E
sup 0≤t ≤T ′
0≤t ≤T
0 t
2 ¯ |f (Y (s), r¯ (s))| ds + 7E sup
0
t 2 g (Y¯ (s), r¯ (s))dW (s) ′ 0
t 2 sup R¯ 4 (s)dW (s) . ′
0 ≤t ≤T
0
(3.4)
H. Li et al. / Journal of Computational and Applied Mathematics 237 (2013) 5–17
9
For the second part of Eq. (3.4), by the LG condition, Eq. (2.9) and Y (tk ) = Y¯ (tk ) = Ytk , we have
E
T′
|f (Y¯ (s), r¯ (s))| ds ≤ L 2
sup 0 ≤t ≤T ′
t
sup |Y¯ (s)|2
1+E
0
du
0≤s≤u
0 T′
≤L
sup |Y (s)|2
1+E
du.
(3.5)
0≤s≤u
0
Similarly, for the third part of Eq. (3.4), by Doob’s martingale inequality,
E
t 2 t 2 sup g (Y¯ (s), r¯ (s))dW (s) ≤ 4E sup |g (Y¯ (s), r¯ (s))| ds ′ ′
0 ≤t ≤T
0 ≤t ≤T
0
0
T′ ≤ 4L
1+E
sup |Y (s)|
2
du.
(3.6)
0≤s≤u
0
Now, we focus our attention on the remaining parts of (3.4). By Eqs. (2.5), (2.9), the GL condition and the LG condition, we have
E
t
sup
0≤t ≤T ′
2 2 ¯ |R1 (s)| ds ≤ θ LE
0
≤ 2θ LE 2
T′
E
0≤tk ≤tnu
0
T′
E
sup
≤ 2θ L (∆ + ∆)E 2
|f (Ytk , rtk )∆| |Ftk + E 2
0≤tk ≤tnu
0
2 2
2 sup Ytk+1 − Ytk du
T′
|g (Ytk , rtk )∆Wtk | |Ftk du
sup 0≤tk ≤tnu
1 + sup |Ytk | Ftnu du . 0≤tk ≤tnu 2
E 0
2
(3.7)
So by the assumption 0 < ∆ < 1, θ ∈ [0, 1] and Y (tk ) = Y¯ (tk ) = Ytk , we have
E
|R¯ 1 (s)| ds ≤ 4L2 2
sup 0≤t ≤T ′
t
T′
sup |Y (s)|2
1+E
0
du.
(3.8)
0≤s≤u
0
With the same kind of the steps above, from Eq. (2.6), we can easily show that
E
|R¯ 2 (s)| ds ≤ 4L2 2
sup 0≤t ≤T ′
t
T′
sup |Y (s)|2
1+E
0
du.
(3.9)
0≤s≤u
0
Second, by the LG condition, η ∈ [0, 1], Eqs. (2.7), (2.9) and Y (tk ) = Y¯ (tk ) = Ytk , we obtain
E
|R¯ 3 (s)| ds ≤ η2 L 2
sup 0≤t ≤T ′
t
0
T′
1+E T′
du
0≤s≤u
0
≤L
sup |Y¯ (s)|2
1+E
sup |Y (s)|2
du.
(3.10)
0≤s≤u
0
Moreover, by Doob’s martingale inequality, Eq. (2.8) and the GL condition, we have
E
t 2 t 2 ¯ ¯ sup R4 (s)dW (s) ≤ 4E sup |R4 (s)| ds ′ ′
0≤t ≤T
0≤t ≤T
0
≤ 4η LE 2
0
T′
E 0
2 sup Ytk+1 − Ytk du .
(3.11)
0≤tk ≤tnu
Thus, with the similar step in Eqs. (3.7), (3.8) and η ∈ [0, 1], we obtain
E
t 2 sup R¯ 4 (s)dW (s) ≤ 16L2 ′
0≤t ≤T
0
T′
1+E
0
sup |Y (s)|2 0≤s≤u
du.
(3.12)
10
H. Li et al. / Journal of Computational and Applied Mathematics 237 (2013) 5–17
Substituting (3.5), (3.6), (3.8), (3.9), (3.10), (3.12) into (3.4), we get
sup |Y (t )|
1+E
2
0≤t ≤T ′
≤ 1 + 7|y0 |2 + 14L(2 + 8L + T + 4LT ) ×
T′
1+E
0
sup |Y (s)|2
du.
(3.13)
0≤s≤u
Therefore, from the Gronwall inequality we obtain
E
sup |Y (t )|2
≤ (1 + 7|y0 |2 )e14L(2+8L+T +4LT )T .
(3.14)
0≤t ≤T
The proof is complete.
Lemma 3.2. Under conditions (H 1), (H 2), (H 3) and (H 4), we have
E
max |Ytk |2
0≤k≤nT
≤ M,
(3.15)
where M is given by (3.2). This is an immediate result of Lemma 3.1, since Y (tk ) = Y¯ (tk ) = Ytk . Lemma 3.3. Under conditions (H 1), (H 2), (H 3) and (H 4), there exists a constant M which is dependent on y0 , L and T , but independent of ∆, such that 4 l=1
2 Rl,k ≤ ML(4L + T + 4LT )T ∆ + o(∆), E max 1≤n≤nT 0≤k≤n−1
(3.16)
where M is given by (3.2). Proof. By the Cauchy–Schwarz inequality and the GL condition, we can obtain from Eq. (2.5) that
2 R1,k E max 1≤n≤nT 0≤k≤n−1 2 = E max θ {f ( Ytk+1 , rtk ) − f (Ytk , rtk )}∆ 1≤n≤nT 0≤k≤n−1 2 2 |θ ∆| ≤ E max f ( Ytk+1 , rtk ) − f (Ytk , rtk ) 1≤n≤nT
0≤k≤n−1
0≤k≤n−1
≤ θ ∆E 2
∆
0≤k≤nT −1
f ( Yt
2 , rtk ) − f (Ytk , rtk ) k+1
0≤k≤nT −1
≤ θ LT ∆E 2
Yt
k+1
2 − Ytk
0≤k≤nT −1
≤ 2θ LT ∆E 2
2 2 E |f (Ytk , rtk )∆| |Ftk + E |g (Ytk , rtk )∆Wtk | |Ftk .
(3.17)
0≤k≤nT −1
Then by using the Cauchy–Schwarz inequality, the Itô’s isometry, the LG condition, Lemma 3.2 and θ ∈ [0, 1], we get
2 tn T E max R1,k ≤ 2θ 2 L2 T (∆ + 1)∆E E 1 + |Ytnz |2 |Ftnz dz 1≤n≤nT 0 0≤k≤n−1 tn T ≤ 2θ 2 L2 T (∆ + 1)∆ E 1 + max |Ytn |2 dz 0
≤ 2ML T ∆ + o(∆). 2
2
0≤n≤nT
(3.18)
H. Li et al. / Journal of Computational and Applied Mathematics 237 (2013) 5–17
11
With the similar steps as in (3.17) and (3.18), by the GL condition and the LG condition, we have
2 E max R2,k ≤ 2ML2 T 2 ∆ + o(∆). 1≤n≤nT 0≤k≤n−1
(3.19)
It is also easy to show by the LG condition and Lemma 3.2 that
2 E max R3,k ≤ MLT 2 ∆ + o(∆). 1≤n≤nT 0≤k≤n−1
(3.20)
By using Doob’s martingale inequality and Eq. (2.8), we have
2 2 tk + 1 E max R4,k = E max g ( Ytk+1 , rtk ) − g (Ytk , rtk ) dW (z ) η 1≤n≤nT 1≤n≤nT t k 0≤k≤n−1 0≤k≤n−1 2 T G(z )dW (z ) , ≤ 4η2 E 0
where
G(z ) =
g ( Ytk+1 , rtk ) − g (Ytk , rtk ) I[tk ,tk+1 ) (z ),
0≤k≤nT −1
and IΛ be the indicator function for set Λ. Obviously, G(z ) ≤
g ( Ytk+1 , rtk ) − g (Ytk , rtk ) I[0,T ] (z ).
max 0≤k≤nT −1
By Itô’s isometry, the GL condition, with similar steps in (3.17), (3.18) and η ∈ [0, 1], we have
2 T E |G(z )|2 dz R4,k ≤ 4η2 E max 1≤n≤nT 0 0≤k≤n−1 T 2 ≤ 4η2 E max g ( Ytk+1 , rtk ) − g (Ytk , rtk ) I[0,T ] (z )dz 0≤k≤nT −1 0 2 ≤ 4η2 LTE max Ytk+1 − Ytk
0≤k≤nT −1
≤ 4ML2 T ∆ + o(∆). Thus, substituting (3.18), (3.19), (3.20), (3.21) into (3.17), the required assertion follows.
(3.21)
Lemma 3.4. Under conditions (H 1), (H 2), (H 3) and (H 4), for any t ∈ [tk , tk+1 ), we have 3 l =1
E sup t ≤s≤t k
s
2 s 2 R4,k dz + E sup dW (z ) ≤ o(∆). ∆ tk ≤s≤t tk ∆Wtk
R l ,k
tk
(3.22)
The proof of this lemma is similar to that in Lemma 3.3. Now we are in a position to prove Theorem 3.1. Proof of Theorem 3.1. From Eq. (2.10), we have Z (t ) = E
sup |Y (s) − y(s)|2
0≤s≤t
s s = E sup (f (Y¯ (z ), r¯ (z )) − f (y(z ), r (z )))dz + (g (Y¯ (z ), r¯ (z )) − g (y(z ), r (z )))dW (z ) 0≤s≤t 0 0 +
3 l =1
s
R¯ l (z )dz + 0
0
s
2 R¯ 4 (z )dW (z ) .
(3.23)
12
H. Li et al. / Journal of Computational and Applied Mathematics 237 (2013) 5–17
Let n = [t /∆], the integer part of t /∆. Then, by the fundamental inequality, the Hölder inequality, Doob’s martingale inequality and Eq. (2.9), we have Z (t ) = E
sup |Y (s) − y(s)|2
0≤s≤t
t |g (Y¯ (z ), r¯ (z )) − g (y(z ), r (z ))|2 dz |f (Y¯ (z ), r¯ (z )) − f (y(z ), r (z ))|2 dz + 40E 0 0 2 s 2 4 3 R l ,k + 10 E max Rl,k + 10 E sup dz 1≤m≤n tk ≤s≤t tk ∆ l =1 0≤k≤m−1 l =1
t
≤ 10TE
+ 10E sup t ≤s≤t k
s tk
R4,k
∆Wtk
2 dW (z ) .
(3.24)
We focus on the last three parts of the right side. From Lemmas 3.3 and 3.4, we have 10
4 l =1
2 s 2 3 R l ,k + 10 E sup R l ,k dz E max 1≤m≤n ∆ tk ≤s≤t t k l =1 0≤k≤m−1
+ 10E sup t ≤s≤t k
s tk
2 dW (z ) ≤ 10ML(4L + T + 4LT )T ∆ + o(∆). ∆Wtk R4,k
(3.25)
Let t ∈ [tk , tk+1 ). From (2.1)–(2.10), it is easy to know that Y¯ (tk ) and I{r (t )̸=r (tk )} are conditionally independent with respect to the σ -algebra generated by r (tk ). So by the same procedures as in Theorem 3.1 in [10] or Theorem 4.1 in [4], we can obtain t
|f (Y¯ (s), r¯ (s)) − f (y(s), r (s))|2 ds ≤ 2L
E 0
E |Y¯ (s) − y(s)|2 ds + max(−qii ) · 4ML∆ + o(∆),
(3.26)
E |Y¯ (s) − y(s)|2 ds + max(−qii ) · 4ML∆ + o(∆).
(3.27)
i∈S
0 t
|b(Y¯ (s), r¯ (s)) − b(y(s), r (s))|2 ds ≤ 2L
E
t
0
t
i∈S
0
Substituting (3.25), (3.26), (3.27) into (3.24) shows that
E
sup |Y (s) − y(s)|2
≤ 20L(4 + T )
0≤s≤t
t
E |Y¯ (s) − y(s)|2 ds 0
+ 10ML(4L + T + 4LT )T + max(−qii ) · 40ML(4 + T ) ∆ + o(∆).
(3.28)
i∈S
Note that E |Y¯ (s) − y(s)|2 ≤ 2E |Y (s) − y(s)|2 + 2E |Y¯ (s) − Y (s)|2 .
(3.29)
Suppose tk ≤ s < tk+1 , by (2.4), (2.10), Lemmas 3.3 and 3.4, we can then show in the same way as above that E |Y (s) − Y¯ (s)|2 ≤ E |Ytk+1 − Ytk |2 ≤ 6ML∆ + o(∆).
(3.30)
Substituting (3.29), (3.30), into (3.28) immediately shows that
E
sup |Y (s) − y(s)|2
0≤s≤t
t ≤ 40L(4 + T ) (E |Y (s) − y(s)|2 + E |Y¯ (s) − Y (s)|2 )ds 0 + 10ML(4L + T + 4LT )T + max(−qii ) · 40ML(4 + T ) ∆ + o(∆) i∈S t 2 ≤ 40L(4 + T ) E sup |Y (s) − y(s)| du 0
0≤s≤u
+ 240ML2 (4 + T )T + 10ML(4L + T + 4LT )T + max(−qii ) · 40ML(4 + T ) ∆ + o(∆). i∈S
(3.31)
H. Li et al. / Journal of Computational and Applied Mathematics 237 (2013) 5–17
13
Therefore, from the Gronwall inequality we obtain that
E
sup |Y (t ) − y(t )|
2
≤ Ce40L(4+T )T ∆ + o(∆),
0 ≤t ≤T
where C = 240ML2 (4 + T )T + 10ML(4L + T + 4LT )T + max(−qii ) · 40ML(4 + T ). i∈S
The proof is complete.
4. The general multi-dimensional case The results derived in Section 3 can be easily generalized to the multi-dimensional case, we just summarize the related numerical schemes and the convergence results in this section. Consider the solution y = {(y1 (t ), . . . , yd (t ))T , t ≥ 0} of the d-dimensional SDEwMS y(t ) = y(0) +
t
f (y(s), r (s))ds +
0
m j=1
t
g j (y(s), r (s))dW j (s), 0
for t ≥ 0. Here y(0) ∈ Rd and r (0) = r0 ∈ S denotes the deterministic initial values, W j = {W j (t ), t ≥ 0}, j ∈ {1, 2, . . . , m} is a Wiener process. r = {r (t ), t ≥ 0} is a Markov chain. The function f : Rd × S → Rd has the kth component f k (·, ·). The function g j : Rd × S → Rd , j ∈ {1, 2, . . . , m} has the kth component g k,j (·, ·). We define the function f¯η for η = {η1 ∈ [0, 1], . . . , ηd ∈ [0, 1]} with the kth component f¯ηk (x, i) = f k (x, i) − ηk
m d
g i ,j 1 ( x , i )
j1 ,j2 =1 i=1
∂ g i,j2 (x, i) ∂ xi
for (x, i) ∈ R × S, which satisfies the following GL condition and the LG condition. (H 3′ )GL. There exists a constant L5 > 0, for all (x, i), (y, i) ∈ Rd × S, such that |f¯η (x, i) − f¯η (y, i)|2 ≤ L5 |x − y|2 . d
(H 4′ )LG. There exists a constant L6 > 0, for all (x, i) ∈ Rd × S, such that |f¯η (x, i)|2 ≤ L6 (1 + |x|2 ). The kth component of the proposed family of strong PCEM schemes is given by the predictor
Ytk
n+1
= Ytkn + f k (Ytn , rtn )∆ +
m
g k,j (Ytn , rtn )∆Wtn , j
j =1
and by the corrector Ytkn+1 = Ytkn + {θk f¯ηk ( Ytn+1 , rtn ) + (1 − θk )f¯ηk (Ytn , rtn )}∆ +
m {ηk g k,j ( Ytn+1 , rtn ) + (1 − ηk )g k,j (Ytn , rtn )}∆Wtn , j =1
for θk , ηk ∈ [0, 1], k ∈ {1, 2, . . . , d}. Hence we can define the approximation solution Y similarly to Eq. (2.10); then we can derive the following theorem analogically. Theorem 4.1. Assume the SDEwMS (1.1) defined on (Ω , F , {Ft }t ≥0 , P ) satisfying (a) W , r are independent, (b) f (·, ·), g (·, ·), f¯η (·, ·) satisfy conditions (H 1), (H 2), (H 3′ ) and (H 4′ ), with initial conditions y(0) = y0 ∈ Rd and r (0) = i0 ∈ S; then the unique strong solution y and the numerical solution Y satisfy:
E
sup |Y (t ) − y(t )|
0 ≤t ≤T
2
≤ C ∆ + o(∆),
where C is a positive constant dependent on y0 , Q , L and T , but independent of ∆. 5. Numerical stability In this section we consider numerical stability issues, extending the analysis by Platen and Shi [24] to the Markovian switching case. When simulating discrete time approximations of solutions of SDEwMSs, numerical stability is clearly as
14
H. Li et al. / Journal of Computational and Applied Mathematics 237 (2013) 5–17
important as numerical efficiency. There have been various efforts made in the literature trying to study numerical stability for a given scheme approximating solutions of SDEs; see, for instance, Bruti-Liberati and Platen [13], Higham [25], Hofmann and Platen [26], Saito and Mitsui [27], Platen and Shi [24], and Schurz [17–19] among others. Generally, for analyzing numerical stability, some specifically designed test equations are necessary to be introduced to allow insight into the stability behavior of the numerical methods. In deterministic theory the test equation is generally scalar and linear, and indicates that the stability theory has relevance to the behavior near equilibrium on a nonlinear problem. In stochastic theory, the situation is more complex, it is not very clear to which extend the study of a scalar linear test problem is relevant to the nonlinear equations. In this section, we only specialize to a linear test problem and show the connections of stability properties between SDEs and SDEwMSs. For the purpose of comparison, we first consider that the test SDE which has been used by Platen and Shi [24] is a linear and scalar Itô-SDE with multiplicative noise defined as dX (t ) =
1−
3 2
α λX (t )dt + α|λ|X (t )dW (t ),
t ≥ 0,
(5.1)
where X (0) = x0 ≥ 0, λ < 0 and α ∈ [0, 1). Eq. (5.1) can be derived by linearizing about the fixed point X (t ) ≡ 1 of the following SDE
3
dX (t ) = − 1 −
2
α λX (t )(1 − X (t ))dt −
α|λ|X (t )(1 − X (t ))dW (t ),
t ≥ 0,
(5.2)
which is a normalized version of a stochastic population model (see Higham [25]). It is worth noting that the test Eq. (5.1) is just a special case of more general bilinear systems discussed in [17] when investigating the mean square stability and A-stability of numerical methods for SDEs. Obviously, Eq. (5.1) has an explicit solution of the form
X (t ) = x0 exp (1 − α)λt +
α|λ|W (t ) ,
t ≥ 0.
(5.3)
As a unified criterion, Platen and Shi [24] proposed the concept of p-stability criterion, which means that a process is p-stable if in the long run its pth moment vanishes. Hence, for p > 0, λ < 0, p-stability in Eq. (5.1) may be characterized by lim E (|X (t )|p ) = 0
t →∞
0≤α<
iff
1 1+
p 2
.
Since for different combinations of values of λ∆, α and p with given time step size ∆, a discrete time approximation Y and the original continuous process X have different stability properties. To explore these differences, the concept of stability region is introduced. The stability region, denoted by Γ , is determined by those triplets (λ∆, α, p) ∈ (−∞, 0) × [0, 1) × (0, ∞) for which the discrete time approximation Y is p-stable with time step size ∆, when applied to the test Eq. (5.1). By defining the random variable
Y (tn+1 ) , Gn+1 (λ∆, α) = Y (tn )
which is called the transfer function of the approximation Y at time tn , Platen and Shi [24] derived the following useful result which can determine the stability regions for given schemes by the following theorem. Theorem 5.1. A discrete time approximation is for given λ < 0, α ∈ [0, 1) and p > 0, p-stable if and only if E ((Gn+1 (λ∆, α))p ) < 1. In the spirit of Platen and Shi [24], the stability regions for a range of schemes of SDEwMSs are discussed in this paper. Here we consider that the test process Xr = {X (t , r (t )), t ≥ 0} satisfies the linear and scalar Itô-SDEwMS with multiplicative noise dX (t ) =
1−
3 2
α(r (t )) λ(r (t ))X (t )dt + α(r (t ))|λ(r (t ))|X (t )dW (t ),
(5.4)
for t ≥ 0, where r = {r (t ), t ≥ 0} is a Markov chain taking values in a finite state space S = {1, 2, . . . , N }, r (0) = i0 ∈ S , X (0) = x0 ≥ 0 and (α(r (t )), λ(r (t ))) ∈ {(α(i), λ(i)), α(i) ∈ [0, 1), λ(i) < 0, i = 1, 2, . . . , N }. It is well known that the explicit solution of the test equation (5.4) is (see Mao and Yuan [4]) X (t ) = x0 exp
t
(1 − α(r (t )))λ(r (t ))dt + 0
t α(r (t ))|λ(r (t ))|dW (t ) , 0
for t ≥ 0. Now, we introduce the following stability criterion.
(5.5)
H. Li et al. / Journal of Computational and Applied Mathematics 237 (2013) 5–17
15
Definition 5.1. Given p > 0, for all y0 ∈ R and r0 ∈ S, a process Yr = {Y (t , r (t )), t ≥ 0} is called state-p-stable if for each i = 1, 2, . . . , N , Yi = {Y (t , i), t ≥ 0} satisfies lim E (|Y (t , i)|p ) = 0.
t →∞
Remarks 5.1. Actually, if the process Yr is state-p-stable, then it must be p-stable, that is lim E (|Y (t , r (t ))|p ) = 0,
t →∞
∀y0 ∈ R, r0 ∈ S ,
but the reverse is not true; see Pang et al. [28]. For the convenience of numerical comparison, we use this state-p-stability concept in this paper. For more details about p-stability of SDEwMSs, we refer the reader to Pang et al. [28] and references therein. Definition 5.2. The state-stability region Γs is determined by those triplets (λ∆, α, p) ∈ (−∞, 0) × [0, 1) × (0, ∞) for which the discrete time approximation Yr is state-p-stable, when applied to the test Eq. (5.4), for each i = 1, 2, . . . , N, (α(i), λ(i)∆, p) ∈ Γs with time step size ∆. Then we can obtain the following conclusion. Theorem 5.2. The state-stability region Γs , which is generated by one algorithm applied to the test equation (5.4), is the same as the stability region Γ , which is generated by the same algorithm applied to the test equation (5.1). Proof. Γ ⊆ Γs is obvious. To prove Γs ⊆ Γ , suppose for all i = 1, 2, . . . , N, (α(i), λ(i)∆, p) ∈ Γs . By Definition 5.2, we can see that Yr is state-p-stable. And by Definition 5.1, we know for each i = 1, 2, . . . , N , Yi = {Y (t , i), t > 0} is p-stable. So (α(i), λ(i)∆, p) ∈ Γ , for all i ∈ {1, 2, . . . , N }. Immediately, we have Γs ⊆ Γ . Remarks 5.2. By the conclusions derived by Platen and Shi [24], we can also see that the PCEM methods are more efficient than the EM method under these conditions in SDEwMSs case when the parameters selection are reasonable. 6. Numerical examples In this section, we discuss two numerical examples to illustrate our theory established in the previous sections. Let us now consider several combinations of parameters θ and η in Eq. (2.2), the names of the methods listed below are similar to those used by Bruti-Liberati and Platen [13] and Platen and Shi [24]. (1) θ = 0, η = 0 (called the EM scheme), (2) θ = 12 , η = 21 (called the symmetric PCEM scheme), (3) (4) (5) (6)
θ θ θ θ
= 12 , η = 0 (called the semi-drift-implicit PCEM scheme), = 1, η = 0 (called the drift-implicit PCEM method), = 0, η = 21 (called the semi-diffusion-implicit PCEM scheme), = 1, η = 1 (called the fully implicit PCEM scheme).
For a given problem, we will compare the simulated paths for these different degrees of implicitness. If these paths differ significantly from each other, then some numerical stability problem is likely to be present and one needs to make an effort in providing extra numerical stability for further research. Example 6.1. Let W be a Wiener process. Let r be a right-continuous Markov chain taking values in S = {1, 2} with generator
Q =
−0.5 0.5
0.5 . −0.5
And W and r are assumed to be independent. Consider a one-dimensional linear SDEwMS dy(t ) = y(t )a(r (t ))ds + y(t )b(r (t ))dW (t ),
t ≥ 0,
(6.1)
where a(i), b(i), for all i ∈ S, are constants, and the initial value y0 ∈ R, r (0) = i0 ∈ S. It is well known that Eq. (6.1) has an explicit solution y(t ) = y0 exp
t
a(r (s))ds + 0
t
b(r (s))dW (s) − 0
1 2
t
b (r (s))ds , 2
t ≥ 0.
(6.2)
0
For simulation reason, it is convenient to transform (6.2) into the following recursion form with y(t0 ) = y0 ,
√
y(tk+1 ) ≈ y(tk ) exp (tk+1 − tk )a(rtk ) + b(rtk ) tk+1 − tk ξk+1 −
1 2
(tk+1 − tk )b2 (rtk ) .
(6.3)
16
H. Li et al. / Journal of Computational and Applied Mathematics 237 (2013) 5–17
Table 1 Estimation of the errors between the numerical and exact solutions.
E suptk ∈[0,10] |Y(θi ,ηi ) (tk ) − y(tk )|2 ∆ \ θi , ηi
0, 0
1 2
,
0.1 0.05 0.01 0.005 0.001
4.0e−03 1.0e−03 4.7e−05 1.4e−05 1.3e−06
5.8e−06 5.0e−07 4.2e−09 7.8e−10 2.4e−11
1 2
1 2
,0
2.1e−04 7.5e−05 1.0e−05 4.9e−06 9.5e−07
1, 0
0,
4.7e−03 1.1e−03 4.7e−05 1.4e−05 1.3e−06
3.6e−03 8.8e−04 3.4e−05 8.6e−06 3.4e−07
1, 0
0,
5.4e+07 1.3e+07 5.2e+05 1.3e+05 5.3e+03
2.9e+07 9.6e+06 4.9e+05 1.3e+05 5.2e+03
1 2
1, 1 5.2e−03 1.2e−03 5.0e−05 1.5e−05 1.3e−06
Table 2 Estimation of the errors between the numerical and exact solutions.
E suptk ∈[0,10] |Y(θi ,ηi ) (tk ) − y(tk )|2 ∆ \ θi , ηi
0, 0
1 2
,
0.1 0.05 0.01 0.005 0.001
3.0e+07 9.9e+06 5.1e+05 1.3e+05 5.5e+03
6.1e+04 5.6e+03 3.8e+01 6.3e+00 1.7e−01
1 2
1 2
,0
1.5e+05 3.0e+04 2.3e+03 9.5e+02 1.6e+02
1 2
1, 1 5.4e+07 1.3e+07 5.2e+05 1.3e+05 5.3e+03
Here each ξk+1 is an independent standard √ normal random variable, so that the increment ∆Wtk = W (tk+1 ) − W (tk ) of the Wiener process W can be generated by tk+1 − tk ξk+1 . Notice that y(tk+1 ) in (6.3) is not the exact value of y(t ) at the division points tk+1 , because r (s) is not necessarily constant on s ∈ [tk , tk+1 ]. However, since P {r (tk+1 ) = i|r (tk ) = i} = 1 + qii (tk+1 − tk ) + o(tk+1 − tk ) → 1 as ∆ ↓ 0, for sufficiently small ∆, it is reasonable to use (6.3) as an approximation of the exact solution of y. Choose initial values y0 = 10, i0 = 1, T is fixed at 10, a(1) = −0.5, a(2) = −1,
b(1) = 0.1 b(2) = 0.1.
By applying the previously described procedure in Section 2, the trajectory of the approximate solution Y with given step size ∆ can be constructed. To carry out the numerical simulation we successively choose the step size ∆ as in the following Table 1, and for each ∆, we repeatedly simulate and compute suptk ∈[0,10] |Y(θi ,ηi ) (tk ) − y(tk )|2 , (i = 1, 2, 3, 4, 5, 6) for 106 times, then calculate the
sample mean E [suptk ∈[0,10] |Y(θi ,ηi ) (tk ) − y(tk )|2 ] (i = 1, 2, 3, 4, 5, 6). The results are listed in the following Table 1. Clearly, we can easily show that Example 6.1 satisfies the conditions of Eq. (5.4), so it will be state-p-stable, and the simulation results listed in Table 1 just illustrate the stability properties in a certain extent. On the one hand, the numerical method reveals that the numerical solution Y defined by the strong PCEM methods converge to the exact solution y in L2 as step size ∆ ↓ 0. On the other hand, when the parameters’ selection is reasonable, the PCEM methods are much more efficient than the EM method, which strongly demonstrate our results. Example 6.2. Consider the same 1-dimensional linear SDEwMS dy(t ) = y(t )a(r (t ))ds + y(t )b(r (t ))dW (t ),
t ≥ 0.
(6.4)
This time, choose initial values y0 = 10, r0 = 1, T = 10, and
a(1) = 0.5, a(2) = 1,
b(1) = 0.1 b(2) = 0.1
By applying the previously described procedure, the trajectory of the approximate solution Y with given step size ∆ can be constructed. To carry out the numerical simulation we successively choose the step size ∆ as in the following Table 2, and for each ∆, we repeatedly simulate and compute suptk ∈[0,10] |Y(θi ,ηi ) (tk ) − y(tk )|2 , (i = 1, 2, 3, 4, 5, 6) for 106 times; then calculate the
sample mean E [suptk ∈[0,10] |Y(θi ,ηi ) (tk ) − y(tk )|2 ] (i = 1, 2, 3, 4, 5, 6). The results are listed in the following Table 2. Since we can easily show that Example 6.2 does not satisfy the conditions of Eq. (5.4), it will not be state-p-stable, even not p-stable, the computer simulation results illustrate this point to some extent. Nevertheless, when the parameters θ and η are selected reasonably, for example, θ = η = 12 in this example, some PCEM methods are much more efficient than the EM method to a certain extent.
H. Li et al. / Journal of Computational and Applied Mathematics 237 (2013) 5–17
17
7. Concluding remarks In this paper, we investigated the PCEM methods based on jump-adapted time discretizations to handle scenario simulation of the exact solutions of SDEwMSs under global Lipschitz conditions. We showed that these numerical schemes give strong L2 -convergence with order 0.5. Furthermore, we established the numerical state p-stability criterion of the PCEM methods for a specifically designed linear test equation. It has been shown that such schemes exhibit good computational efficiency when the parameters are selected. There are several interesting problems that deserve further investigation. For example, it would clearly be of interest to extend the strong convergence theory to the case where coefficients are not globally Lipschitz continuous, it remained an open question whether the approximation using the PCEM methods also converges in the strong sense if the coefficients of the equations are not globally Lipschitz continuous. Another important problem is how to develop and analyze the various numerical stability criteria for such schemes when making decisions about the selection of an appropriate combination of parameters for strong PCEM methods for a given simulation task. Since for each possible combination of these parameters, the schemes are of strong order 0.5. Nevertheless, the numerical results show that the symmetric PCEM scheme exhibits good numerical properties. The most common method for clarifying this issue is likely to study the shapes of stability regions for related approximating schemes. It is worthwhile to compare the observed stability regions with these numerical methods, which are known to show good numerical stability properties. We hope to address at least some of the specific issues in subsequent work. Acknowledgments The author gratefully acknowledges many useful comments and suggestions from the referees and editors. References [1] C.T. Cheng, C.P. Ou, K.W. Chau, Combining a fuzzy optimal model with a genetic algorithm to solve multi-objective rainfall-runoff model calibration, J. Hydrol. 268 (1–4) (2002) 72–86. [2] M.K. Ghosh, A. Arapostathis, S.I. Marcus, Ergodic control of switching diffusions, SIAM J. Control Optim. 35 (6) (1997) 1952–1988. [3] A. Jobert, L.C.G. Rogers, Option pricing with Markov-modulated dynamics, SIAM J. Control Optim. 44 (6) (2006) 2063–2078. [4] X. Mao, C. Yuan, Stochastic Differential Equations with Markovian Switching, Imperial College Press, 2006. [5] T. Rolski, H. Schmidili, V. Schmidt, J. Teugels, Stochastic Processes for Insurance and Finance, Wiley Online Library, 2000. [6] D.R. Smith, Markov-switching and stochastic volatility diffusion models of short-term interest rates, J. Bus. Econom. Statist. 20 (2) (2002) 183–197. [7] C.L. Wu, K.W. Chau, J.S. Huang, Modelling coupled water and heat transport in a soil-mulch-plant-atmosphere continuum (SMPAC) system, Appl. Math. Model. 31 (2) (2007) 152–169. [8] H. Yang, G. Yin, Ruin probability for a model under Markovian switching regime, in: Probability, Finance and Insurance: Proceedings of a Workshop at the University of Hong Kong, Hong Kong, 15–17 July 2002, World Scientific, 2004, pp. 206–217. [9] M.Y. Zhao, C.T. Cheng, K.W. Chau, G. Li, Multiple criteria data envelopment analysis for full ranking units associated to environment impact assessment, Int. J. Environ. Pollut. 28 (3) (2006) 448–464. [10] C. Yuan, X. Mao, Convergence of the Euler–Maruyama method for stochastic differential equations with Markovian switching, Math. Comput. Simulation 64 (2) (2004) 223–235. [11] G. Yin, Q.S. Song, Z. Zhang, Numerical solutions for jump-diffusions with regime switching, Stochastics 77 (1) (2005) 61–79. [12] N. Bruti-Liberati, E. Platen, Strong approximations of stochastic differential equations with jumps, J. Comput. Appl. Math. 205 (2) (2007) 982–1001. [13] N. Bruti-Liberati, E. Platen, Srong predictor–corrector Euler methods for stochastic differetial equations, Stoch. Dyn. 8 (3) (2008) 561–581. [14] Z. Hou, J. Tong, Z. Zhang, Convergence of jump-diffusion non-linear differential equation with phase semi-Markovian switching, Appl. Math. Model. 33 (9) (2009) 3650–3660. [15] R. Li, Y. Huo, Convergence and stability of numerical solutions to SDDEs with Markovian switching, Appl. Math. Comput. 175 (2) (2006) 1080–1091. [16] A. Rathinasamy, K. Balachandran, Mean square stability of semi-implicit Euler method for linear stochastic differential equations with multiple delays and Markovian switching, Appl. Math. Comput. 206 (2) (2008) 968–979. [17] H. Schurz, Asymptotical mean square stability of an equilibrium point of some linear numerical solutions with multiplicative noise, Stoch. Anal. Appl. 14 (3) (1996) 313–353. [18] H. Schurz, Stability, Stationarity, and Boundedness of Some Implicit Numerical Methods for Stochastic Differential Equations and Applications, LogosVerlag, Berlin, 1997. [19] H. Schurz, The invariance of asymptotic laws of stochastic systems under discretization, Z. Angew. Math. Mech. 79 (6) (1999) 375–382. [20] J.R. Klauder, W.P. Petersen, Numerical integration of multiplicative-noise stochastic differential equations, SIAM J. Numer. Anal. 22 (1985) 1153–1166. [21] W.P. Petersen, Some experiments on numerical simulations of stochastic differential equations and a new algorithm, J. Comput. Phys. 113 (1) (1994) 75–81. [22] G.N. Milstein, E. Platen, H. Schurz, Balanced implicit methods for stiff stochastic systems, SIAM J. Numer. Anal. 35 (3) (1998) 1010–1019. [23] P.E. Kloeden, E. Platen, Numerical Solution of Stochastic Differential Equations, Springer-Verlag, New York, 1999. [24] E. Platen, L. Shi, On the numerical stability of simulation methods for SDEs, Research Paper Series, Quantitative Finance Research Centre, University of Technology, Sydney, 2008, 234, pp. 1–25. [25] D.J. Higham, Mean-square and asymptotic stability of numerical methods for stochastic ordinary differential equations, SIAM J. Numer. Anal. 38 (3) (2000) 753–769. [26] N. Hofmann, E. Platen, Stability of weak numerical schemes for stochastic differential equations, Comput. Math. Appl. 28 (10–12) (1994) 45–57. [27] Y. Saito, T. Mitsui, Stability analysis of numerical schemes for stochastic differential equations, SIAM J. Numer. Anal. 33 (6) (1996) 2254–2267. [28] S. Pang, F. Deng, X. Mao, Almost sure and moment exponential stability of Euler–Maruyama discretizations for hybrid stochastic differential equations, J. Comput. Appl. Math. 213 (1) (2008) 127–141.