A High Order Accurate and Effective Scheme for Solving Markovian Switching Stochastic Models

: In this paper, we propose a new weak order 2.0 numerical scheme for solving stochastic differential equations with Markovian switching (SDEwMS). Using the Malliavin stochastic analysis, we theoretically prove that the new scheme has local weak order 3.0 convergence rate. Combining the special property of Markov chain, we study the effects from the changes of state space on the convergence rate of the new scheme. Two numerical experiments are given to verify the theoretical results.


Introduction
In recent years, stochastic differential equations with Markovian switching (SDEwMS) has attracted the attention of many scholars. Comparing with stochastic differential equations with jumps (SDEwJ), SDEwMS are not only applied in finance, but they also have many applications in other fields, which include control system, biomathematics, chemistry, and mechanics. Zhou and Mamon [1] proposed an accessible implementation of interest rate models with Markov-switching. Siu [2] proposed a high-order Markov-switching model for risk measurement. He, Qi, and Kao [3] studied the HMM-based adaptive attack-resilient control for Markov jump system and application to an aircraft model. In this paper, we consider the following q-dimensional stochastic differential equations with Markovian switching: dX s = a(s, r s , X s )ds + b(s, r s , X s )dW s , 0 ≤ s ≤ T, where r s ∈ S, the Brownian motion W s = (W 1 s , W 2 s , . . . , W m s ) T s≥0 , is independent of the Markov chain r t (t ≥ 0), a : [0, T] × S × R q −→ R q and b : [0, T] × S × R q −→ R q×m . Qualitative theory of the existence and uniqueness of the solution for SDEwMS has been studied for the past years (see [4,5]). Many scholars have studied the stability of SDEwMS, for example, stability of linear or semi-linear type of Equation (1) has been studied by Basak [6] and Ji [7]. Mao [4] discussed the exponential stability for general nonlinear stochastic differential equations with Markovian switching. Yang, Yin and Li [8] focused on stability analysis of numerical solutions to jump diffusion and jump diffusion with Markovian switching. Ma and Jia [9] considered the stability analysis of linear Itô stochastic differential equations with countably infinite Markovian switchings.
Generally, most of SDEwMS do not have explicit solutions and hence require numerical solutions. Yuan and Mao [10] firstly developed a Euler-Maruyama scheme for solving SDEwMS and estimated the errors between the numerical and exact solutions under Lipschitz conditions. Yuan and Mao [11] proved that the strong rate of convergence of the Euler-Maruyama scheme is equal to 0.5 under non-Lipschitz conditions. Then, many scholars extended to the semi-implicit Euler scheme [12][13][14]. Furthermore, Li [15] and Chen [16] developed an Euler scheme for solving stochastic differential equations with • We propose a new weak order 2.0 scheme and approximate multiple stochastic integral by utilizing the compound Poisson process. • By integration-by-parts formula of Malliavin calculus theory [25], we rigorously prove that the new scheme has local weak order 3.0 convergence rate. We also prove that the convergence rate is related to the maximum state difference and upper bounds of state values. • We give two numerical experiments to confirm our theoretical convergence results and the convergence rate effects from Markov chain space. • In the experiments, we make comparisons with other schemes such as Euler scheme and Runge-Kutta scheme and verify the new scheme is effective and accurate.
Some notations to be used later are listed as follows: • | · | is the norm for vector or matrix defined by |A| 2 = trace(A T A). • C k p (R q , R) is the set of k times continuously differentiable functions which, together with their partial derivatives of order up to k, have at most polynomial growth.
C is a generic constant depending only on the upper bounds of derivatives of a, b and g, and C can be different from line to line.
The paper is organized as follows. In Section 2, we introduce some basic concepts. In Section 3, we propose a new scheme and obtain its local order 3.0 convergence results. In Section 4, two numerical examples are given to verify the theoretical results.

Markov Chain
In this paper, let (Ω, F , {F t } t≥0 , P) be a complete probability space with a filtration {F t } t≥0 satisfying the usual conditions. Let x T denote the transpose of a vector or matrix x. Let {r t , t ≥ 0} be a right-continuous Markov chain on the probability space taking values in a finite state space S = {1, 2, ..., M} with generator Q = (q ij ) M×M given by where δ > 0, q ij ≥ 0 and, for i = j, q ii = − ∑ i =j q ij . For the above Markov chain, N(de, dt) is a Poisson random measure with intensity λ(de)dt on E × [0, T], in which λ is the Lebesgue measure on E , where E = R \ {0} is equipped with its Borel field E. dr t = E l(r t− , e)N(de, dt) with l(i, e) = j − i f or e ∈ ∆ ij and l(i, e) = 0 f or e / ∈ ∆ ij , where ∆ ij are the intervals having length q ij satisfying and so on.
In [17], the discrete Markov chain {r nh , n = 0, 1, ...} given a step size h > 0 is simulated as follows: compute the one-step transition probability matrix Let r 0 = i 0 and generate a random number ξ 1 which is uniformly distributed in [0, 1]. Define where we set 0 ∑ j=1 P i 0 ,j (h) = 0 as usual. Generate independently a new random number ξ 2 , which is again uniformly distributed in [0, 1], and then define
Assumption 1. Assume that there exist two positive constants L and K such that

Malliavin Calculus
Suppose that H is a real separable Hilbert space with scalar product denoted by ·, · H . The norm of an element h ∈ H is denoted by h H . Let B = {W(h), h ∈ H} denote an isonormal Gaussian process associated with the Hilbert space H on (Ω, F , {F t } t≥0 , P).
We call a row vector α = (j 1 , j 2 , ...j l ) with j i ∈ {−1, 0, 1, ..., m} for i ∈ {1, 2, ..., l} a multi-index of length l : l(α) ∈ {1, 2, ..., m} and denote by v the multi-index of length zero (l(v) := 0). Let M be the set of all multi-indices, i.e., A random variable F is Malliavin differentiable if and only if F ∈ D 1,2 , where the space D 1,2 ⊂ L 2 (P) is defined by completion with respect to the norm · 1,2 . The details about the Malliavin derivative of Poisson process can be found in [25]. For any integer p ≥ 1, D k,p is the domain of D k (k ∈ N) in L p (Ω), i.e., D k,p is the closure of the class of smooth random variables F with respect to the norm Lemma 1 (Product rule). Let F, G ∈ D 1,2 , then FG ∈ D 1,2 , and we have Lemma 2 (Chain rule and Duality formula). Let F ∈ D 1,2 and u(t) ∈ D 1,2 , 0 ≤ t ≤ T and let ϕ be a real continuous function on R. Then,

Main Results
First, we give the equi-step time division; assume ∆t = T/N, t n = n∆t for n = 0, 1, 2, . . . , N. For simple representation, we assume X n,i Then, we propose the following weak order 2.0 numerical scheme for solving SDEsMS.
i, X n ) with l(i, e) = j − i f or e ∈ ∆ ij and l(i, e) = 0 f or e / ∈ ∆ ij .
Remark 1. For the multiple stochastic integral in Scheme 1, we simulate t n+1 t n E L −1 e a i kÑ (de, dt) and t n+1 t n t t n E L −1 e b i k N(de, ds)dW t . by the following steps.
Step 1. Mark a discrete Markov chain i ∈ S from the definition of Markov process.
Step 2. Generate ∆N n and τ m , where ∆N n is the number of jumps and τ m is the switching time.
Step 4. Solve those multiple stochastic integrals by the following equations.

Convergence Theorem
In this section, we prove our new weak order 2.0 scheme in the condition of local weak convergence. We should firstly define weak convergence before proving the theorem.
with d ∈ N + , that is X N converges weakly with order β to X T as ∆t → 0. Theorem 1. (Local weak convergence) Suppose X n,i t n+1 and X n+1,i (1 ≤ n ≤ N) satisfy Equations (11) and (17), respectively. If Assumption 1 holds and X n,i t n+1 , X n+1,i ∈ D 2,4 , g ∈ C 4 p (R q , R) is a smooth function, we have where C is a generic constant depending only on K, the maximum state difference max Proof. For ease of proof, using multi-dimensional Taylor formula, we have is the kth component of explicit solution X n,i t n+1 . To simplify the notions, we assume r t n = i, a i k = a k (t n , i, X n ), b i k = b k (t n , i, X n ), h n,i k,α = h k,α (t n , i, X n ). By Itô-Taylor expansion, we can get the explicit solution and we have the numerical solution X n+1,i subtracting (17) from (16), we deduce I α h k,α (·, i · , X n,i · ) t n ,t n+1 .
Using the duality formula and Equation (18), we have where By taking Malliavin derivative with respect to X n+1,i k , we obtain for t n < t ≤ t n+1 , which by combining chain rule yields for t n < t ≤ t n+1 , where Φ 1 (t n , X n , i, ∆t, ∆W n , ∆Ñ n ) is a function not depending only on time t. Since for t n < t ≤ t n+1 . Similarly, for t n < t ≤ t n+1 , by chain rule we can obtain that Then, by Lemma 2 and Inequality (27), we have Using Lemma 2 (chain rule and duality formula), we have Combining Equation (19), J n,i 11 = 0, and Inequalities (28) and (29), we obtain that where C is a constant depending on K, the maximum state difference max  = (1, 1, 1), applying the Itô isometry formula, we have For α = (1, 1, 1), we similarly obtain ∑ α∈A 3 α =(1,1,1) Combining Inequalities (31) and (32), we have From Inequalities (30) and (33), we finally obtain The proof is completed.

Remark 2.
If Assumption 1 holds, under the conditions X n,i t n+1 , X n+1,i ∈ D 2,4 and g ∈ C 4 p (R q , R) we can obtain order 3.0 local weak convergence rate of Scheme 1. However, under a weaker regularity condition on the coefficients of a, b, we can only obtain lower order of local weak convergence of Scheme 1. For example, if X n,i t n+1 , X n+1,i ∈ D 1,2 and g ∈ C 2 p (R q , R), in Theorem 1, we can deduce

Numerical Experiments
In this section, before studying numerical examples, we simply simulate the generation of Markov chains and Brownian motions. Let the state space Markov chain r t be in S = {1, 2, 3} with a transition probability matrix as  In our numerical experiments, we consider two one-dimensional SDEwMS. We choose the sample number N sp = 5000, where N sp is the number of sample paths. The errors of global weak convergence can be measured by: and the average errors of local weak convergence can be measured by: where N = 1/∆t, ∆t are 2 −3 , 2 −4 , 2 −5 , 2 −6 , 2 −7 . We let ϕ(X j i ) = sin(X j i ). Assume X j i and X i,t j are the numerical solution and explicit solutions at the time t j , respectively, where j ∈ {1, 2, ..., N}.
Example 1. We consider the Ornstein-Uhlenbeck (O-U) process for SDEwMS: where W t is a one-dimensional Brownian motion. The Markov chain r t is in S = {1, 2, 3} and W t and r t are assumed to be independent. The group coefficients f and g are given by It is well known that Equation (34) has the explicit solution In Table 1, we set the time t ∈ [0, 1], where the terminal time is T = 1. CR stands for the convergence rate with respect to time step ∆t. We compute global errors and average local errors (Avg.local errors) of the new scheme. We can get that global convergence rate (Glo.CR) and average local convergence rate (Avg.local CR) have the order 2.0 and order 3.0, respectively. For more intuitive display, we show the numerical results of the local and global errors in Figure 2 (left). In Table 2, we simply utilize the Euler scheme and Runge-Kutta scheme to make a comparative experiment with the new scheme. It is easily to know that the new scheme has the most precision, but it also takes longer CPU time than Runge-Kutta scheme. For more intuitive display, the CPU Time and result of errors of all schemes are shown in Figure 2 (right). In addition, Figure 3 shows the mean-square stability of the new scheme with three kinds of time steps.

Example 2.
We consider the following geometrical model for SDEwMS: where W t is a one-dimensional Brownian motion. The Markov chain r t is in S = {1, 2} and W t and r t are assumed to be independent. The coefficients f and g are given by f (1) = 1.72, g(1) = 0.35, f (2) = 1.6, g(2) = 0.2. It is well known that Equation (35) has the explicit solution In Table 3, we find that Example 2 can also obtain that Glo.CR and Avg.local CR have the order 2.0 and order 3.0 of the new scheme. In Figure 4 (left), we clearly show that the global errors are different between a single state and switching states. In Table 4, for all three schemes, we fix a rational upper bounds of state values and find that the convergence rate of the new scheme will tend to order 1.0 when all state values become bigger. Meanwhile, the convergence rates of Milstein scheme and Euler scheme stay stable at order 1.0. In Table 5, for all three schemes, we fix the state difference and find that the convergence rate of the new scheme decline rapidly when all state values become bigger. The convergence rate of the new scheme is even negative when the state values are extremely big. Meanwhile, the convergence rates of Milstein scheme and Euler scheme are very low when the state values are extremely big. For more intuitive display, Figure 4 (right) shows the change process of convergence rate with the change times of state values (CTSV).

Conclusions
In this paper, we propose a new weak order 2.0 scheme for solving SDEwMS. By using trapezoidal rule and the integration-by-parts formula, we theoretically prove the new scheme has order 3.0 convergence under local weak condition. Two numerical experiments are given to confirm theoretical results. In addition, the first numerical experiment is also given to compare the New Scheme with other schemes, such as Euler scheme and Runge-Kutta scheme on convergence rate and accuracy, and the second numerical experiment is also given to explain that the change of Markov chain has some effects on convergence rate of all schemes. According to above experiments, we can obtain that the new scheme has the most precision of all schemes but costs longer time. Moreover, we also find that the maximum state difference and the upper bounds of state values have some effects on all schemes.