New Simplified High-Order Schemes for Solving SDEs with Markovian Switching Driven by Pure Jumps

: New high-order weak schemes are proposed and simplified to solve stochastic differential equations with Markovian switching driven by pure jumps (PJ-SDEwMs). Using Malliavin calculus theory, it is rigorously proven that the new numerical schemes can achieve a high-order convergence rate. Some numerical experiments are provided to show the efficiency and accuracy.


Introduction
Let Ω, F , {F t } t≥0 , P be a complete probability space with a filtration {F t } t≥0 gen- erated by a Poisson process.In this paper, we mainly study the second-order weak schemes of the following Equations (PJ-SDEwMs) on the probability space Ω, F , {F t } t≥0 , P : X t = X 0 + t 0 a(s, X s , r s )ds + t 0 E b(s, X s , r s , e) Ñ(de, ds) with initial value X 0 ∈ R d , where r s is Markov chain, Ñ(de, ds) is a compensated Poisson measure, and E = R d \ {0} is equipped with its Borel field E. The drift coefficient is denoted by a : [0, T] × R d × S → R d , and the jump diffusion coefficient is represented by Recently, the study of SDEs with Markovian switching driven by pure jumps has attracted increasing interest.PJ-SDEwMs can be seen as a generalization of the SDEs with jump.It is also possible to think of it as a generalization of SDEs with Markovian switching, of course.It is not only used in finance but also has a wide range of applications in control systems, bio-mathematics, chemistry and mechanics (see [1][2][3]).The authors [4] study mode coupling in a multimode step-index microstructured polymer optical fibers for potential sensing and communication applications.Ji and Chizeck [5] focused on the control problem for systems with continuous-time Markovian jump parameters.Mao [6] discussed the exponential stability for general nonlinear SDEwMs.Similar to SDEs with jump, it is difficult to obtain an explicit solution for SDEwMs.Therefore, we need effective schemes which are accurate and computationally convenient to approximate the true solutions.Yuan and Mao [7] discovered the convergence of the Euler-Maruyama scheme, which is used to obtain the stationary distribution of SDEwMs in [8].Mao and Yuan [9] gave the systematic presentation of the theory of SDEs with Markovian switching.Then, the existence and uniqueness of solutions for neutral SDEwMs were proven in [10] under non-Lipschitz conditions, and Euler approximate solutions were provided for solving SDEwMs.Common numerical schemes for solving SDEs with jumps or SDEwMs include Euler-Maruyama scheme [7,9,11,12], Milstein scheme [13,14], and jump-adapted scheme [15,16].The authors of [17] studied the balanced implicit numerical methods for solving SDEs driven by Poisson jumps.Zhou and Mamon [18] expanded three short-rate models, integrating the switching of economic regimes via a discrete-time finite-state Markov chain.
A high-order model incorporating drift and volatility modulation by a discrete-time weak Markov chain is introduced in [19].Yang, Yin and Li [20] utilized stochastic approximation techniques to analyze the stability of numerical solutions for jump diffusions with Markovian switching.Furthermore, the authors [21] study the convergence of SDEs differential equations containing delay and Markovian switching, and an Euler scheme for solving SDEwMs under non-Lipschitz conditions is given in [22].Given the application requirements in finance and related fields, there is an increasing interest in studying high-order numerical schemes for solving SDEs.For instance, Fan [23] developed a strong approximation order 1.5 scheme for solutions of SDEwMs.Liu and Li [24] proved the convergence of the weak stochastic Taylor scheme with an appropriate order.In [25], a new weak scheme for solving SDEwMs driven by Brownian motion was proposed and achieved the convergence rate of second-order.Additionally, the numerical results of several weak schemes are presented, with a focus on the second-order weak stochastic Taylor scheme and the extrapolation of the Euler scheme.We refer to the high-order numerical methods of solving SDEs in [24][25][26][27] to propose novel numerical schemes of pure jump SDEs, which can achieve a weak second-order convergence rate.
The primary contributions of this paper can be succinctly highlighted as follows: • For PJ-SDEwMs with mark-dependent jump coefficient b = b(t, X t , r t , e), we first propose Scheme 1 using Wagner-Platen expansion.However, Scheme 1 contains multiple stochastic integrals, which are not easily computed.Thus, to avoid the use of some double integrals, we propose another new Scheme 2, by employing the trapezoidal rule to approximate the following multiple stochastic integrals L 1 e a(s, X s , r s ) Ñ(de, ds) dt and Furthermore, we can use the definition of compound Poisson process to compute L e 1 b(s, X s , r s , e 2 ) Ñ(de 1 , ds) Ñ(de 2 , dt), which has no high-accuracy based on the truncation approximation.• Especially, for PJ-SDEwMs with mark-independent jump coefficient b = b(t, X t , r t ), we propose Scheme 3 by using the trapezoidal rule and duality formula, which does not involve multiple stochastic integrals.Moreover, Scheme 3 is not a special case of Scheme 2. Using Malliavin calculus theory, it is strictly proven that Scheme 3 has a local weak order-3.0convergence rate.The greatest state difference and the upper bound of the state value are connected to the convergence rate.

•
The convergence and stability results of Schemes 2 and 3 are validated through numerical experiments, which are also compared with the Euler scheme to verify its effectiveness and accuracy.Scheme 3 is simpler and faster than Scheme 2 in the case of mark-independent PJ-SDEwMs.
The following is a list of some notations to be used later: In Section 2, we give the introduction of fundamental concepts, encompassing the Markov Chain, Itô-Taylor expansion, and Malliavin stochastic calculus which include duality formula and chain rule.Section 3 presents our novel weak second-order numerical schemes, accompanied by a rigorous proof establishing their local weak convergence order of 3.0.In Section 4, we give the practical application of our proposed new schemes, where we present some numerical examples to validate the effectiveness and accuracy of them.The paper concludes with Section 5, providing a succinct summary of our work.
The following notations are listed for future reference: ) is the set of functions which have at most polynomial growth.♦ C is a generic constant depending only on the upper bounds of derivatives of a, b , g and the largest state difference.

Markov Chain
On the probability space Ω, F , {F t } t≥0 , P , we assume that {r t , t ≥ 0} is a rightcontinuous Markov chain and takes values in a finite state space S = {1, 2, . . ., M} with generator where δ > 0, q ij ≥ 0 and, for i ̸ = j, q ii = − ∑ i̸ =j q ij .Let E = R\{0} be the mark set equipped with its Borel field B(E ).Now, on E × [0, T], we consider a given intensity measure of the form λ(de) := γ(e)de with kernel function γ(e) ≥ 0 for all e ∈ E and λ(de) := 0 for e / ∈ E , and suppose that the total intensity λ E := E γ(e)de < ∞.Moreover, dr t = E h(r t− , e)N(de, dt) with h(i, e) = j − i for e ∈ ∆ ij and h(i, e) = 0 for e / ∈ ∆ ij , which the intervals ∆ ij have length q ij , that is and so on (see [9]).

Malliavin Stochastic Calculus
Suppose the operator D t,e is the Malliavin derivative of Poisson process at (t, e).A random variable U is Malliavin differentiable if and only if U ∈ D l,m .Here, the stochastic Sobolev spaces D l,m consist of all F T −measurable U ∈ L 2 (P) with the norm where the Malliavin derivative D α s 1 ...s l ,e is defined as with especially D 0 s j ,e = 1 for 1 ≤ j ≤ l.

Main Results
First, we consider a regular time uniform discretization: , which is the k−th component of Using the classical Wagner-Platen expansion, we can derive the following Scheme 1 for solving SDEwMs with mark-dependent jump coefficient.

Scheme 1 (Wagner-Platen expansion). Given the initial condition X
From the Wagner-Platen expansion and trapezoidal rule we obtain where the truncation error with Here, we write a n,i k for a k (t n , X n,i , i) and b n,i k,e for b k (t n , X n,i , i, e).Then, by Equation ( 12), we propose the following second-order new scheme.

Scheme 2. Given the initial condition X
with ∆t = t n+1 − t n .
Remark 1.If the jump coefficient function b = b(t, X t , r t , e), we generate the compound Poisson process and the multiple stochastic integral in Scheme 2 can be computed by where the pairs (τ k , ξ k ) of k-th jump time and marks are independent uniformly distributed in (N In the special case of a mark-independent jump coefficient b(t, X t , r t , e) = b(t, X t , r t ), we use the following discrete-time approximation with the truncation error Based on Equation ( 16), we propose the following simplified scheme for solving markindependent PJ-SDEwMs. where Remark 2. Here, γ(e) is a kernel function, which may be symmetric, i.e., γ(e) = γ(−e); or non-symmetric, i.e., or singular, i.e.

Local Weak Convergence Theorems
In this section, using Malliavin stochastic analysis and the Wagner-Platen expansion, we rigorously prove and obtain the local weak order-3.0convergence of Schemes 1-3.
Theorem 1. (Local weak convergence) Suppose that X n,i t n+1 and X n+1,i (0 where m ∈ N + is a generic constant, which can vary from line to line. Proof of Theorem 1: Subtracting Equation (18) from Equation (12) yields Then, by the mean value formula of integrals and duality formula, we have with Under the conditions of this theorem, we finally obtain the inequality (21).□ Theorem 2. Suppose that X n,i t n+1 and X n+1,i (0 ≤ n ≤ N − 1) satisfy Equation (12) and Scheme 2, respectively.If the functions a, b where m ∈ N + is a generic constant, which can vary from line to line.
Proof of Theorem 2. Using the multi-dimensional Taylor formula, for ease of proof, we assume where is the k-th component of explicit solution X n,i t n+1 .Then, it follows from the Itô-Taylor expansion that which by the fact t n E L 0 b i,e k Ñ(de, ds) k Ñ(de 1 , ds 1 ) Ñ(de 2 , ds 2 ). ( It follows from the definition of the operator L 1 e 1 that By the duality formula in Lemma 3 and from Equations ( 29) and (30), we deduce where By taking Malliavin derivative with respect to X n+1,i k , we obtain for t n < t ≤ t n+1 and e 2 ∈ E .Furthermore, the chain rule (10) gives Then, it follows from Taylor formula and Lemma 5 that where Y n+1 i,e 2 = X n+1,i + b i,e 2 + 1 2 ∆t L 0 b i,e 2 + ∑ j∈S L i,j a i and By the duality formula, we have which by using the fact (37) Note that and we deduce from Lemma 3 that Using the duality formula in Lemma 3, we conclude that Combining the inequalities (37)-( 40), we obtain For α = (1, 1, 1), applying the Itô isometry Formula (4), we have For α ̸ = (1, 1, 1), we obtain By the inequalities ( 42) and ( 43), we deduce From the inequalities (41) and (44), we finally obtain □ Theorem 3. Assume that X n,i t n+1 and X n+1,i (0 ≤ n ≤ N − 1), respectively, satisfy Equation ( 16) and Scheme 3.
where m ∈ N + is a generic constant, which could change line by line.
Proof of Theorem 3: Using multi-dimensional Taylor formula, to make the proof easier, we have where (47) is the k-th component of explicit solution X n,i t n+1 .Then, from the Itô-Taylor expansion, we can obtain which yields By the duality formula in Lemma 3 and from Equation (49), we deduce where For t n < s ≤ t ≤ t n+1 , by using Malliavin derivative in relation to X n+1,i k , we get which by combining chain rule (10) gives where the functions Φ(t n , X n , ∆t, ∆ Ñn ) and Ψ(t n , X n , ∆t, ∆ Ñn ) do not depending only on t and e.Furthermore, from Lemma 3, using the notation λ E := E λ(de) = E γ(e)de, we have which by using the fact we deduce Using Itô's formula, we can obtain by the duality formula we have Applying the duality formula in Lemma 3, we have Combining Equation (50), ϵ 1 g,k = ϵ 2 g,k = ϵ 3 g,k = 0, and inequality (57), we have For α = (1, 1, 1), by using the Itô isometry Formula (4), we obtain For α ̸ = (1, 1, 1), we obtain Combining inequalities (59) and (60), we have (61) From inequalities (58) and (61), we finally obtain □

Numerical Experiments
Assume that the state space Markov chain r t is in S = {1, 2, 3}, and the transition probability matrix is We choose N sp = 5000 as the sample size for our numerical experiments, where N sp is the total number of sample pathways.We can measure the average errors of local weak convergence and the errors of global weak convergence as follows: .We let φ X j i = sin X j i .Let X n,i represent numerical solution and X t n represent explicit solutions at the time t n , where j ∈ {1, 2, ..., N}.Now, we give three numerical examples, including mark-dependent PJ-SDEwMs, markindependent PJ-SDEwMs (Ornstein-Uhlenbeck type) and PJ-SDEwMs (geometrical type).
Example 1.We consider the following mark-dependent PJ-SDEwMs: where µ is a constant and the Markov chain r t is in S = {1, 2, 3}.The group coefficients g are given by g(1) = 0.35, g(2) = 0.3, g(3) = 0.25.We use the Itô formula to obtain the explicit solution of Equation (62) which is t n E e Ñ(de, ds) = Then, for solving the PJ-SDEwMs (62), we have Scheme 2 with the form t n E e Ñ(de, dt) which by the computation of compound Poisson process (see [27]) and the trapezoidal rule gives where E e 2 λ(de 2 ) = 1 0 e 2 2 de 2 = 1 3 , We use CR to represent the rate of convergence over the time step ∆t.To evaluate the performance of Scheme 2, we calculate their global errors and average local errors.It is gratifying to find that the global convergence rate (Glo.CR) has order 2.0, while the average local convergence rate (Avg.localCR) has order 3.0.This means that Scheme 2 has excellent convergence properties and shows the great accuracy in numerical calculations (see Table 1).At the same time, we draw trajectories of numerical and analytical solutions in Figure 1.By comparing different state values, it was found that regardless of how the state values change, the simulation efficacy demonstrated by Scheme 2 remains highly favorable.Example 2. Consider the Ornstein-Uhlenbeck (O-U) PJ-SDEwMs as follows: where µ is a constant and the Markov chain r t is in S = {1, 2}, and the group coefficients g are given by g(1) = 0.3 and g(2) = 0.25.It is evident that Equation (63) possesses an explicit solution: e −µ(t−s) g(r s ) Ñ(de, ds).
According to the evaluation results presented in Table 2, a detailed analysis was conducted regarding the performance of Schemes 2 and 3.The evaluation included the calculation of their global error and average local error.Interestingly, both Schemes 2 and 3 demonstrated a second-order convergence rate during the global convergence process, indicating a relatively rapid approach towards the optimal solution.Furthermore, in terms of local convergence performance, both schemes exhibited a higher level of convergence with a third-order average convergence rate.This finding suggests that both Schemes 2 and 3 exhibit favorable performance.In Table 3, we compare the global errors and convergence rate of the Euler scheme, Schemes 2 and 3.It is obvious that Scheme 3 makes it simpler and more convenient for us to calculate in comparison to Scheme 2. Equally gratifying is that Scheme 3 greatly reduces the computing time.Scheme 3 takes 3.597104 seconds, while Scheme 2 takes 10.024912 seconds.Example 3. We consider the following PJ-SDEwMs with mark-independent jump coefficient: where the Markov chain r t is in S = {1, 2}, and Ñ(de, dt) represents a compensated Poisson measure in one dimension.Assuming that N t and r t are independent, with the group coefficients f and g provided by We utilize the Itô formula to derive the explicit solution for Equation (64) as The time t ∈ [0, T] with T = 1 is set in Table 4, and we use Scheme 3 to solve the PJ-SDEwMs (64).Scheme 3 has the second-order global convergence rate (Glo.CR) and the third-order average local convergence rate (Avg.localCR).In Tables 5 and 6, we compare the Euler Scheme with Scheme 3 with different state values.For clearer display, we draw two pictures (errors of the new scheme and CPU times) according to the two tables in Figure 1.In Figure 2 (left), we clearly demonstrate that there are differences in global errors between individual states and transition states.For a more intuitive presentation, Figure 3 (right) illustrates the fluctuation of convergence rates as the number of times changes, showcasing the process of convergence rate variation with the change in state value times (CTSV).These visualizations allow us to gain a clearer understanding of the impact of state changes on convergence rates.The observation enhances our comprehension of the dynamic nature of the system.

Discussion
In this work, we mainly study stochastic differential equations with Markovian switching driven by pure jumps (PJ-SDEwMs) and give three numerical schemes.In general, PJ-SDEwMs contains mark-dependent jump coefficient b = b(t, X t , r t , e), which we can solve using Schemes 1 and 2. Compared to the Itô-Taylor expansion scheme (Scheme 1), Scheme 2 is easier to calculate by using the trapezoidal rule to approximate the following multiple stochastic integrals: L 1 e a(s, X s , r s ) Ñ(de, ds) dt and t n+1 t n E t t n L 0 b(s, X s , r s , e)ds Ñ(de, dt).
In particular, PJ-SDEwMs contains mark-independent jump coefficient b = b(t, X t , r t ), we can compute it in Schemes 2 and 3.Because multidimensional random integrals are avoided, Scheme 3 is simpler and more convenient (Example 2 demonstrates this very well).In addition, by using Malliavin calculus theory, we strictly proved that the proposed new schemes have local weak order-3.0convergence rates.However, through Example 3, we find that as the upper bound of the state values gradually increases, the simulation effect of Scheme 3 is still good, but the convergence rates will gradually decrease.

Conclusions
In this paper, we propose three new weak second-order numerical schemes to solve stochastic differential equations with Markovian switching driven by pure jumps.By using the Malliavin stochastic analysis method, the new schemes are strictly analyzed theoretically, and the second-order convergence rate is proven.Finally, the correctness and effectiveness of the second-order schemes are verified by three numerical experiments.In addition, we find that as the upper bound of the state values increases, the global convergence rate of Scheme 3 gradually decreases to the first-order.Besides, the maximum state difference and the variation of Markov chains have a certain impact on the convergence rate.

Figure 2 .
Figure 2. (Left) The results of the global errors and the average local errors for Scheme 3 in Example 3. (Right) The correlations for global errors and CPU time of all schemes.

Figure 3 .
Figure 3. (Left) The convergence rates of different states.(Right) The correlations for global convergence rates and the variation of state values.

Table 1 .
The results of convergence rates and errors for Scheme 2 in Example 1.

Table 2 .
The results of convergence rates and errors for Scheme 3 in Example 2.

Table 3 .
The results of global convergence rates and errors for three schemes in Example 2.

Table 4 .
The results of convergence rates and errors for Scheme 3 in Example 3.

Table 5 .
Global convergence rates of multiple groups of different state values for two schemes with X 0 = 0.5.