A New Simpliﬁed Weak Second-Order Scheme for Solving Stochastic Differential Equations with Jumps

: In this paper, we propose a new weak second-order numerical scheme for solving stochastic differential equations with jumps. By using trapezoidal rule and the integration-by-parts formula of Malliavin calculus, we theoretically prove that the numerical scheme has second-order convergence rate. To demonstrate the effectiveness and the second-order convergence rate, three numerical experiments are given.


Introduction
Many studies are interested in the Lévy process such as the classical stochastic differential equations (SDEs). SDEs have been widely applied to economics and finance [1]. Nowadays, stochastic differential equations with jumps (SDEwJs) are paid more attention by many scholars (see [2][3][4]). As is known, the geometrical and Ornstein-Uhlenbeck (O-U) models are very important in finance, e.g., American option pricing about asset price follows the geometrical model and can be a jump diffusion process (see [4,5]). Chockalingam and Muthuraman [6] studied pricing options when asset prices jump. Li and Zhang [7] employed additive subordination to construct pure jump models for volatility index. Dang, Nguyen and Sewell [8] considered the pricing of Asian options under models that incorporate both regime-switching and jump-diffusions. The O-U process is used to calculate the short term interest rate, and it can also used to account for the mean reversion of prices in modeling commodity processes (see [9]). In this paper, we consider the SDEwJs of the following form: with initial value X 0 ∈ R q . Here, E = R q \ {0} is equipped with its Borel field E,Ñ(de, ds) is compensated poisson measure, and the operators a : [0, t] × R q → R q and b : [0, t] × E × R q → R q are the drift coefficient and the jump coefficient, respectively. Qualitative theory of the existence and uniqueness of the solution for SDEwJs is studied in [10]. Generally, most of SDEwJs do not have explicit solutions and hence require numerical solutions. It is important to apply appropriate numerical schemes in practice. About numerical solutions of SDEwJs, Platen and Bruti-Liberati [9] systematically introduced the weak schemes and strong schemes with respect to SDEwJs. There are Euler-Maruyama schemes [11], Milstein schemes [12], and jump-adapted schemes [13,14] for solving SDEwJs. Higham and Kloeden [15] presented and analyzed two implicit methods for SDEs and proved in [16] that implicit methods share the same finite time strong convergence rate as the explicit Euler scheme. Hu and Gan [17] studied numerical stability of balanced methods for solving SDEwJs.
Mil'shtein [18] studied the second-order accuracy integration of stochastic differential equations. Similarly, for the higher order numerical methods of SDEwJs, Gardoń [19] proposed the order 1.5 approximation for solutions of jump-diffusion equations. There are [9] the weak or strong order 2.0 Taylor schemes and jump-adapted order 2.0 Taylor schemes. Moreover, the higher order of Runge-Kutta methods for jump-diffusion differential equations can be found in [20]. However, these numerical schemes include multiple stochastic integrals, which are difficult to accurately compute and simulate. Therefore, the simplified weak order 2.0 scheme does not contain the multiple stochastic integrals; particularly, it is convenient and easily calculated for solving SDEwJs in multi-dimensional case. Liu and Li [21] proposed the weak stochastic Taylor order 2.0 (WST2) scheme of SDEwJs with three-point distribution random variables, obtained the convergence rate of the Itô-Taylor scheme, i.e. using the Itô-Taylor expansion and product rule.
In this paper, we propose a new simplified numerical scheme by using the compound Poisson process with a sequence of pairs (τ k , ξ k ) of jump time τ k and marks ξ k , where (τ k , ξ k ) are uniformly distributed in the square [0, 1] × [0, 1]. We rigorously prove and obtain its weak second-order convergence rate by using the Malliavin stochastic analysis.
The important contributions of this paper can now be highlighted as follows: • Our new scheme is simplified without multiple integrals, employs the compound Poisson process with uniformly distribution, and can be easily used to compute the pure jump stochastic models. • Using a new proof method, i.e. the trapezoidal rule and integration-by-parts formula of Malliavin calculus theory, we rigorously prove that the new scheme has secondorder convergence rate. • We give three experiments to demonstrate the effectiveness and the accuracy of our new scheme, which are consistent with our theoretical results.
Some notations to be used later are listed as follows: 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 it can be different from line to line.
The paper is organized as follows. In Section 2, we introduce some needed preliminaries. In Section 3, we propose a new numerical scheme and obtain its second-order convergence error estimate. In Section 4, three numerical examples are given to verify the theoretical results.

Preliminaries
Let there be a filtered probability space (Ω, F , P) with the filtration {F t } t≥0 , and assume the filtration {F t } t≥0 satisfies the usual hypotheses of completeness, i.e., F 0 contains all sets of P-measure zero, and maintains right continuity, i.e., F t = F t+ . Moreover, the filtration is assumed to be Poisson random measure N(A, t) on E × [0, t], where E = R \ {0} is equipped with its Borel field E. Let λ(de) be a σ-finite measure on (E , E) satisfying λ(de) := γ(e)de with γ(e) ≥ 0 for all e ∈ E and λ(de) := 0 for e / ∈ E , and λ E = E γ(e)de < ∞. The compensated Poisson random measure can be represented bỹ N(de, dt) = N(de, dt) − λ(de)dt.
For Equation (1) and x ∈ R q , we define operators L 0 ϕ, Then, the Itô formula can be presented as:

The Weak Schemes for Solving SDEwJs
The simple numerical method for the approximate solution of the SDEwJs is the Euler scheme. Thus, we introduce the Euler scheme for solving SDEwJs. Thanks to Itô-Taylor expansion, we can obtain weak convergence schemes for solving SDEsJs as follows.
When accuracy and efficiency are required, it is important to construct numerical methods with higher order of convergence. With respect to weak second-order scheme, Liu and Li [21] proposed the weak stochastic Taylor (WST2) scheme and approximated the Poisson jump measure by using three point distribution.
and T n obeys three-point distribution, which has the following probability law: Then, Buckwar and G. Riedler [20] proposed the Runge-Kutta second-order implicit scheme for solving SDEwJs.
In addition, the discrete-time approximations considered are divided into regular and jump-adapted schemes. Regular schemes employ time discretizations that do not include the jump times of the Poisson jump measure and Jump-adapted time discretizations include these jump times. Platen and Bruti-Liberati [9] proposed the following jump-adapted weak second-order scheme.

Scheme 4. (Jump-adapted weak 2-order Scheme)
Assume the initial condition X 0 . For 0 ≤ n ≤ M − 1, we consider a jump-adapted time discretization 0 = t 0 < t 1 < · · · < t M = T, which is constructed by a superposition of jump times {τ 1 , τ 2 , · · ·} of the compensated Poisson measureÑ(de, dt) and equidistant time discretization with step size ∆t, as given in Section 2.1. For convenience, we set X t n = X n in this section and define in the almost sure limit. We have the weak second-order scheme where ∆Ñ n =Ñ t n+1 −Ñ t n and ∆t = t n+1 − t n (0 ≤ n ≤ M − 1).

Remark 1.
Comparing the above methods to solve SDEwJs, the coefficients of Runge-Kutta implicit scheme for solving SDEwJs are not easily determined; the WST2 scheme in [21] needs to consider the probability law of three point distribution and Poisson jump marks, and Jump-adapted weak second-order scheme needs to consider time discretizations that include these jump times.

Malliavin Stochastic 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 the operator D k be the Malliavin derivative of order k with respect to the lévy process. 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 Malliavin derivative of Poisson process (see [22] for details) has the following definition:  By induction, it follows that, if F ∈ D 1,2 , then D t,e (F n ) = (F + D t,e F) n − F n .
Let F ∈ D 1,2 and ϕ be a real continuous function on E . Suppose ϕ(F) ∈ L 2 (P) and ϕ(F + D t,e F) ∈ L 2 (P ×Ñ × λ). Then, ϕ(F) ∈ D 1,2 and D t,e ϕ(F) = ϕ(F + D t,e F) − ϕ(F). (13) According to Lemma 3, by chain rule (13), taking Malliavin derivatives with respect to ∆Ñ n and (∆Ñ n ) 2 , we obtain where Here, we write a k for a k (t n , X n ) and b k for b k (t n , X n ). Then, we have the following scheme.

Remark 2.
Unlike jump-adapted scheme, Scheme 5 employs regular time discretization, and it is a simplified scheme, which does not involve multiple stochastic integrals. If the jump coefficient function b = b(t, X t , e), we generate the compound Poisson process for simulating compensated Poisson process, where the pairs (τ k , ξ k ) are uniformly distributed in the square [0, 1] × [0, 1]. We found that Scheme 5 has a little impact of the random generation of the compound Poisson process in Example 3. From the numerical experiments, comparing with the three point distribution in WST2 Scheme in [17], our scheme can be implemented in programming and costs less run time. Scheme 5 is an explicit scheme different from Runge-Kutta implicit scheme in [20]. Moreover, in the next subsection, using the integration-by-parts formula of Malliavin calculus, we theoretically prove that our scheme has second-order convergence rate.

Local Weak Convergence Theorem
In this section, using Malliavin stochastic analysis and Itô-Taylor expansion, we obtain the local weak order-2.0 convergence of our new weak second-order scheme.
Theorem 1. (Local weak convergence) Suppose X n t n+1 and X n+1 (0 ≤ n ≤ N − 1) satisfy Equation (1) and Scheme 5, respectively. Under Assumption 1, if the functions a, b ∈ C k p (R q , R) and g ∈ C 3 b , then where r ∈ N + is a generic constant which can vary from line to line.
Proof. Using multi-dimensional Taylor formula, for ease of proof, we have By the duality formula in Lemma 2 and from Equation (19), we deduce where By taking Malliavin deriavtive with respect to X n+1 k we obtain D t,e X n+1 which by combining chain rule (13) gives where φ(t n , X n , ∆t, ∆Ñ n ) is a function not depending only on t.

Numerical Experiments
In this section, we consider three one-dimensional stochastic differential equations with pure jump models to verify the accuracy and effectiveness of the theoretical results. We compare our new scheme in precision as well as the time of computing with other schemes. Assume T = 1 is terminal time, N sp = 5000 is the number of sample paths in the numerical experiment, and N ∈ {2 3 , 2 4 , 2 5 , 2 6 , 2 7 }. The errors of global weak convergence can be measured by e global ∆t and the errors of local weak convergence can be measured by where ϕ(X j i ) = sin(X j i ) and X j i and X i,t j are the numerical solution and explicit solution at the time t j , respectively. Example 3. Consider the O-U process with pure jump: whereÑ(de, ds) is one-dimensional compensated poisson process with γ(e) = 1, E = [0, 1] and λ E = 1 0 de = 1. Let a = 1.5, h = 0.01, and e ∈ U(0, 1). Applying Itô formula to e at X t , Equation (29) has the explicit solution X t = X 0 e −at + he −at t 0 E e as eÑ(de, ds).
For the model (29), by using Itô Taylor expansion, the new Scheme 5 is

Scheme 5 becomes
where τ k is the kth jump time, ξ k is the kth jump mark, and the pairs (τ k , ξ k ) are uniformly distributed in the square [0, 1] × [0, 1]. In this example, the right of Figure 1 shows that, even if the jumps number ∆N n increase in [t n , t n+1 ], the convergence rate of Scheme 5 can still achieve second order. Meanwhile, the left of Figure 1 presents that the global errors may have some difference with different ∆N n ; a higher number of jumps accompanies less global error in numerical computing in [t n , t n+1 ]. Then, we conclude that there is a little impact of the random generation of the compound Poisson process.
whereÑ(de, ds) is a one-dimensional compensated Poisson process with λ E = 1. Let a = 0.8 and h = 0.01; the equation has the explicit solution by using Itô formula For the model (30), the new Scheme 5 is X n+1 =X n + aX n ∆t + hX n ∆Ñ n + 1 2 a 2 X n (∆t) 2 + ahX n ∆Ñ n ∆t + 1 2 In this pure jump geometrical model, we apply Scheme 5 to compute E[X N ] using 5000 sample trajectories. In Table 1, we obtain the scheme errors with ∆t = {2 −3 , 2 −4 , 2 −5 , 2 −6 , 2 −7 }. In addition, we show the global errors and average local errors of the scheme have the convergence rates of global second-order and average local third-order, respectively. Meanwhile, we intuitively display the accuracy of Scheme 5 from the left of Figure 2, where the two coordinate axes are log 2 ∆t and log 2 |ϕ(X t N ) − ϕ(X N )|. On the other hand, from the right of Figure 2, we verify the result of mean-square convergence stability with To illustrate the advantages of Scheme 5 in computational efficiency, we also apply Euler Scheme, Runge-Kutta Scheme, WST2 Scheme, and Jump-adapted Scheme to solve Equation (30). The global errors and convergence rates of the schemes are displayed in Table 2 and the left of Figure 3 displays the sample global errors and the corresponding CPU time of the Euler Scheme, Scheme 5, Runge-Kutta Scheme, WST2 Scheme, and Jump-adapted Scheme with X 0 = 0.5, a = 1.5, h = 0.001.

Conclusions
In this paper, we propose a new simplified weak second-order numerical scheme for solving stochastic differential equations with jumps. By using trapezoidal rule and the integration-by-parts formula of Malliavin calculus, we theoretically prove that the numerical scheme has second-order convergence rate. In numerical experiments, we used three examples to verify the convergence of our scheme. Meanwhile, we compared the computational time, programming complexity and precision of our scheme with other schemes, such as Euler Scheme, Runge-Kutta Scheme, WST2 Scheme and Jump-adapted Scheme and found our new scheme can be easily computed and consume less run time.