Second-order numerical schemes for decoupled forward-backward stochastic differential equations with jumps
aa r X i v : . [ m a t h . NA ] J u l Second-order numerical schemes for decoupled forward-backwardstochastic differential equations with jumps
Weidong Zhao
Department of Mathematics, Shandong University, Jinan, Shandong 250100, ChinaEmail: [email protected]
Wei Zhang
Department of Mathematics, Beijing University of Technology, Beijing 100022, ChinaEmail:[email protected]
Guannan Zhang
Department of Computational and Applied Mathematics, Oak Ridge National Laboratory, Oak Ridge,Tennessee 37831, United StatesEmail: [email protected]
Abstract
We propose new numerical schemes for decoupled forward-backward stochastic differ-ential equations (FBSDEs) with jumps, where the stochastic dynamics are driven by a d -dimensional Brownian motion and an independent compensated Poisson random measure.A semi-discrete scheme is developed for discrete time approximation, which is constitutedby a classic scheme for the forward SDE [17, 25] and a novel scheme for the backwardSDE. Under some reasonable regularity conditions, we prove that the semi-discrete schemecan achieve second-order convergence in approximating the FBSDEs of interest; and suchconvergence rate does not require jump-adapted temporal discretization. Next, to add inspatial discretization, a fully discrete scheme is developed by designing accurate quadraturerules for estimating the involved conditional mathematical expectations. Several numericalexamples are given to illustrate the effectiveness and the high accuracy of the proposedschemes. Mathematics subject classification:
Key words:
Decoupled FBSDEs with L`evy jumps, backward Kolmogorov equation, non-linear Feynman-Kac formula, second-order convergence, error estimates.
1. Introduction
In this work, we study numerical solution of decoupled forward-backward stochastic differ-ential equations (FBSDEs) with jumps, where the underlying stochastic jump processes arecharacterized by Poisson random measures. The term “decoupled” refers to the fact that theforward SDE is independent of the solution of the backward SDE. This work is motivated bya wide variety of applications offered by FBSDEs. In finance and insurance, FBSDEs-basedapproaches [23, 29] have gained a great attention by both academics and practitioners, becauseFBSDEs provide us a unified framework to describe the mathematical problems which arisein option pricing [13], portfolio hedging [14], market utility maximization [2] and risk mea-sures [24,26], etc. Moreover, in the presence of jump behaviors in many financial problems [25],L`evy jump processes have been incorporated into FBSDEs [11, 14], so as to accurately capture and properly interpret event-driven stochastic phenomena, such as corporate defaults, opera-tional failures, insured events, etc. In mathematics, one can relate FBSDEs with jumps to aclass of nonlinear partial integro-differential equations (PIDEs), based on the extension of thenonlinear Feynman-Kac theory studied in [1]. As such, FBSDEs become a powerful probabilistictechnique for studying analytical and numerical solutions and properties of the PIDEs, wherethe nonlocal integral operators of the PIDEs are characterized by Poisson random measures inthe framework of FBSDEs. In engineering science, a particular application of the PIDEs is tomodel anomalous diffusion [19], i.e., super-diffusion and sub-diffusion, that has been verified ex-perimentally to be present in various applications, e.g., contaminant transport in groundwaterand plasma physics. In this setting, FBSDEs-based probabilistic numerical schemes have beendeveloped in [31] to solve the governing PIDEs, which illustrated effectiveness of the FBSDEsmodels.There are many theoretical results on FBSDEs with jumps over the past two decades. Theexistence and uniqueness were proved by Tang and Li [29] for backward stochastic differentialequations with Poisson jumps and Lipschitzian coefficients, which was then extended, by Rongin [27], to the case of non-Lipschitzian coefficients. In [1], Barles, Buckdahn and Pardouxestablished a comparison theorem for decoupled FBSDEs with jumps as well as the link betweensuch FBSDEs and PIDEs, which generalized the results in [21, 22] to the case of a naturalfiltration associated with a Brownian motion and a Poisson random measure. After that, inthe context of FBSDEs with jumps, Øksendal and Sulem [20] established maximum principles,and Royer [28] introduced nonlinear expectations. For a general overview of related topics,see [9, 11] and the references therein.The obstacle of applying FBSDEs with jumps to real-world engineering and finance prob-lems results from the challenge of solving FBSDEs analytically or numerically. Since it is typi-cally difficult to obtain analytical solutions, numerical solutions are highly desired in practicalapplications. Numerical methods for FBSDEs without jumps have been well studied in the lit-erature [3,7,8,10,12,15,32–34,36], nevertheless, there are very few numerical schemes developedfor FBSDEs with jumps, and most of those schemes only focused on temporal discretization.For instance, a Picard’s iterative method was provided in [18], and numerical schemes of back-ward SDE were studied in [4, 5]. Due to the aforementioned applications of FBSDEs withjumps, it is of great significance to develop high-order temporal-spatial discretization schemesfor solving not only the FBSDEs but also the PIDEs and related engineering problems.In this paper, we propose novel numerical schemes for decoupled FBSDEs driven by a d -dimensional Brownian motion and an independent compensated Poisson random measure. Ingeneral, the approximation of the FBSDEs under consideration includes two steps, i.e., con-structing a semi-discrete scheme for temporal discretization, and extending it to a fully discretescheme by incorporating effective spatial discretization. By imposing appropriate regularityconditions on the coefficients, the generator and the terminal condition, we rigorously prove the second-order convergence of the semi-discrete scheme with respect to ∆ t . In spatial discretiza-tion, a carefully designed quadrature rule is critical to approximate all the involved conditionalmathematical expectations which are, in this case, multiple integrals with respect to both theBrownian motion and the Poisson random measure. The integrals with respect to the Brownianmotion is estimated by the Gauss-Hermite rule. For the integrals with respect to the Poissonrandom measure, we propose a general quadrature rule for the case that the jump componenthas finite activities. A specific form of the quadrature rule can be determined based on thetype of the underlying L`evy measure. For the numerical experiments in §
6, the L`evy measureis defined as a uniform distribution on bounded domains, so that Gauss-Legendre rule is anappropriate choice. Moreover, to avoid the explosion of the total number of quadrature pointswith the increase of time steps, we construct a piecewise Lagrange interpolating polynomialson a pre-determined spatial mesh, which are used to evaluate the integrand at all quadraturepoints.The main contributions of this paper are as follows: • propose a second-order discrete time approximation (semi-discrete) scheme for decoupledFBSDEs with jumps. • rigorously analyze the convergence rate of the proposed discrete time approximationscheme with respect to ∆ t . • propose a fully discrete scheme by developing new quadrature rules for estimating involvedconditional mathematical expectations.The outline of the paper is organized as follows. In §
2, we introduce the mathematicaldescription of the FBSDEs under consideration. In §
3, we propose the semi-discrete scheme,i.e., time discretization, for the FBSDEs of interest. Rigorous error analysis for the proposedsemi-discrete scheme is conducted in §
4. The fully discrete scheme for the case of Poissonrandom measures with finite activities is proposed in §
5. Numerical examples are given in § §
2. Preliminaries
Let (Ω , F , {F t } ≤ t ≤ T , P ) be a stochastic basis satisfying the usual hypotheses of complete-ness, i.e., F contains all the sets of P -measure zero and possesses right continuity, i.e., F t = F t + .The filtration {F t } ≤ t ≤ T is assumed to be generated by two mutually independent processes,i.e., one d -dimensional Brownian motion W t = ( W t , . . . , W dt ) ⊤ and one Poisson random measure µ ( A, t ) on E × [0 , T ] where E = R q \{ } is equipped with its Borel field E . The compensator of µ and the resulting compensated Poisson random measure are denoted by ν ( de, dt ) = λ ( de ) dt and˜ µ ( de, dt ) = µ ( de, dt ) − λ ( de ) dt , respectively, such that { ˜ µ ( A × [0 , t ]) = ( µ − ν )( A × [0 , t ]) } ≤ t ≤ T is a martingale for all A ∈ E . λ ( de ) is assumed to be a σ -finite measure on ( E, E ) satisfying Z E (1 ∧ | e | ) λ ( de ) < + ∞ , where | · | denotes the standard Euclidean norm in Euclidean spaces.In the probability space (Ω , F , {F t } ≤ t ≤ T , P ), we introduce the following forward-backwardstochastic differential equation with jumps X t = X + Z t b ( s, X s ) ds + Z t σ ( s, X s ) dW s + Z t Z E c ( s, X s − , e )˜ µ ( de, ds ) ,Y t = ξ + Z Tt f ( s, X s , Y s , Z s , Γ s ) ds − Z Tt Z s dW s − Z Tt Z E U s ( e )˜ µ ( de, ds ) , (2.1)where the quadruplet ( X t , Y t , Z t , U t ) is the unknown, b : [0 , T ] × R q → R q is referred to asthe drift coefficient, σ : [0 , T ] × R q → R q × d is referred to as the local diffusion coefficient, c : [0 , T ] × R q × E → R q is referred to as the jump coefficient, f : [0 , T ] × R q × R p × R p × d × R p → R p is referred to as the generator of the FBSDE, and the process Γ s is defined byΓ s = R E U s ( e ) η ( e ) λ ( de ) for a given bounded function η : E → R , i.e., sup e ∈ E | η ( e ) | < + ∞ . Theterminal condition ξ is an F T -measurable random vector in R q . A quadruplet ( X t , Y t , Z t , U t )is called an L -adapted solution if it is an {F t } -adapted, square integrable processes satisfyingthe FBSDEs in (2.1).Under standard assumptions on the given data b , σ , f , ϕ and c (see [1] for details), thereexists a unique solution ( Y t , Z t , U t ) ∈ S × L ( W ) × L (˜ µ ) for the backward SDE in (2.1),where S is the set of {F t } -adapted c`adl`ag processes { Y t , ≤ t ≤ T } such that k Y k S := E h(cid:0) sup ≤ t ≤ T | Y t | (cid:1) i < ∞ , L ( W ) the set of F t -progressively measurable q × d dimensionalprocesses { Z t , ≤ t ≤ T } such that k Z k L ( W ) := E hR T | Z t | dt i < ∞ , and L (˜ µ ) the set ofmappings U : Ω × [0 , T ] × E → R such that k U k L (˜ µ ) := E hR T R E U t ( e ) λ ( de ) dt i < ∞ . Now we introduce a class of nonlinear partial integro-differential equations (PIDEs) thatwill be related to the FBSDEs in (2.1) later. We consider the unique viscosity solution u ( t, x ) ∈C ([0 , T ] × R q ) of the following nonlinear PIDE, i.e., ∂u∂t ( t, x ) + e L [ u ]( t, x ) + f ( t, x, u, σ ∇ u, B [ u ]) = 0 , for ( t, x ) ∈ [0 , T ) × R q ,u ( T, x ) = ϕ ( x ) , for x ∈ R q , (2.2)where ϕ ( x ) is the terminal condition at the time t = T , e L is the second-order integral-differentialoperator of the form e L [ u ]( t, x ) = q X i =1 b i ( t, x ) ∂u∂x i ( t, x ) + 12 q X i,j =1 ( σσ ⊤ ) i,j ( t, x ) ∂ u∂x i ∂x j ( t, x )+ Z E u ( t, x + c ( t, x, e )) − u ( t, x ) − q X i =1 ∂u∂x i ( t, x ) c ( t, x, e ) ! λ ( de ) , (2.3)and B is an integral operator defined as B [ u ]( t, x ) = Z E (cid:2) u ( t, x + c ( t, x, e )) − u ( t, x ) (cid:3) η ( e ) λ ( de ) . For ( t, x ) ∈ [0 , T ] × R q , let E xt [ · ] denote the mathematical expectation under the condition that X t = x , i.e., E xt [ · ] := E [ ·| X t = x ]. To relate the FBSDEs in (2.1) with the PIDE in (2.2), weconsider the FBSDEs of the following form X t,xs = x + Z st b ( r, X t,xr ) dr + Z st σ ( r, X t,xr ) dW r + Z st Z E c ( r, X t,xr − , e )˜ µ ( de, dr ) ,Y t,xs = ξ + Z Ts f ( r, X t,xr , Y t,xr , Z t,xr , Γ t,xr ) dr − Z Ts Z t,xr dW r − Z Ts Z E U t,xr ( e )˜ µ ( de, dr ) , (2.4)where the solution is ( X t,xs , Y t,xs , Z t,xs , U t,xs ) and Γ t,xs = R E U t,xs ( e ) η ( e ) λ ( de ) for t ≤ s ≤ T . Notethat the superscripts in (2.4) indicate the fact that the forward SDE in (2.4) starts from thetime-space point ( t, x ) ∈ [0 , T ] × R q .According to Theorem 3.4 in [1], if the terminal condition ξ of the FBSDEs is a functionof X t,xT , defined by ξ = ϕ ( X t,xT ) ( ϕ ( · ) is the terminal condition of the PIDE), then the triple( Y t,xs , Z t,xs , U t,xs ) for t ≤ s ≤ T can be represented by the unique viscosity solution u ( t, x ) ofthe PIDE (2.2) as follows: Y t,xs = u ( s, X t,xs ) ,Z t,xs = σ ( s, X t,xs ) ∇ u ( s, X t,xs ) ,U t,xs = u ( s, X t,xs − + c ( s − , X t,xs − , e )) − u ( s, X t,xs − ) , (2.5)where ∇ u denotes the gradient of u with respect to x and the function Γ t,xs is defined byΓ t,xs = B [ u ]( s, X t,xs ). Particularly, when s = t , we have u ( t, x ) = Y t,xt = E [ Y t | X t = x ].
3. The semi-discrete scheme for FBSDEs with jumps
In this section, we propose a numerical scheme for discrete-time approximation of the FB-SDEs under consideration. Instead of the FBSDEs (2.1), we will use the conditional repre-sentation of the FBSDEs given in (2.4) throughout this section. Specifically, discretizations ofthe forward SDE and backward SDE are discussed in § § § , T ]: T := { t < t < · · · < t N = T } (3.1)with ∆ t n := t n +1 − t n and ∆ t := max ≤ n ≤ N − ∆ t n . We assume that the time partition T has thefollowing regularity: max ≤ n ≤ N − ∆ t n min ≤ n ≤ N − ∆ t n ≤ c , (3.2)where c ≥ T is not a jump-adapted partition. Due to the decoupling of the FBSDEs in (2.4), the forward SDE can be discretized separately.Here we briefly recall some classic numerical schemes and their properties discussed in [25]. Anyof these schemes can serve as the approximation of the forward SDE in our schemes for theFBSDE. By setting t = t n , s = t n +1 and x = X n in (2.4), the forward SDE can be written as X t n ,X n t n +1 = X n + Z t n +1 t n b ( s, X t n ,X n s ) ds + Z t n +1 t n σ ( s, X t n ,X n s ) dW s + Z t n +1 t n Z E c ( s, X t n ,X n s − , e )˜ µ ( de, ds ) , (3.3)where we assume that the solution X t n ,X n s starts at the time instant t = t n and spatial location X t n ,X n t n = X n . By using the Itˆo-Taylor expansion, numerical schemes of strong-order β (or theweak-order β ) [25] can be represented in a general form, i.e., X n +1 = X n + Φ( t n , t n +1 , X n , I J ∈A β ) , (3.4)where Φ is the incremental, A β is a hierarchical set such that the convergence rate of the schemeis β in a strong or weak sense. Details of the index set I J ∈A β and the definition of A β can befound in [25] (pp. 196 and pp. 290). The scheme (3.4) has the following properties: • Stability : for an integer r >
0, there exists a constant C ∈ (0 , ∞ ) such thatmax ≤ n ≤ N E [ | X n | r ] ≤ C (1 + E [ | X | r ]) . (3.5) • Approximation error: there exist positive real numbers r , r , r , α, β, γ such that for anyfunction g ∈ C β +2 P , we have (cid:12)(cid:12)(cid:12) E X n t n [ g ( X t n ,X n t n +1 ) − g ( X n +1 )] (cid:12)(cid:12)(cid:12) ≤ C (1 + | X n | r )(∆ t ) β +1 , (cid:12)(cid:12)(cid:12) E X n t n [( g ( X t n ,X n t n +1 ) − g ( X n +1 ))∆ ˜ W ⊤ t n +1 ] (cid:12)(cid:12)(cid:12) ≤ C (1 + | X n | r )(∆ t ) γ +1 , (cid:12)(cid:12)(cid:12) E X n t n [( g ( X t n ,X n t n +1 ) − g ( X n +1 ))∆˜ µ ∗ t n +1 ] (cid:12)(cid:12)(cid:12) ≤ C (1 + | X n | r )(∆ t ) α +1 , (3.6)where C β +2 P is the set of 2 β +2 times continuously differentiable functions which, togetherwith their derivatives of order up to 2 β + 2, have at most polynomial growth. Accordingto Theorem 6.4.1 and Theorem 12.3.4 in [25], it is easy to derive that α = β = γ forstrong and weak Taylor schemes. In this paper, we prove in Theorem 4.2 that the second-order convergence of the proposed semi-discrete scheme for the FBSDE in (2.1) requires α = β = γ = 2. Now we study the discretization of the backward SDE in (2.4) driven by the process X t n ,X n s in (3.3) for s ∈ [ t n , t n +1 ]. Within the interval [ t n , t n +1 ], the backward SDE can be rewritten as Y t n ,X n t n = Y t n ,X n t n +1 + Z t n +1 t n f t n ,X n s ds − Z t n +1 t n Z t n ,X n s dW s − Z t n +1 t n Z E U t n ,X n s ( e )˜ µ ( de, ds ) , (3.7)where f t n ,X n s denotes f ( s, X t n ,X n s , Y t n ,X n s , Z t n ,X n s , Γ t n ,X n s ) for notational simplicity. Due to therelation between Γ t,xs and U t,xs , in what follows, all the numerical schemes for the backwardSDE will be proposed to approximate ( Y t n ,X n s , Z t n ,X n s , Γ t n ,X n s ). Since there are three unknownstochastic processes involved in (3.7), we now construct three discretized reference equations for Y t n ,X n t n , Z t n ,X n t n and Γ t n ,X n t n in § § § § Y t n ,X n t n Taking the conditional mathematical expectation E X n t n [ · ] on both sides of (3.7), we obtain Y t n ,X n t n = E X n t n h Y t n ,X n t n +1 i + Z t n +1 t n E X n t n h f t n ,X n s i ds, (3.8)due to the fact that R t n +1 t n Z t n ,X n s dW s and R t n +1 t n R E U t n ,X n s ( e )˜ µ ( de, ds ) for t > t n are martingales.Note that the integrand E X n t n [ f t n ,X n s ] is a deterministic function of s ∈ [ t n , t n +1 ] under the σ -algebra F t n . Thus, numerical integration approaches can be used to approximate the temporalintegral in (3.8). In this effort, we use the Crank-Nicolson scheme, i.e., the trapezoidal rule,such that Z t n +1 t n E X n t n h f t n ,X n s i ds = 12 ∆ t n f t n ,X n t n + 12 ∆ t n E X n t n h f t n ,X n t n +1 i + R ny , (3.9)where the residual R ny is R ny := Z t n +1 t n (cid:26) E X n t n h f t n ,X n s i − f t n ,X n t n − E X n t n h f t n ,X n t n +1 i (cid:27) ds. (3.10)Substituting (3.9) into (3.8), we obtain the reference equation for solving Y t n ,X n t n : Y t n ,X n t n = E X n t n h Y t n ,X n t n +1 i + 12 ∆ t n f t n ,X n t n + 12 ∆ t n E X n t n h f t n ,X n t n +1 i + R ny . (3.11) Z t n ,X n t n To proceed, we introduce a new Gaussian process ∆ ˜ W s defined by∆ ˜ W s = 2∆ W s − t n Z st n ( r − t n ) dW r , ∀ s ∈ [ t n , t n +1 ] , (3.12)where ∆ W s = W s − W t n is the d -dimensional standard Brownian motion in the FBSDEs in(2.4). It is easy to see that ∆ ˜ W s = (∆ ˜ W s , ∆ ˜ W s , · · · , ∆ ˜ W ds ) ⊤ is also a d -dimensional Gaussianprocess with the properties E X n t n [∆ ˜ W s ] = 0, E X n t n [∆ ˜ W is ∆ ˜ W js ] = 0 for i = j , and E X n t n [(∆ ˜ W is ) ] = E X n t n h (2∆ W is − t n Z st n ( r − t n ) dW ir ) i = 4( s − t n ) − s − t n ) ∆ t n + 3( s − t n ) ∆ t n , for i = 1 , . . . , d. In the case of s = t n +1 , we have E X n t n [∆ ˜ W it n +1 ] = 0 and E X n t n [(∆ ˜ W it n +1 ) ] = ∆ t n for i = 1 , . . . , d .Multiplying (3.7) by the transpose of ∆ ˜ W t n +1 in (3.12), and taking the conditional mathe-matical expectation E X n t n [ · ] on both sides, we obtain0 = E X n t n h Y t n ,X n t n +1 ∆ ˜ W ⊤ t n +1 i + Z t n +1 t n E X n t n h f t n ,X n s ∆ ˜ W ⊤ t n +1 i ds − E X n t n (cid:20)Z t n +1 t n Z t n ,X n s dW s · ∆ ˜ W ⊤ t n +1 (cid:21) . (3.13)Then, the right endpoint rule is used to discretize the first temporal integral in (3.13), suchthat Z t n +1 t n E X n t n h f t n ,X n s ∆ ˜ W ⊤ t n +1 i ds = ∆ t n E X n t n h f t n ,X n t n +1 ∆ ˜ W ⊤ t n +1 i + R nz, , (3.14)where R nz, := R t n +1 t n E X n t n [ f t n ,X n s ∆ ˜ W ⊤ t n +1 ] ds − ∆ t n E X n t n [ f t n ,X n t n +1 ∆ ˜ W ⊤ t n +1 ] is the residual. For thesecond temporal integral in (3.13), based on the properties of ∆ ˜ W ⊤ t n +1 , we discretize it by − E X n t n (cid:20)Z t n +1 t n Z t n ,X n s dW s · ∆ ˜ W ⊤ t n +1 (cid:21) = −
12 ∆ t n Z t n ,X n t n + R nz, , (3.15)where the residual is R nz, := ∆ t n Z t n ,X n t n − E X n t n [ R t n +1 t n Z t n ,X n s dW s · ∆ ˜ W ⊤ t n +1 ] . Substituting (3.14)and (3.15) into (3.13), we obtain the reference equation for Z t n ,X n t n , i.e.,12 ∆ t n Z t n ,X n t n = E X n t n h Y t n ,X n t n +1 ∆ ˜ W ⊤ t n +1 i + ∆ t n E X n t n h f t n ,X n t n +1 ∆ ˜ W ⊤ t n +1 i + R nz , (3.16)where R nz := R nz, + R nz, . Γ t n ,X n t n Similar to the definition of ∆ ˜ W s , by using the compensated Poisson random measure ˜ µ ( de, ds )in (2.4), we define a new stochastic process ∆˜ µ ∗ s as∆˜ µ ∗ s = Z st n Z E (cid:18) − t − t n )∆ t n (cid:19) η ( e )˜ µ ( de, dt ) , ∀ s ∈ [ t n , t n +1 ] . (3.17)Then, multiplying (3.7) by ∆˜ µ ∗ t n +1 and taking the conditional mathematical expectation E X n t n [ · ]on both sides, we obtain0 = E X n t n h Y t n ,X n t n +1 ∆˜ µ ∗ t n +1 i + Z t n +1 t n E X n t n h f t n ,X n s ∆˜ µ ∗ t n +1 i ds − E X n t n (cid:20)Z t n +1 t n Z E U t n ,X n s ( e )˜ µ ( de, ds )∆˜ µ ∗ t n +1 (cid:21) . (3.18)Analogous to the reference equation (3.16), we also discretize the first temporal integral in(3.18) using the right endpoint rule, such that Z t n +1 t n E X n t n h f t n ,X n s ∆˜ µ ∗ t n +1 i ds = ∆ t n E X n t n h f t n ,X n t n +1 ∆˜ µ ∗ t n +1 i + R n Γ , , (3.19)where R n Γ , := R t n +1 t n E X n t n [ f t n ,X n s ∆˜ µ ∗ t n +1 ] ds − ∆ t n E X n t n [ f t n ,X n t n +1 ∆˜ µ ∗ t n +1 ] is the residual. For thesecond temporal integral in (3.18), we have − E X n t n (cid:20)Z t n +1 t n Z E U t n ,X n s ( e )˜ µ ( de, ds )∆˜ µ ∗ t n +1 (cid:21) = − E X n t n (cid:20)Z t n +1 t n Z E U t n ,X n t n ( e )˜ µ ( de, ds )∆˜ µ ∗ t n +1 (cid:21) + R n Γ , = − E X n t n (cid:20) Γ t n ,X n t n Z t n +1 t n (cid:18) − s − t n )∆ t n (cid:19) ds (cid:21) + R n Γ , = −
12 ∆ t n Γ t n ,X n t n + R n Γ , , (3.20)where R n Γ , := ∆ t n Γ t n ,X n t n − E X n t n [ R t n +1 t n R E U t n ,X n s ( e )˜ µ ( de, ds )∆˜ µ ∗ t n +1 ] . By (3.18), (3.19) and(3.20), we obtain the reference equation for Γ t n ,X n t n , i.e.,12 ∆ t n Γ t n ,X n t n = E X n t n h Y t n ,X n t n +1 ∆˜ µ ∗ t n +1 i + ∆ t n E X n t n h f t n ,X n t n +1 ∆˜ µ ∗ t n +1 i + R n Γ , (3.21)where R n Γ := R n Γ , + R n Γ , . Now we combine the approximation in (3.4) of the forward SDE and the three referenceequations in (3.11), (3.16) and (3.21) to propose our semi-discrete scheme (temporal discretiza-tion) for the FBSDEs in (2.1). Let ( X n +1 , Y n , Z n , Γ n ) denote the approximation to the exactsolution ( X t n ,X n t n +1 Y t n ,X n t n , Z t n ,X n t n , Γ t n ,X n t n ) of the FBSDEs in (2.1) for n = N − , . . . ,
0. Basedon the partition T of the time interval [0 , T ], the approximate solution ( X n +1 , Y n , Z n , Γ n ) isconstructed following the procedure in Scheme 1. Scheme 1.
Given the initial condition X for the forward SDE in (2.1) and the terminal condi-tion ( Y N , Z N , Γ N ) for the backward SDE in (2.1) , solve the approximate solution ( X n +1 , Y n , Z n , Γ n ) ,for n = N − , . . . , , by X n +1 = X n + Φ( t n , t n +1 , X n , I J ∈A β ) , (3.22)12 ∆ t n Z n = E X n t n h Y n +1 ∆ ˜ W ⊤ t n +1 i + ∆ t n E X n t n h f n +1 ∆ ˜ W ⊤ t n +1 i , (3.23)12 ∆ t n Γ n = E X n t n h Y n +1 ∆˜ µ ∗ t n +1 i + ∆ t n E X n t n h f n +1 ∆˜ µ ∗ t n +1 i , (3.24) Y n = E X n t n (cid:2) Y n +1 (cid:3) + 12 ∆ t n f n + 12 ∆ t n E X n t n (cid:2) f n +1 (cid:3) , (3.25) where f n +1 = f ( t n +1 , X n +1 , Y n +1 , Z n +1 , Γ n +1 ) , f n = ( t n , X n , Y n , Z n , Γ n ) , ∆ ˜ W t n +1 and ∆˜ µ ∗ t n +1 are defined according to (3.12) and (3.17) by setting s = t n +1 , respectively. From the dependence of (3.22)–(3.25), we can see that the scheme in (3.22) for X n +1 isindependent of the other three schemes, so that, at each time step, X n +1 is always firstlydetermined. Then, by observing that (3.23) and (3.24) are explicit schemes, we can solve Z n and Γ n by substituting X n +1 into (3.23) and (3.24), respectively. Next, since (3.25) includes f n that depends on Y n , Z n and Γ n , it is an implicit scheme for Y n . If the generator f n is nonlinear and Lipschitz continuous with respect to Y n , then Y n can be obtained by substituting X n +1 , Z n and Γ n into (3.25) and solving a nonlinear equation.In addition, we would like to discuss the application of Scheme 1 to discrete time approxi-mations of the PIDE in (2.2) when the terminal condition ξ of the FBSDEs is a function of X T ,i.e., ξ = ϕ ( X T ). The goal is to construct an approximate solution u ( t n , x ) for n = N − , . . . , , x ∈ R q . Specifically, based on the relation between u ( t, x ) and ( X t , Y t , Z t , Γ t ) in (2.5), thediscrete time approximation, denoted by u n ( x ), is defined by u n ( x ) := E [ Y n | X n = x ] ≈ u ( t n , x ) = E [ Y t n | X t n = x ] for x ∈ R q . (3.26)It is easy to see that both Y t n and Y n are deterministic values under the conditions X t n = x and X n = x , respectively. Moreover, the convergence of ( X n +1 , Y n , Z n , Γ n ) to ( X t n +1 , Y t n , Z t n , Γ t n )0as N → ∞ will ensure the convergence of u n ( x ) to u ( t n , x ). Hence, Scheme 1 can be viewedas an effective probabilistic scheme for the PIDE in (2.2). Moreover, Z n and Γ n provideapproximations of σ ∇ u and B [ u ] which enable accurate characterization of local and nonlocaldiffusive fluxes in practical engineering problems.
4. Error estimates for the semi-discrete scheme
In this section, we estimate the truncation error of Scheme 1. Since error estimates forthe scheme for the forward SDEs have been well established in the literature (see [25] andthe references therein), we focus on analyzing the approximation error of ( Y n , Z n , Γ n ) for n = 0 , . . . , N −
1. The general procedure of our analysis is similar to that for classic time-steppingschemes. We first construct an upper bound of the global truncation error of ( Y n , Z n , Γ n ) byrecursively accumulating local truncation errors. Then, we estimate all the local truncationerrors in the upper bound, which relates the global truncation error to the maximum time stepsize ∆ t := max ≤ n ≤ N ∆ t n .To proceed, we need to specify the definition of the approximation error of ( Y i , Z i , Γ i ). For i = 0 , . . . , N , the errors of Y i , Z i , Γ i and f i := f ( t i , X i , Y i , Z i , Γ i ) are respectively defined by e iy := Y t i ,X i t i − Y i , e iz := Z t i ,X i t i − Z i ,e i Γ := Γ t i ,X i t i − Γ i , e if := f t i ,X i t i − f i , (4.1)where Y t i ,X i t i := E [ Y t i | X t i = X i ] and likewise for Z t i ,X i t i , Γ t i ,X i t i , f t i ,X i t i . It should be noted that Y t i ,X i t i +1 and Y t i +1 ,X i +1 t i +1 for 0 ≤ i ≤ N are usually different stochastic processes because of thedifference between X t i ,X i t i +1 and X i +1 . This fact can be easily shown with the use of the solution u ( t, x ) of the PIDE in (2.2). According to the relationship in (2.5), it is easy to see that Y t i ,X i t i +1 = u ( t i +1 , X t i ,X i t i +1 ) , Y t i +1 ,X i +1 t i +1 = u ( t i +1 , X i +1 ) , where X t i ,X i t i +1 and X i +1 are obtained by (3.3) and (3.22), respectively. As such, we introducethe following residual notations that will be used later: R iy := E X i t i h Y t i ,X i t i +1 − Y t i +1 ,X i +1 t i +1 i , R iy := E X i t i h f t i ,X i t i +1 − f t i +1 ,X i +1 t i +1 i , R iz := E X i t i h(cid:16) Y t i ,X i t i +1 − Y t i +1 ,X i +1 t i +1 (cid:17) ∆ ˜ W ⊤ t i +1 i , R iz := E X i t i h(cid:16) f t i ,X i t i +1 − f t i +1 ,X i +1 t i +1 (cid:17) ∆ ˜ W ⊤ t i +1 i , R i Γ := E X i t i h(cid:16) Y t i ,X i t i +1 − Y t i +1 ,X i +1 t i +1 (cid:17) ∆˜ µ ∗ t i +1 i , R i Γ := E X i t i h(cid:16) f t i ,X i t i +1 − f t i +1 ,X i +1 t i +1 (cid:17) ∆˜ µ ∗ t i +1 i (4.2)for i = 0 , . . . , N −
1. Note that the above residuals represent the local weak approximations ofthe scheme (3.22) for solving forward SDEs.In the following theorem, we construct an upper bound of the errors e ny , e nz and e n Γ (0 ≤ n ≤ N −
1) with the use of the residuals R iy , R iz and R i Γ in (3.11), (3.16) and (3.21), respectively,1for i = n, . . . , N , as well as the residuals defined in (4.2). Theorem 4.1.
Based on the partition T in (3.1) of the time interval [0 , T ] , if the generator f ( t, x, y, z, γ ) is Lipschitz continuous with respect to x , y , z and γ where Lipschitz constant isdenoted by L , then with sufficiently small time step ∆ t := max ≤ i ≤ N ∆ t i , the errors e ny , e nz and e n Γ in (4.1) for n = 0 , . . . , N − can be bounded by E [ | e ny | ] + ∆ t N − X i = n (cid:18) C ∆ t − C ∆ t (cid:19) i − n E [ | e iz | + | e i Γ | ] ≤ C ′ (cid:16) E [ | e Ny | ] + ∆ t E [ | e Nz | + | e N Γ | ] (cid:17) + N − X i = n (cid:18) C ∆ t − C ∆ t (cid:19) i − n C − C ∆ t E " t (cid:0) | R iy | + | R iz | + | R i Γ | (cid:1) + 1∆ t (cid:0) |R iy | + |R iz | + |R i Γ | (cid:1) + ∆ t (cid:0) |R iy | + |R iz | + |R i Γ | (cid:1) , (4.3) where C is a positive constant depending on L and c defined in (3.2) , C ′ is a positive constantdepending on c , T and L , the residuals R iy , R iz and R i Γ for i = n, . . . , N are defined in (3.11) , (3.16) and (3.21) , respectively, and R iy , R iy , R iz , R iz , R i Γ , R i Γ are defined in (4.2) .Proof. This proof consists of four steps. Step 1, 2 and 3 are dedicated to the estimating e ny , e nz and e n Γ , respectively, and those estimates are combined together at Step 4 that completesthe proof. • Step 1 : Estimating the error e ny = Y t n ,X n t n − Y n . Subtracting the scheme (3.25) from the reference equation (3.11), we have e ny = E X n t n h Y t n ,X n t n +1 − Y n +1 i + 12 ∆ t n (cid:16) f t n ,X n t n − f n (cid:17) + 12 ∆ t n E X n t n h f t n ,X n t n +1 − f n +1 i + R ny = E X n t n h Y t n ,X n t n +1 − Y t n +1 ,X n +1 t n +1 + Y t n +1 ,X n +1 t n +1 − Y n +1 i + ∆ t n (cid:16) f t n ,X n t n − f n (cid:17) + ∆ t n E X n t n h f t n ,X n t n +1 − f t n +1 ,X n +1 t n +1 + f t n +1 ,X n +1 t n +1 − f n +1 i + R ny = E X n t n (cid:2) e n +1 y (cid:3) + ∆ t n e nf + ∆ t n E X n t n [ e n +1 f ] + R ny + ∆ t n R ny + R ny . Then under the conditions of the theorem, we have the estimate (cid:12)(cid:12) e ny (cid:12)(cid:12) ≤ (cid:12)(cid:12)(cid:12) E X n t n [ e n +1 y ] (cid:12)(cid:12)(cid:12) + ∆ t n (cid:12)(cid:12)(cid:12) E X n t n [ e n +1 f ] (cid:12)(cid:12)(cid:12) + ∆ t n | e nf | + |R ny | + ∆ t n |R ny | + | R ny |≤ (cid:12)(cid:12)(cid:12) E X n t n [ e n +1 y ] (cid:12)(cid:12)(cid:12) + ∆ t n L E X n t n (cid:2) | e n +1 y | + | e n +1 z | + | e n +1Γ | (cid:3) + ∆ t n L ( | e ny | + | e nz | + | e n Γ | ) + |R ny | + ∆ t n |R ny | + | R ny | . For the squared error | e ny | , given any positive real number γ and positive integer m , by usingthe inequalities ( a + b ) ≤ (1 + γ ∆ t ) a + (1 + γ ∆ t ) b and ( P mn =1 a n ) ≤ m P mn =1 a n , we have2the following estimate. | e ny | ≤ (1 + γ ∆ t ) (cid:12)(cid:12)(cid:12) E X n t n [ e n +1 y ] (cid:12)(cid:12)(cid:12) + (cid:18) γ ∆ t (cid:19) ( ∆ t n L (cid:0) | e ny | + | e nz | + | e n Γ | (cid:1) + ∆ t n L E X n t n (cid:2) | e n +1 y | + | e n +1 z | + | e n +1Γ | (cid:3) + |R ny | + ∆ t n |R ny | + | R ny | ) ≤ (1 + γ ∆ t ) (cid:12)(cid:12)(cid:12) E X n t n [ e n +1 y ] (cid:12)(cid:12)(cid:12) + ( t L (cid:0) | e ny | + | e nz | + | e n Γ | (cid:1) (4.4)+ 15∆ t L E X n t n (cid:2) | e n +1 y | + | e n +1 z | + | e n +1Γ | (cid:3) + 5 |R ny | + 5∆ t |R ny | + 5 | R ny | ) + 1 γ ( tL | e ny | + | e nz | + | e n Γ | )+ 15∆ tL E X n t n [ | e n +1 y | + | e n +1 z | + | e n +1Γ | ] ) + 1 γ ∆ t ( |R ny | + 5∆ t |R ny | + 5 | R ny | ) . • Step 2 : Estimating the error e nz = Z t n ,X n t n − Z n . Subtracting the scheme (3.23) from the reference equation (3.16), we obtain∆ t n e nz = E X n t n h ( Y t n ,X n t n +1 − Y n +1 )∆ ˜ W ⊤ t n +1 i + ∆ t n E X n t n h ( f t n ,X n t n +1 − f n +1 )∆ ˜ W ⊤ t n +1 i + R nz . (4.5)Substituting the identities E X n t n h ( Y t n ,X n t n +1 − Y n +1 )∆ ˜ W ⊤ t n +1 i = R nz + E X n t n h e n +1 y ∆ ˜ W ⊤ t n +1 i , E X n t n h ( f t n ,X n t n +1 − f n +1 )∆ ˜ W ⊤ t n +1 i = R nz + E X n t n h e n +1 f ∆ ˜ W ⊤ t n +1 i , into (4.5), | e nz | can be estimated by | e nz | = (cid:12)(cid:12)(cid:12)(cid:12) t n E X n t n h e n +1 y ∆ ˜ W ⊤ t n +1 i + 2 E X n t n h e n +1 f ∆ ˜ W ⊤ t n +1 i + 2∆ t n R nz + 2 R nz + 2∆ t n R nz (cid:12)(cid:12)(cid:12)(cid:12) ≤ t n (cid:12)(cid:12)(cid:12) E X n t n h e n +1 y ∆ ˜ W ⊤ t n +1 i(cid:12)(cid:12)(cid:12) + 2 (cid:12)(cid:12)(cid:12) E X n t n h e n +1 f ∆ ˜ W ⊤ t n +1 i(cid:12)(cid:12)(cid:12) + 2∆ t n |R nz | + 2 |R nz | + 2∆ t n | R nz | . (4.6)By H¨older’s inequality and the inequality ( a + b ) ≤ (1 + ε ) a + (1 + ε ) b for any positive real3number ε , and from (4.6), we deduce | e nz | ≤ (1 + ε ) (cid:18) t n (cid:19) (cid:12)(cid:12)(cid:12) E X n t n h e n +1 y ∆ ˜ W ⊤ t n +1 i(cid:12)(cid:12)(cid:12) (4.7)+ (cid:18) ε (cid:19) ( (cid:12)(cid:12)(cid:12) E X n t n h e n +1 f ∆ ˜ W ⊤ t n +1 i(cid:12)(cid:12)(cid:12) + 2 |R nz | ∆ t n + 2 |R nz | + 2 | R nz | ∆ t n ) (4.8) ≤ (1 + ε ) (cid:18) t n (cid:19) (cid:12)(cid:12)(cid:12) E X n t n h e n +1 y ∆ ˜ W ⊤ t n +1 i(cid:12)(cid:12)(cid:12) + 16 (cid:18) ε (cid:19) ( E X n t n h | e n +1 f | i E X n t n h | ∆ ˜ W ⊤ t n +1 | i (4.9)+ (cid:18) t n (cid:19) |R nz | + |R nz | + (cid:18) t n (cid:19) | R nz | ) . (4.10)(4.11)By the equality E X n t n [ | ∆ ˜ W t n +1 | ] = ∆ t n , and the estimates of E X n t n [ | e n +1 f | ] and | E X n t n [ e n +1 y ∆ ˜ W ⊤ t n +1 ] | ,i.e., E X n t n h | e n +1 f | i ≤ E X n t n h(cid:12)(cid:12) L ( | e n +1 y | + | e n +1 z | + | e n +1Γ | ) (cid:12)(cid:12) i ≤ L E X n t n h | e n +1 y | + | e n +1 z | + | e n +1Γ | i , (4.12)and (cid:12)(cid:12)(cid:12) E X n t n h e n +1 y ∆ ˜ W ⊤ t n +1 i(cid:12)(cid:12)(cid:12) = (cid:12)(cid:12)(cid:12) E X n t n h ( e n +1 y − E X n t n [ e n +1 y ])∆ ˜ W ⊤ t n +1 i(cid:12)(cid:12)(cid:12) ≤ E X n t n h | ∆ ˜ W ⊤ t n +1 | i E X n t n h ( e n +1 y − E X n t n [ e n +1 y ]) i = ∆ t n (cid:26) E X n t n (cid:2) | e n +1 y | (cid:3) − (cid:12)(cid:12)(cid:12) E X n t n [ e n +1 y ] (cid:12)(cid:12)(cid:12) (cid:27) , (4.13)into ( ?? ), and dividing both sides of the resulting inequality by (1+ ε ) t , we obtain an estimateof | e nz | , i.e.,∆ t ε ) | e nz | ≤ c (cid:26) E X n t n [ | e n +1 y | ] − (cid:12)(cid:12)(cid:12) E X n t n [ e n +1 y ] (cid:12)(cid:12)(cid:12) (cid:27) + 6 L ε ∆ t E X n t n h | e n +1 y | + | e n +1 z | + | e n +1Γ | i + 2∆ tε (cid:26) |R nz | (∆ t n ) + |R nz | + | R nz | (∆ t n ) (cid:27) . (4.14) • Step 3 : Estimating the error e n Γ = Γ t n ,X n t n − Γ n . Subtracting the scheme (3.24) and the reference equation (3.21), we have12 ∆ t n e n Γ = E X n t n h ( Y t n ,X n t n +1 − Y n +1 )∆˜ µ ∗ t n +1 i + ∆ t n E X n t n h ( f t n ,X n t n +1 − f n +1 )∆˜ µ ∗ t n +1 i + R n Γ . (4.15)Substituting E X n t n h ( Y t n ,X n t n +1 − Y n +1 )∆˜ µ ∗ t n +1 i = E X n t n h ( Y t n ,X n t n +1 − Y t n +1 ,X n +1 t n +1 )∆˜ µ ∗ t n +1 i + E X n t n h e n +1 y ∆˜ µ ∗ t n +1 i = R n Γ + E X n t n h e n +1 y ∆˜ µ ∗ t n +1 i E X n t n h ( f t n ,X n t n +1 − f n +1 )∆˜ µ ∗ t n +1 i = E X n t n h ( f t n ,X n t n +1 − f t n +1 ,X n +1 t n +1 )∆˜ µ ∗ t n +1 i + E X n t n h e n +1 f ∆˜ µ ∗ t n +1 i = R n Γ + E X n t n h e n +1 f ∆˜ µ ∗ t n +1 i into (4.15), we obtain an expression of e n Γ as e n Γ = 2∆ t n E X n t n h e n +1 y ∆˜ µ ∗ t n +1 i + 2 E X n t n h e n +1 f ∆˜ µ ∗ t n +1 i + 2∆ t n R n Γ + 2 R n Γ + 2∆ t n R n Γ , (4.16)and consequently we obtain an upper bound of | e n Γ | , i.e., | e n Γ | ≤ t n (cid:12)(cid:12)(cid:12) E X n t n [ e n +1 y ∆˜ µ ∗ t n +1 ] (cid:12)(cid:12)(cid:12) + 2 (cid:12)(cid:12)(cid:12) E X n t n [ e n +1 f ∆˜ µ ∗ t n +1 ] (cid:12)(cid:12)(cid:12) + 2∆ t n |R n Γ | + 2 |R n Γ | + 2∆ t n | R n Γ | . (4.17)By H¨older’s inequality and the inequality ( a + b ) ≤ (1 + ε ) a + (1 + ε ) b for any positive realnumber ε , we obtain the following inequality from (4.17), i.e., | e n Γ | ≤ (1 + ε ) (cid:18) t n (cid:19) (cid:12)(cid:12)(cid:12) E X n t n h e n +1 y ∆˜ µ ∗ t n +1 i(cid:12)(cid:12)(cid:12) + (cid:18) ε (cid:19) (cid:26) (cid:12)(cid:12)(cid:12) E X n t n h e n +1 f ∆˜ µ ∗ t n +1 i(cid:12)(cid:12)(cid:12) + 2 |R n Γ | ∆ t n + 2 |R n Γ | + 2 | R n Γ | ∆ t n (cid:27) ≤ (1 + ε ) (cid:18) t n (cid:19) (cid:12)(cid:12)(cid:12) E X n t n h e n +1 y ∆˜ µ ∗ t n +1 i(cid:12)(cid:12)(cid:12) + 16 (cid:18) ε (cid:19) ( E X n t n h | e n +1 f | i E X n t n h | ∆˜ µ ∗ t n +1 | i + (cid:18) t n (cid:19) |R n Γ | + |R n Γ | + (cid:16) t n (cid:17) | R n Γ | ) . (4.18)By substituting the identity E X n t n [ | ∆˜ µ ∗ t n +1 | ] = ∆ t n R E η ( e ) λ ( de ), the estimate of E X n t n [ | e n +1 f | ]given in (4.12), and the estimate | E X n t n [ e n +1 y ∆˜ µ ∗ t n +1 ] | into (4.18), and dividing both sides of theresulting inequality by ε ) R E ρ ( e ) λ ( de )∆ t , we deduce∆ t ε ) R E η ( e ) λ ( de ) | e n Γ | ≤ c (cid:26) E X n t n (cid:2) | e n +1 y | (cid:3) − (cid:12)(cid:12)(cid:12) E X n t n (cid:2) e n +1 y (cid:3)(cid:12)(cid:12)(cid:12) (cid:27) + 6 L (∆ t ) ε E X n t n (cid:2) | e n +1 y | + | e n +1 z | + | e n +1Γ | (cid:3) + 2∆ tε R E η ( e ) λ ( de ) "(cid:18) t n (cid:19) |R n Γ | + |R n Γ | + (cid:18) t n (cid:19) | R n Γ | . (4.19)5 • Step 4 : Combining the estimates from Steps 1-3.
Now we add (4.4) multiplied by the constant c , (4.14), and (4.19) together to obtain theinequality c | e ny | + ∆ t ε ) | e nz | + ∆ t ε ) R E η ( e ) λ ( de ) | e n Γ | (4.20) ≤ c (cid:20) (cid:18) γ + 15 L γ + 15 L ∆ t L ∆ tc ε (cid:19) ∆ t (cid:21) E X n t n (cid:2) | e n +1 y | (cid:3) + (cid:20) c γ + (cid:18) c ε (cid:19) ∆ t (cid:21) L ∆ t E X n t n (cid:2) | e n +1 z | (cid:3) + (cid:20) c γ + (cid:18) c ε (cid:19) ∆ t (cid:21) L ∆ t E X n t n (cid:2) | e n +1Γ | (cid:3) + (cid:18) c γ + 15 c ∆ t (cid:19) L ∆ t (cid:0) | e ny | + | e nz | + | e n Γ | (cid:1) + 5 c (cid:18) γ ∆ t (cid:19) (cid:26) |R ny | + 14 ∆ t |R ny | + | R ny | (cid:27) + 2∆ tε (cid:26) t n ) |R nz | + |R nz | + 1(∆ t n ) | R nz | (cid:27) + 2∆ tε R E η ( e ) λ ( de ) (cid:26) t n ) |R n Γ | + |R n Γ | + 1(∆ t n ) | R n Γ | (cid:27) . By taking expectation E [ · ] on both sides of (4.20), we deduce c (1 − C ∆ t ) E (cid:2) | e ny | (cid:3) + C ∆ t E (cid:2) | e nz | (cid:3) + C ∆ t E (cid:2) | e n Γ | (cid:3) (4.21) ≤ c (1 + C ∆ t ) E (cid:2) | e n +1 y | (cid:3) + C ∆ t E (cid:2) | e n +1 z | (cid:3) + C ∆ t E (cid:2) | e n +1Γ | (cid:3) + C ∆ t E (cid:20) |R ny | + 14 ∆ t |R ny | + | R ny | (cid:21) + 2∆ tε E "(cid:18) t n (cid:19) |R nz | + |R nz | + (cid:18) t n (cid:19) | R nz | + 2∆ tε R E η ( e ) λ ( de ) E (cid:20) t n ) |R n Γ | + |R n Γ | + 1(∆ t n ) | R n Γ | (cid:21) , where the constants C , C , C , C , C , C are defined by C = (cid:18) γ + 15∆ t (cid:19) L , C = (cid:18) γ + 15 L γ + 15 L ∆ t L ∆ tc ε (cid:19) ,C = 18(1 + ε ) − (cid:18) c γ + 15 c ∆ t (cid:19) L , C = (cid:20) c γ + (cid:18) c ε (cid:19) ∆ t (cid:21) L ,C = 5 c γ ∆ tγ , C = 18(1 + ε ) R E η ( e ) λ ( de ) − (cid:18) c γ + 15 c ∆ t (cid:19) L . (4.22)Now we set ε = 1, γ large enough and ∆ t sufficiently small, such that if 0 < ∆ t ≤ ∆ t then C ≤ C , C ≤ C , C ≤ C , 1 − C ∆ t >
0, and C − C > C ∗ > C − C > C ∗ > C and C ∗ are two positive constants depending on c and L . Then for 0 < ∆ t ≤ ∆ t , we deduce6from (4.21) c (1 − C ∆ t ) E (cid:2) | e ny | (cid:3) + C ∆ t E (cid:2) | e nz | (cid:3) + C ∆ t E (cid:2) | e n Γ | (cid:3) ≤ c (1 + C ∆ t ) E (cid:2) | e n +1 y | (cid:3) + C ∆ t E (cid:2) | e n +1 z | (cid:3) + C ∆ t E (cid:2) | e n +1Γ | (cid:3) + C E (cid:20) t |R ny | + ∆ t |R ny | + 1∆ t | R ny | (cid:21) + C E (cid:20) t |R nz | + ∆ t |R nz | + 1∆ t | R nz | (cid:21) + C E (cid:20) t |R n Γ | + ∆ t |R n Γ | + 1∆ t | R n Γ | (cid:21) . Dividing both sides of the upper inequality by (1 − C ∆ t ), we easily get c E [ | e ny | ] + C ∆ t E [ | e nz | ] + C ∆ t E [ | e n Γ | ] ≤ C ∆ t − C ∆ t (cid:16) c E [ | e n +1 y | ] + C ∆ t E [ | e n +1 z | ] + C ∆ t E [ | e n +1Γ | (cid:17) + C − C ∆ t E (cid:20) t |R ny | + ∆ t |R ny | + 1∆ t | R ny | (cid:21) + C − C ∆ t E (cid:20) t |R nz | + ∆ t |R nz | + 1∆ t | R nz | (cid:21) + C − C ∆ t E (cid:20) t |R n Γ | + ∆ t |R n Γ | + 1∆ t | R n Γ | (cid:21) . (4.23)From the inequality (4.23), by recursively inserting e iy , i = n + 1 , . . . , N −
1, we deduce c E (cid:2) | e ny | (cid:3) + C ∆ t N − X i = n (cid:18) C ∆ t − C ∆ t (cid:19) i − n E (cid:2) | e iz | (cid:3) + C ∆ t N − X i = n (cid:18) C ∆ t − C ∆ t (cid:19) i − n E (cid:2) | e i Γ | (cid:3) ≤ (cid:18) C ∆ t − C ∆ t (cid:19) N − n c E (cid:2) | e Ny | (cid:3) + C ∆ t N X i = n +1 (cid:18) C ∆ t − C ∆ t (cid:19) i − n E (cid:2) | e iz | (cid:3) + C ∆ t N X i = n +1 (cid:18) C ∆ t − C ∆ t (cid:19) i − n E [ | e i Γ | ]+ N − X i = n (cid:18) C ∆ t − C ∆ t (cid:19) i − n C − C ∆ t E (cid:20) t |R iy | + ∆ t |R iy | + 1∆ t | R iy | (cid:21) + N − X i = n (cid:18) C ∆ t − C ∆ t (cid:19) i − n C − C ∆ t E (cid:20) t |R iz | + ∆ t |R iz | + 1∆ t | R iz | (cid:21) + N − X i = n (cid:18) C ∆ t − C ∆ t (cid:19) i − n C − C ∆ t E (cid:20) t |R i Γ | + ∆ t |R i Γ | + 1∆ t | R i Γ | (cid:21) , which immediately leads to c E (cid:2) | e ny | (cid:3) + C ∗ ∆ t N − X i = n (cid:18) C ∆ t − C ∆ t (cid:19) i − n E (cid:2) | e iz | (cid:3) + C ∗ ∆ t N − X i = n (cid:18) C ∆ t − C ∆ t (cid:19) i − n E (cid:2) | e i Γ | (cid:3) ≤ (cid:18) C ∆ t − C ∆ t (cid:19) N − n c E (cid:2) | e Ny | (cid:3) + C ∆ t (cid:18) C ∆ t − C ∆ t (cid:19) N − n E (cid:2) | e Nz | (cid:3) + C ∆ t (cid:18) C ∆ t − C ∆ t (cid:19) N − n E (cid:2) | e N Γ | (cid:3) + N − X i = n (cid:18) C ∆ t − C ∆ t (cid:19) i − n C − C ∆ t E (cid:20) t |R iy | + ∆ t |R iy | + 1∆ t | R iy | (cid:21) + N − X i = n (cid:18) C ∆ t − C ∆ t (cid:19) i − n C − C ∆ t E (cid:20) t |R iz | + ∆ t |R iz | + 1∆ t | R iz | (cid:21) + N − X i = n (cid:18) C ∆ t − C ∆ t (cid:19) i − n C − C ∆ t E (cid:20) t |R i Γ | + ∆ t |R i Γ | + 1∆ t | R i Γ | (cid:21) . The proof is completed.
Remark 1.
It is worth to note that Theorem 4.1 also implies that Scheme 1 is stable withrespect to the terminal condition, that is, for any ε > , there exists a positive number δ > ,such that, if E [ | Y N − Y N | ] < δ , E [ | Z N − Z N | ] < δ and E [ | Γ N − Γ N | ] < δ , where ( Y N , Z N , Γ N ) and ( Y N , Z N , Γ N ) are two different terminal conditions, then for ≤ n ≤ N − , it holds that E [ | Y n − Y n | ] + ∆ t N X i = n E [ | Z i − Z i | ] + ∆ t N X i = n E [ | Γ i − Γ i | ] < ε. The next task is to estimate all the residual terms in (4.3), and the main technique used isthe Itˆo-Taylor expansion [25]. Under some reasonable regularity conditions on the data b , σ , c , f and ϕ in the FBSDE, we now derive estimates for the local truncation errors R ny , R nz and R n Γ defined in (3.11), (3.16) and (3.21), respectively. To proceed, we need the following standardassumption. Assumption 4.1.
Under the condition that X is A -measurable as well as E [ | X | ] < ∞ , weassume that b , σ and c are jointly L -measurable in ( t, x ) ∈ [0 , T ] × R q , and there exist realconstants L > and K > such that | b ( t, x ) − b ( t, x ′ ) | ≤ L | x − x ′ | , | σ ( t, x ) − σ ( t, x ′ ) | ≤ L | x − x ′ | , Z E | c ( t, x, e ) − c ( t, x ′ , e ) | λ ( de ) ≤ L | x − x ′ | , (4.24) and | b ( t, x ) | ≤ K (1 + | x | ) , | σ ( t, x ) | ≤ K (1 + | x | ) , Z E | c ( t, x, e ) | λ ( de ) ≤ K (1 + | x | ) , (4.25) for all t ∈ [0 , T ] and x, x ′ ∈ R q . Under Assumption 4.1, if E [ | X | m ] < ∞ for some integer m ≥
1, the solution of the forwardSDE in (2.1) has the estimate E X n t n h | X t n ,X n s | m i ≤ (cid:16) E X n t n (cid:2) | X n | m (cid:3)(cid:17) e C ( s − t n ) , (4.26)8where s ∈ [ t n , T ] and C is a positive constant depending only on the constants K , L and m .For the sake of presentation simplicity, in the following lemmas and theorems, we only con-sider the one-dimensional case ( q = d = 1), but our results can be extended to multidimensionalcases without any essential difficulty. To proceed, we define three partial integro-differentialoperators: L v ( t, x ) := ∂v∂t ( t, x ) + b ( t, x ) ∂v∂x ( t, x ) + 12 σ ( t, x ) ∂ v∂x ( t, x )+ Z E (cid:20) v ( t, x + c ( t, x, e )) − v ( t, x ) − ∂v∂x ( t, x ) c ( t, x, e ) (cid:21) λ ( de ) ,L v ( t, x ) := σ ( t, x ) ∂v∂x ( t, x ) ,L − v ( t, x ) := v ( t, x + c ( t, x, e )) − v ( t, x ); (4.27)and introduce the following notation: C ( k ,...,k J ) b ( D × · · · × D J ):= ( φ : J Y j =1 D j → R (cid:12)(cid:12)(cid:12)(cid:12) ∂ α · · · ∂ α J φ∂ α x · · · ∂ α J x J is bounded and continuousfor α j ≤ k j , j = 1 , . . . , J, where ~α := ( α , . . . , α J ) ∈ N J ) , where J ∈ N + , D × · · · × D J ⊂ R J .Now we give the estimates of R ny , R nz, , R nz, , R n Γ , and R n Γ , in the following Lemma 1. Lemma 1.
Under Assumption 4.1, if the data of the FBSDEs in (2.1) satisfy the followingregularity conditions: f ( t, x, y, z, γ ) ∈ C (2 , , , , b ([0 , T ] × R ) , b ( t, x ) ∈ C (2 , b ([0 , T ] × R ) , σ ( t, x ) ∈C (2 , b ([0 , T ] × R ) , ϕ ( x ) ∈ C αb ( R ) with α ∈ (0 , , and c ( t, x, e ) ∈ C (2 , , ∞ ) b ([0 , T ] × R ) , then forsufficiently small ∆ t = max n ∆ t n , we have the estimates E (cid:2) | R ny | (cid:3) ≤ C (cid:0) E (cid:2) | X n | (cid:3)(cid:1) (∆ t ) , E (cid:2) | R nz, | (cid:3) ≤ C (1 + E (cid:2) | X n | (cid:3) )(∆ t n ) , E (cid:2) | R nz, | (cid:3) ≤ C (1 + E (cid:2) | X n | (cid:3) )(∆ t n ) , E (cid:2) | R n Γ , | (cid:3) ≤ C (1 + E (cid:2) | X n | (cid:3) )(∆ t n ) , E (cid:2) | R n Γ , | (cid:3) ≤ C (1 + E (cid:2) | X n | (cid:3) )(∆ t n ) , where R ny , R nz, , R nz, , R n Γ , and R n Γ , are defined in (3.10) , (3.14) , (3.15) , (3.19) and (3.20) ,respectively, and C is a positive constant depending only on T , K and upper bounds of thederivatives of b , σ , c , f and ϕ .Proof. Based on the relation shown in (2.5) between the solution ( Y t , Z t , Γ t ) of the BSDE in(2.4) and the solution u ( t, x ) of the PIDE in (2.2), it is easy to prove that u ( t, x ) ∈ C (2 , b ([0 , T ] × R q ) under the regularity conditions given in [1] on f, b, σ, c and ϕ . Then, for t ≤ s ≤ T , thefunction F = F ( t, x ) defined by F ( s, x ) := f ( s, x, u ( s, x ) , ∇ u ( s, x ) σ ( s, x ) , Γ( s, x )) , (4.28)is in the space C (2 , b ([0 , T ] × R q ). Setting x = X t n ,X n s in (4.28) and applying the Itˆo-Taylor9expansion to F ( s, X t n ,X n s ), we obtain F ( s, X t n ,X n s ) = F ( t n , X n ) + Z st n L F ( r, X t n ,X n r ) dr + Z st n L F ( r, X t n ,X n r ) dW r + Z st n Z E L − F ( r, X t n ,X n r − )˜ µ ( de, dr ) , (4.29)where the operators L , L and L − are defined in (4.27).Taking the conditional expectation E X n t n [ · ] on both sides, we have that E X n t n h Z st n L F ( r, X t n ,X n r ) dW r i = 0 and E X n t n h Z st n Z E L − F ( r, X t n ,X n r − )˜ µ ( de, dr ) i = 0 . Then following the same arguments used in the proof of Lemmas 4.2-4.4 in [36], we obtain theestimates of the lemma. The proof is completed.From Lemma 1 and the definitions of R nz in (3.16) and R n Γ in (3.21), we easily get theestimates of R nz and R n Γ , stated in the following lemma. Lemma 2.
Under the conditions of Lemma 1, for sufficiently small time step size ∆ t =max n ∆ t n , we have that E (cid:2) | R nz | (cid:3) ≤ C (cid:0) E (cid:2) | X n | (cid:3)(cid:1) (∆ t ) for ≤ n ≤ N − , E (cid:2) | R n Γ | (cid:3) ≤ C (cid:0) E (cid:2) | X n | (cid:3)(cid:1) (∆ t ) for ≤ n ≤ N − , where R nz , R n Γ are defined in (3.16) , (3.21) , C is a positive constant depending on T , K andthe upper bounds of the derivatives of b , σ , c , f and ϕ . Now combining Theorem 4.1, Lemma 1 and Lemma 2 as well as the estimates given in (3.5)and (3.6), we obtain the convergence rate of Scheme 1 in the following theorem.
Theorem 4.2.
Under Lemma 1 and Lemma 2, if (3.5) and (3.6) hold for the scheme (3.4) for the forward SDE, then, for sufficiently small time step size ∆ t = max n ∆ t n , the errors e ny , e nz and e n Γ in (4.1) for n = 0 , . . . , N − can be bounded by E [ | e ny | ] + ∆ t N − X i = n (cid:18) C ∆ t − C ∆ t (cid:19) i − n E [ | e iz | + | e i Γ | ] ≤ C ( E [ | e Ny | ] + ∆ t E [ | e Nz | + | e N Γ | ]) + C (cid:16) (∆ t ) α + (∆ t ) β + (∆ t ) γ + (∆ t ) (cid:17) , where α , β , γ are defined in (3.6) , C > depends on c and L , C > depends on c , T and L , C > depends on c , T , L , K , X and the upper bounds of the derivatives of b , σ , c , f and ϕ .Proof. From the definitions of R iy , R iy , R iz , R iz , R i Γ and R i Γ in Theorem 4.1, under the0conditions of the theorem, we easily get the estimates E [ | X i | ] ≤ C (1 + E [ | X | ]) , E [ |R iy | ] ≤ C (1 + E [ | X i | r ])(∆ t ) β +2 ≤ C (1 + E [ | X | r ])(∆ t ) β +2 , E [ |R iy | ] ≤ C (1 + E [ | X i | r ])(∆ t ) β +2 ≤ C (1 + E [ | X | r ])(∆ t ) β +2 , E [ |R iz | ] ≤ C (1 + E [ | X i | r ])(∆ t ) γ +2 ≤ C (1 + E [ | X | r )(∆ t ) γ +2 , E [ |R iz | ] ≤ C (1 + E [ | X i | r ])(∆ t ) γ +2 ≤ C (1 + E [ | X | r ])(∆ t ) γ +2 , E [ |R i Γ | ] ≤ C (1 + E [ | X i | r )(∆ t ) α +2 ≤ C (1 + E [ | X | r )(∆ t ) α +2 , E [ |R i Γ | ] ≤ C (1 + E [ | X i | r ])(∆ t ) α +2 ≤ C (1 + E [ | X | r ])(∆ t ) α +2 (4.30)for i = 0 , , . . . , N −
1. By Lemmas 1, 2 and 2, and inequality (3.5), for 0 ≤ i ≤ N −
1, we have E [ | R iy | ] ≤ C (1 + E [ | X | ])(∆ t ) , E [ | R iz | ] ≤ C (1 + E [ | X | ])(∆ t ) , E [ | R i Γ | ] ≤ C (1 + E [ | X | ])(∆ t ) . (4.31)By (4.30) and (4.31), we deduce N − X i = n (cid:16) C ∆ t − C ∆ t (cid:17) i − n C E [ |R iy | + (∆ t ) |R iy | + | R iy | ]∆ t (1 − C ∆ t ) ≤ C (1 + E [ | X | r ] + E [ | X | ])((∆ t ) β + (∆ t ) ) , (4.32) N − X i = n (cid:16) C ∆ t − C ∆ t (cid:17) i − n C E [ |R iz | + (∆ t ) |R iz | + | R iz | ]∆ t (1 − C ∆ t ) ≤ C (1 + E [ | X | r ] + E [ | X | ])((∆ t ) γ + (∆ t ) ) , (4.33)and N − X i = n (cid:16) C ∆ t − C ∆ t (cid:17) i − n C E [ |R i Γ | + (∆ t ) |R i Γ | + | R i Γ | ]∆ t (1 − C ∆ t ) ≤ C (1 + E [ | X | r ] + E [ | X | ])((∆ t ) α + (∆ t ) ) . (4.34)By applying (4.32), (4.33) and (4.34) to Theorem 4.1, we complete the proof.
5. The fully discrete scheme
In this section, we will develop a fully discrete scheme based on Scheme 1 by assuming thatthe jump process of X t in (2.1) has finite activity . This means the Poisson random measure˜ µ ( de, dt ) can be represented by˜ µ ( de, dt ) = µ ( de, dt ) − λρ ( e ) dedt, (5.1)where 0 < λ < ∞ is the jump intensity and ρ ( e ) de is the probability measure of each jump sizesatisfying R E ρ ( e ) de = 1. For jump processes with infinite activities, i.e., λ = ∞ , substantialefforts are needed to construct new spatial discretization approaches, which is out of scope ofthis paper and will be considered in our future works.To proceed, we first introduce a partition the q -dimensional Euclidean space R q by S =1 S × S × · · · × S q , where S k for k = 1 , . . . , q is a partition of the one-dimensional space R , i.e., S k = n x ki (cid:12)(cid:12)(cid:12) x ki ∈ R , i ∈ Z , x ki < x ki +1 , lim i → + ∞ x ki = + ∞ , lim i →−∞ x ki = −∞ o , where ∆ x k := max i ∈ Z {| x ki − x ki − |} and ∆ x = max ≤ k ≤ q ∆ x k . For each multi-index i =( i , i , . . . , i q ) ∈ Z q , the corresponding grid point in S is denoted by x i = ( x i , . . . , x qi q ).Recalling (2.5) and (3.26), we can see that, if the terminal condition ξ is a function of X T ,( Y t,xt , Z t,xt , Γ t,xt ) can be treated as functions of t and X t,xt = x for 0 ≤ t ≤ T . Analogously, thesemi-discrete solution ( Y n , Z n , Γ n ) can be treated as functions of X n . Thus, in this section, wealso write ( Y n , Z n , Γ n ) as functions of x ∈ R q , i.e., Y n ( x ) := E [ Y n | X n = x ] , Z n ( x ) := E [ Z n | X n = x ] , Γ n ( x ) := E [Γ n | X n = x ] . Our objective is to approximate the exact solution ( Y t n ,x i t n , Z t n ,x i t n , Γ t n ,x i t n ) by constructing ( Y n i , Z n i , Γ n i ),such that Y n i ≈ Y n ( x i ) ≈ Y t n ,x i t n , Z n i ≈ Z n ( x i ) ≈ Z t n ,x i t n , Γ n i ≈ Γ n ( x i ) ≈ Γ t n ,x i t n , for n = 0 , . . . , N − i ∈ Z q .To this end, it is critical to develop effective quadrature rules for approximating the condi-tional mathematical expectations E X n t n [ · ] in (3.23)-(3.25). For instance, at each time-space point( t n , x i ) ∈ T × S , approximating Y n ( x i ) using (3.25) requires quadrature rules for E x i t n [ Y n +1 ]and E x i t n [ f n +1 ]. In what follows, we take E x i t n [ Y n +1 ] as an example to propose our new quadraturerule. Slight modifications are needed for approximating E x i t n [ Y n +1 ∆ ˜ W ⊤ t n +1 ] and E x i t n [ Y n +1 ∆˜ µ ∗ t n +1 ];all the proposed quadrature rules can be directly used to estimate the expectations of f n +1 , f n +1 ∆ ˜ W ⊤ t n +1 and f n +1 ∆˜ µ ∗ t n +1 .It is observed that E x i t n [ Y n +1 ] is defined with respect to the probability measure of theincremental stochastic process ∆ X n +1 := X n +1 − x i = Φ( t n , t n +1 , x i , I J ∈A β ) starting from( t n , x i ), where Φ is determined by the selected scheme for the forward SDE. In this section, forthe sake of simplicity, we choose the forward Euler method for the scheme in (3.22), i.e., X n +1 = x i + b ( t n , x i )∆ t n + σ ( t n , x i )∆ W t n +1 + N tn +1 X k = N tn +1 c ( t n , x i , e k ) , (5.2)where N t for t ∈ [ t n , t n +1 ] is the underlying Poisson process. Since (5.2) only achieves first-order convergence in the weak sense, the overall convergence of Scheme 1 will be of first order.High-order schemes for the forward SDE [25], such as order-2.0 weak Taylor scheme, can alsobe used, but the corresponding quadrature rules for approximating E x i t n [ · ] will be dramaticallydifferent from the case of using (5.2). Since the jump intensity λ in (5.1) is finite, the number ofjumps of X t within ( t n , t n +1 ] follows a compensated Poisson distribution N t n +1 − N t n − λ ∆ t n ,where the size of each jump, i.e., c ( t n , x i , e ), follows the distribution ρ ( e ) de . Next, we observethat ∆ W it n +1 = p ∆ t n ξ i for i = 1 , . . . , d, (5.3)where ξ i follows the standard normal distribution N (0 , ̺ ( ξ i ) theprobability density function of ξ i , and by ̺ d ( ξ ) the joint probability density function of ξ =( ξ , . . . , ξ d ) ⊤ .2 Now, we can write out the expression of E x i t n [ Y n +1 ] as E x i t n (cid:2) Y n +1 (cid:3) = ∞ X m =0 P n N t n +1 − N t n = m o E (cid:20) Y n +1 (cid:16) x i + b ( t n , x i )∆ t n + σ ( t n , x i ) p ∆ t n ξ + m X k =1 c ( t n , x i , e k ) (cid:17)(cid:21) = ∞ X m =0 e − λ ∆ t n ( λ ∆ t n ) m m ! E (cid:20) Y n +1 (cid:16) x i + b ( t n , x i )∆ t n + σ ( t n , x i ) p ∆ t n ξ + m X k =1 c ( t n , x i , e k ) (cid:17)(cid:21) = e − λ ∆ t n Z R d Y n +1 (cid:16) x i + b ( t n , x i )∆ t n + σ ( t n , x i ) p ∆ t n ξ (cid:17) ̺ d ( ξ ) dξ + ∞ X m =1 e − λ ∆ t n ( λ ∆ t n ) m m ! ( Z R d Z E · · · Z E Y n +1 (cid:16) x i + b ( t n , x i )∆ t n + σ ( t n , x i ) p ∆ t n ξ + m X k =1 c ( t n , x i , e k ) (cid:17) ̺ d ( ξ ) m Y k =1 ρ ( e k ) dξ de · · · de m ) , (5.4)where c ( t n , x i , e k ) = c ( t n , x i , e k , . . . , e qk ) for k = 1 , . . . , m is the size of the k -th jump and { c ( t n , x i , e k ) } mk =1 follows the joint distribution Q mk =1 ρ ( e k ).Now we study how to approximate E x i t n [ Y n +1 ] in (5.4). First, we observe that the probabilityof having m jumps within ( t n , t n +1 ] is of order O (( λ ∆ t n ) m ), thus the sum of the infinite sequencein (5.4) can be approximated by the sum of a finite sequence by retaining finite number of jumps.We denote by E x i t n ,M y [ Y n +1 ] the approximation of E x i t n [ Y n +1 ] by retaining the first M y jumpswithin ( t n , t n +1 ]. Then, it is easy to see that the error introduced by the truncation is of order O (( λ ∆ t n ) M y +1 ), so that M y = 2 is necessary to match the local truncation error introducedby the semi-discrete scheme in (3.25). An analogous notation E x i t n ,M f [ f n +1 ] is used to representthe approximation of E x i t n [ f n +1 ] by retaining the first M f jumps, where M f = 1 is sufficient tomatch the local truncation error in (3.25).Next, we also need to approximate a d -dimensional integral with respect to ξ for m = 0,and an m × q + d dimensional integral with respect to ( ξ, e , . . . , e m ) for m = 1 , . . . , M y . Thiscan be accomplished by selecting an appropriate quadrature rule based on the properties of ̺ ( ξ ), ρ ( e ) and the smoothness of Y n +1 ( x ) with respect to x . A straightforward choice is touse Monte Carlo methods by drawing samples from ̺ d ( ξ ) and Q mk =1 ρ ( e k ), but they are overallinefficient because of the slow convergence. When Y n +1 ( x ) is sufficiently smooth with respectto x , an alternative way is to use the tensor product of high-order one-dimensional quadraturerules, e.g., Newton-Cotes rules and Gaussian rules, etc. For example, the integrals with respectto ξ in (5.4) can be approximated using the tensor product of the Gauss-Hermite rule [35]. Forthe integrals with respect to ( e , . . . , e m ), the Gauss-Legendre rule is a good choice when ρ ( e ) iscompactly supported, e.g., e follows a uniform distribution; the Gauss-Laguerre rule is appro-priate when ρ ( e ) is the density of an exponential distribution. Without loss of generality, for m = 0 , . . . , M y , we denote by { w mi , s mi } S m i =1 and { v mj , q mj } Q m j =1 to represent the chosen quadraturerule for estimating the integrals in (5.4) with respect to ξ and ( e , . . . , e m ), respectively, where { w mi } S m i =1 , { v mj } Q m j =1 are quadrature weights and { s mi } S m i =1 , { q mj } Q m j =1 are quadrature points. Notethat q mj for j = 1 , . . . , Q m has m components, denoted by { q mj, , . . . , q mj,m } , which correspondto the quadrature abscissa for ( e , . . . , e m ). Then, the approximation of E x i t n [ Y n +1 ], denoted by b E x i t n ,M y [ Y n +1 ], is represented by3 b E x i t n ,M y (cid:2) Y n +1 (cid:3) := e − λ ∆ t n S X i =1 w i Y n +1 (cid:16) x i + b ( t n , x i )∆ t n + σ ( t n , x i ) p ∆ t n s i (cid:17) + M y X m =1 e − λ ∆ t n ( λ ∆ t n ) m m ! " S m X i =1 Q m X j =1 w mi v mj Y n +1 (cid:16) x i + b ( t n , x i )∆ t n + σ ( t n , x i ) p ∆ t n s mi + m X k =1 c ( t n , x i , q mj,k ) (cid:17) . (5.5)Analogously, the approximation of E x i t n [ f n +1 ], denoted by b E x i t n ,M f [ f n +1 ], can be obtained byreplacing Y n +1 with f n +1 in (5.5). For the conditional expectations E x i t n [ Y n +1 ∆ ˜ W ⊤ t n +1 ] and E x i t n [ f n +1 ∆ ˜ W ⊤ t n +1 ], we observe that each component of ∆ ˜ W t n +1 = (∆ ˜ W t n +1 , . . . , ∆ ˜ W dt n +1 ) de-fined in (3.12) can be represented by∆ ˜ W it n +1 = √ ∆ t n (cid:16) ξ i + √ ξ i (cid:17) for i = 1 , . . . , d, where ξ i , ˜ ξ i are independent random variables following standard normal distribution, and ξ i is the same random variable as in (5.3). As such, another d -dimensional integral withrespect to ˜ ξ = ( ˜ ξ , . . . , ˜ ξ d ) is needed in (5.4) and (5.5) to define b E x i t n ,M y [ Y n +1 ∆ ˜ W ⊤ t n +1 ] and b E x i t n ,M f [ f n +1 ∆ ˜ W ⊤ t n +1 ]. For the conditional expectations E x i t n [ Y n +1 ∆˜ µ ∗ t n +1 ] and E x i t n [ f n +1 ∆˜ µ ∗ t n +1 ],we can see that ∆˜ µ ∗ t n +1 defined in (3.17) can be represented by∆˜ µ ∗ t n +1 = Z t n +1 t n Z E (cid:18) − t − t n )∆ t n (cid:19) η ( e ) (cid:2) µ ( de, dt ) − λ ( de ) dt (cid:3) = N tn +1 X k = N tn +1 (cid:18) − τ k − t n )∆ t n (cid:19) η ( e k ) − λ ∆ t n Z E η ( e ) ρ ( e ) de, (5.6)where N t for t ∈ [ t n , t n +1 ] is the standard Poisson process and τ k for k = N t n +1 , . . . , N t n +1 is thejump time instant of the k -th jump within ( t n , t n +1 ]. Hence, E x i t n [ Y n +1 ∆˜ µ ∗ t n +1 ] involves anotherintegral with respect to τ k compared to E x i t n [ Y n +1 ], which requires an additional quadraturerule in (5.5) to construct b E x i t n ,M y [ Y n +1 ∆˜ µ ∗ t n +1 ].Based on the quadrature rules used in (5.5), we observe that it is highly possible the quadra-ture points do not belong to the spatial grid S . In this case, we follow the same strategy asin [33, 35] to resolve this issue, i.e., constructing piecewise Lagrange interpolating polynomialsbased on S to interpolate the integrands at non-grid quadrature points. Again, taking Y n +1 ( x )as an example, it can be approximated by Y n +1 ( x ) ≈ b Y n +1 ( x ) := p +1 X j =1 · · · p +1 X j q =1 " Y n +1( i j ,...,i jq ) q Y k =1 Y ≤ j ≤ p +1 j = jk x k − x ki j x ki jk − x ki j , where b Y n +1 ( x ) is a p -th order tensor-product Lagrange interpolating polynomial and Y n +1( i j ,...,i jq ) is the approximate solution of Y n +1 ( x ) at the spatial point ( x i j , . . . , x qi jq ). For k = 1 , . . . , q ,the interpolation points { x ki j } p +1 j =1 ⊂ S k are the closest p + 1 neighboring points of x k , such that( x i j , . . . , x qi jq ) for j k = 1 , . . . , p + 1 and k = 1 , . . . , q constitute a local tensor-product sub-grid4around x . In summary, the fully discrete scheme of the FBSDEs in (2.1) is given as follows: Scheme 2.
Given initial condition X for the forward SDE in (2.1) and the terminal condition ϕ ( X T ) for the backward SDE in (2.1) , solve the approximate solution ( Y n i , Z n i , Γ n i ) , for n = N − , · · · , and i ∈ Z q , by X n +1 = x i + Φ( t n , t n +1 , x i , I J ∈A β ) , (5.7) Y n i = b E x i t n ,M y h b Y n +1 i + 12 ∆ t n f n i + 12 ∆ t n b E x i t n ,M f h b f n +1 i , (5.8)12 ∆ t n Z n i = b E x i t n ,M y h b Y n +1 ∆ ˜ W ⊤ t n +1 i + ∆ t n b E x i t n ,M f h b f n +1 ∆ ˜ W ⊤ t n +1 i , (5.9)12 ∆ t n Γ n i = b E x i t n ,M y h b Y n +1 ∆˜ µ ∗ t n +1 i + ∆ t n b E x i t n ,M f h b f n +1 ∆˜ µ ∗ t n +1 i , (5.10) where b f n +1 = f ( t n +1 , X n +1 , b Y n +1 , b Z n +1 , b Γ n +1 ) , f n i = ( t n , x i , Y n i , Z n i , Γ n i ) , ∆ ˜ W ⊤ t n +1 and ∆˜ µ ∗ t n +1 are defined in (3.12) and (3.17) with s = t n +1 , respectively. Similar to the semi-discrete scheme, Scheme 2 can be directly used as a fully discrete schemefor the PIDE in (2.2). The solution u ( t n , x i ) of the PIDE is approximated by Y n i for n =0 , . . . , N − i ∈ Z q . We observe that at each grid point ( t n , x i ), the computation of Y n i only depends on ( X n +1 , Y n +1 , Z n +1 , Γ n +1 ) even though an implicit time-stepping scheme isused. This means { Y n i } i ∈ Z q at each time step can be computed independently , so that thedifficulty of solving linear systems with possibly dense matrices, due to the nonlocality of theintegral operator, is completely avoided. This feature makes it straightforward to developmassively parallel algorithms and incorporate adaptive spatial interpolation methods. Remark 2.
It is noted that the total computational cost of the Scheme 2 is dominated by thecost of approximating E x i t n [ · ] at each grid point ( t n , x i ) ∈ T × S using the formula in (5.5) .For example, when solving a three-dimensional problem q = d = 3 and retaining two L`evyjumps M y = M f = 2 , we are facing a large amount of six-dimensional integration problems.In this case, sparse-grid quadrature rules [6, 16, 30] can be used to alleviate the explosion ofcomputational cost due to curse of dimensionality .
6. Numerical examples
In this section, we report on the results of two one-dimensional numerical examples thatillustrate the accuracy and the effectiveness of Schemes 1 and 2. We take uniform partitionsin both temporal and spatial domains with the time and space step sizes denoted by ∆ t and∆ x , respectively. The time step number N is then given by N = T / ∆ t where T is the terminaltime. For the sake of illustration, we only solve FBSDEs on bounded spatial domains. The goalis to test the convergence rates of time discertizaiton and spatial interpolation with respect to∆ t and ∆ x , respectively. To this end, we always set the number of quadrature points to besufficiently large, so that the error contributed by the use of quadrature rules is too small toaffect the convergence rates of interest.5 We consider the following nonlinear
FBSDEs: dX t = dW t + Z E e ˜ µ ( de, dt ) , − dY t = ( ( Y t − · exp( Y t )2 exp (cid:2) sin( X t + t ) + 2 (cid:3) − Z t Y t sin( X t + t ) + 2 − Γ t ) dt − Z t dW t − Z E U ( e )˜ µ ( de, dt ) , (6.1)where the terminal condition is ϕ ( X T ) = sin( X T + T ) + 2. The L`evy measure is defined by λ ( de ) = λρ ( e ) de := X [ − δ,δ ] ( e ) de with δ > , (6.2)where X [ − δ,δ ] ( e ) is the characteristic function of the interval [ − δ, δ ], so that λ = 2 δ is the jumpintensity and ρ ( e ) = δ X [ − δ,δ ] ( e ) is the density function of a uniform distribution on [ − δ, δ ].The exact solution of the FBSDEs is Y t = sin( X t + t ) + 2 ,Z t = cos( X t + t ) , Γ t = cos( X t + t − δ ) − cos( X t + t + δ ) − δ sin( X t + t ) . (6.3)Accordingly, the PIDE corresponding to (6.1) is ∂u∂t + 12 ∂ u∂x + Z E ( u ( t, x + e ) − u ( t, x )) λ ( de )+ ( u − · exp( u )2 exp (cid:2) sin( x + t ) + 2 (cid:3) − u sin( x + t ) + 2 ∂u∂x − B [ u ] = 0 ,u ( T, x ) = sin( x + T ) + 2 , where B [ u ] = cos( x + t − δ ) − cos( x + t + δ ) − δ sin( x + t ).Since the density function ρ ( e ) is uniform with the support [ − δ, δ ], we use the tensor productof the 8-point Gauss-Legendre rule and the 8-point Gauss-Hermite rule to approximate theintegrals involved in E x i t n [ · ].First, we test the convergence rate with respect to ∆ t where the terminal time is T = 1.To this end, we set ∆ x = 0 .
01 and use piecewise cubic Lagrange interpolation to construct b Y n +1 ( x ) for n = 0 , . . . , N −
1, such that the time discretization error dominates the total error.Setting δ = 1 and ∆ t = 2 − , − , − , − , − , the numerical results are shown in Table 6.1.As expected, the convergence rate with respect to ∆ t depends on the number of jumps retainedin b E x i t n ,M y [ · ] and b E x i t n ,M f [ · ]. For example, when M y = M f = 0, i.e., no jump is included, ourscheme fails to converge. In order to achieve second-order convergence, we must set M y ≥ M f ≥ x by setting δ = 1, T = 1, N = 1024, M y = 3, M f = 2, and ∆ x = 2 − , − , − , − , − . The error is measured in L ∞ norm. InTable 6.2, we can see that the spatial discretization error decays as expected, i.e., second-orderand third-order convergence rates for piecewise linear and piecewise quadratic interpolations,respectively.6 Table 6.1: Errors and convergence rates with respect to ∆ t in Example 1, where T = 1, δ = 1,∆ x = 0 .
01, and piecewise cubic Lagrange interpolation is used. k Y ,x − b Y ( x ) k L ∞ ([0 , ∆ t = 2 − ∆ t = 2 − ∆ t = 2 − ∆ t = 2 − ∆ t = 2 − CR M y = 0, M f = 0 5.178E-1 4.786E-1 4.498E-1 4.290E-1 4.090E-1 0.084 M y = 1, M f = 0 4.155E-2 1.778E-2 7.855E-3 3.631E-3 1.691E-3 1.153 M y = 2, M f = 1 4.539E-3 9.878E-4 2.211E-4 5.065E-5 1.144E-5 2.155 M y = 3, M f = 2 3.414E-3 7.305E-4 1.609E-4 3.646E-5 8.087E-6 2.177 k Z ,x − b Z ( x ) k L ∞ ([0 , ∆ t = 2 − ∆ t = 2 − ∆ t = 2 − ∆ t = 2 − ∆ t = 2 − CR M y = 0, M f = 0 1.591E-0 2.155E-0 2.505E-0 2.237E-0 2.534E-0 -0.140 M y = 1, M f = 0 1.629E-1 9.244E-2 5.028E-2 2.156E-2 1.146E-2 0.976 M y = 2, M f = 1 1.846E-2 5.459E-3 1.536E-3 3.293E-4 8.706E-5 1.951 M y = 3, M f = 2 1.475E-2 4.230E-3 1.161E-3 2.450E-4 6.376E-5 1.982 k Γ ,x − b Γ ( x ) k L ∞ ([0 , ∆ t = 2 − ∆ t = 2 − ∆ t = 2 − ∆ t = 2 − ∆ t = 2 − CR M y = 0, M f = 0 5.165E-1 7.206E-1 6.360E-1 5.444E-1 5.270E-1 0.035 M y = 1, M f = 0 3.519E-1 1.805E-1 9.636E-2 4.416E-2 2.240E-2 0.998 M y = 2, M f = 1 2.151E-2 5.364E-3 1.453E-3 3.365E-4 8.458E-5 1.998 M y = 3, M f = 2 1.549E-2 3.856E-3 1.057E-3 2.467E-4 6.212E-5 1.989 Table 6.2: Errors and convergence rates with respect to ∆ x in Example 1, where T = 1, δ = 1, x ∈ [0 , N = 1024, M y = 3 and M f = 2. Linear interpolation∆ x = 2 − ∆ x = 2 − ∆ x = 2 − ∆ x = 2 − ∆ x = 2 − CR k Y ,x − b Y ( x ) k ∞ k Z ,x − b Z ( x ) k ∞ k Γ ,x − b Γ ( x ) k ∞ x = 2 − ∆ x = 2 − ∆ x = 2 − ∆ x = 2 − ∆ x = 2 − CR k Y ,x − b Y ( x ) k ∞ k Z ,x − b Z ( x ) k ∞ k Γ ,x − b Γ ( x ) k ∞ We consider the following nonlinear FBSDE: dX t = sin(2 X t + t ) dt + (cid:2) cos( X t ) + t + 2 (cid:3) dW t + Z E e ˜ µ ( de, dt ) , − dY t = ( − sin(2 X t + t ) Y t Z t [cos( X t ) + t + 2](sin( t ) + 2) exp( − X t ) − . X t ) + t + 2] Y t − Γ t ) dt − Z t dW t − Z E U ( e )˜ µ ( de, dt ) , (6.4)where the terminal condition is ϕ ( X T ) = [sin( T ) + 2] exp( − X T ). The L`evy measure ˜ µ ( de, dt )is defined as in (6.2). The exact solution of the FBSDEs is Y t = (sin( t ) + 2) exp( − X t ) ,Z t = − (cos( X t ) + t + 2)(sin( t ) + 2) exp( − X t ) , Γ t = (sin( t ) + 2) (cid:2) exp( − X t + δ ) − exp( − X t − δ ) − δ exp( − X t ) (cid:3) . (6.5)7Accordingly, the PIDE corresponding to (6.4) is ∂u∂t + sin(2 x + t ) ∂u∂x + 12 (cid:2) cos( x ) + t + 2 (cid:3) ∂ u∂x + Z E ( u ( t, x + e ) − u ( t, x )) λ ( de ) − sin(2 x + t ) u ( t, x )(sin( t ) + 2) exp( − x ) ∂u∂x −
12 [cos( x ) + t + 2] u ( t, x ) − B [ u ] = 0 ,u ( T, x ) = sin( T + 2) exp( − x ) , where B [ u ] = (sin( t ) + 2)[exp( − x + δ ) − exp( − x − δ ) − δ exp( − x )]. Similar to Example 1, weuse the tensor product of the 8-point Gauss-Legendre rule and the 8-point Gauss-Hermite ruleto approximate the integrals involved in E x i t n [ · ]. In this example, we use forward Euler schemein (5.2), such that the overall convergence rate will be expected to be first order.First, we test the convergence rate with respect to ∆ t where the terminal time is T = 1.To this end, we set ∆ x = 0 .
01 and use piecewise cubic Lagrange interpolation to construct b Y n +1 ( x ) for n = 0 , . . . , N −
1, so that the time discretization error dominates the total error.Setting δ = 1 and ∆ t = 2 − , − , − , − , − , the numerical results are shown in Table 6.3.As expected, the convergence rate with respect to ∆ t depends on the number of jumps retainedin constructing b E x i t n ,M y [ · ] and b E x i t n ,M f [ · ]. In this case, we can only achieve, at most, first-orderconvergence with respect to ∆ t due to the use of the forward Euler scheme.Next, we test the convergence rate with respect to ∆ x by setting δ = 1, T = 1, N =1024, M y = 2 and M f = 1. The spatial mesh size is set to ∆ x = 2 − , − , − , − , − forlinear interpolation and ∆ x = 2 − , − , − , − , − for quadratic interpolation. The error ismeasured in L ∞ norm. In Table 6.4, we can see that the spatial discretization error decays asexpected, i.e., second-order and third-order convergence rates for piecewise linear and piecewisequadratic interpolations, respectively. Table 6.3: Errors and convergence rates with respect to ∆ t in Example 2, where T = 1, δ = 1,∆ x = 0 .
01, and piecewise cubic Lagrange interpolation are used. k Y ,x − b Y ( x ) k L ∞ ([0 , ∆ t = 2 − ∆ t = 2 − ∆ t = 2 − ∆ t = 2 − ∆ t = 2 − CR M y = 0, M f = 0 1.331E-1 9.191E-2 7.689E-2 6.763E-2 6.206E-2 0.264 M y = 1, M f = 0 2.996E-2 1.362E-2 5.339E-3 2.226E-3 1.998E-3 1.043 M y = 2, M f = 1 3.071E-2 1.047E-2 3.835E-3 1.613E-3 7.155E-4 1.355 k Z ,x − b Z ( x ) k L ∞ ([0 , ∆ t = 2 − ∆ t = 2 − ∆ t = 2 − ∆ t = 2 − ∆ t = 2 − CR M y = 0, M f = 0 4.757E-1 4.522E-1 5.558E-1 6.583E-1 6.881E-1 -0.161 M y = 1, M f = 0 3.779E-1 1.733E-1 7.756E-2 3.598E-2 1.661E-2 1.128 M y = 2, M f = 1 1.129E-1 5.351E-2 2.714E-2 1.201E-2 5.771E-3 1.074 k Γ ,x − b Γ ( x ) k L ∞ ([0 , ∆ t = 2 − ∆ t = 2 − ∆ t = 2 − ∆ t = 2 − ∆ t = 2 − CR M y = 0, M f = 0 1.379E-1 1.194E-1 8.948E-2 8.391E-2 7.765E-2 0.217 M y = 1, M f = 0 5.125E-2 2.528E-2 1.226E-2 7.598E-3 3.996E-3 0.909 M y = 2, M f = 1 4.789E-2 2.116E-2 1.057E-2 5.373E-3 2.509E-3 1.049
7. Concluding remarks
In this work, we propose new numerical schemes for decoupled forward-backward stochasticdifferential equations with jumps, which feature high-order temporal and spatial convergence8
Table 6.4: Errors and convergence rates with respect to ∆ x in Example 2, where T = 1, δ = 1, x ∈ [0 , N = 1024, M y = 2 and M f = 1. Linear interpolation∆ x = 2 − ∆ x = 2 − ∆ x = 2 − ∆ x = 2 − ∆ x = 2 − CR k Y ,x − b Y ( x ) k ∞ k Z ,x − b Z ( x ) k ∞ k Γ ,x − b Γ ( x ) k ∞ x = 2 − ∆ x = 2 − ∆ x = 2 − ∆ x = 2 − ∆ x = 2 − CR k Y ,x − b Y ( x ) k ∞ k Z ,x − b Z ( x ) k ∞ k Γ ,x − b Γ ( x ) k ∞ rates. This advantage has been verified by both theoretical analysis and numerical experiments.Meanwhile, we also realized that our schemes cannot achieve the desired convergence rates inthe sense that the solution of the FBSDEs does not satisfy the necessary regularity conditions.For example, this may happen in real-world financial problems, such as option pricing. However,the regularity conditions do not limit the applicability of the proposed approach, because ourmethod can be directly employed as a probabilistic scheme for related PIDEs which are widelyused to describe anomalous transport in subsurface flow and plasma physics. In these settings,there is a variety of problems satisfying the regularity conditions, and high-order schemes arehighly desired. Moreover, compared to existing deterministic approaches (e.g., finite elements)for the PIDEs, the ability to completely avoids the solution of dense linear systems, as well asto utilize efficient adaptive approximation, and the potential of massively parallel implementa-tion, make our technique highly advantageous. Our future works will focus on extending theproposed numerical schemes to the case of Poisson random measures with infinite activities,and integrating sparse grid methods for high-dimensional FBSDEs with jumps. Acknowledgments.
This work is partially supported by the National Natural Science Foun-dations of China under grant numbers 91130003 and 11171189; by Natural Science Foundationof Shandong Province under grant number ZR2011AZ002; by the U.S. Department of Energy,Office of Science, Office of Advanced Scientific Computing Research, Applied Mathematicsprogram under contract number ERKJE45; and by the Laboratory Directed Research and De-velopment program at the Oak Ridge National Laboratory, which is operated by UT-Battelle,LLC, for the U.S. Department of Energy under Contract DE-AC05-00OR22725.
References [1] Guy Barles, Rainer Buckdahn, and Etienne Pardoux,
Backward stochastic differential equationsand integral-partial differential equations , Stochastics and Stochastics Reports (1997), no. 1-2,57–83.[2] Dirk Becherer, Bounded solutions to backward SDEs with jumps for utility optimization and in-difference hedging , The Annals of Applied Probability (2006), no. 4, 2027–2054.[3] Christian Bender and Robert Denk, A forward scheme for backward SDEs , Stochastic Processesand their Applications (2007), no. 12, 1793–1812.[4] Bruno Bouchard and Romuald Elie,
Discrete-time approximation of decoupled Forward–BackwardSDE with jumps , Stochastic Processes and their Applications (2008), no. 1, 53–75.[5] Bruno Bouchard, Romuald Elie, and Nizar Touzi,
Discrete-time approximation of BSDEs andprobabilistic schemes for fully nonlinear PDEs , Advanced financial modelling, Walter de Gruyter, Berlin, 2009, pp. 91–124.[6] Hans-Joachim Bungartz and Michael Griebel,
Sparse grids , Acta Numerica (2004), 1–123.[7] Jean-Francois Chassagneux and Dan Crisan, Runge–Kutta schemes for backward stochastic dif-ferential equations , The Annals of Applied Probability (2014), no. 2, 679–720.[8] Chuchu Chen and Jia Hong, Mean-square convergence of numerical approximations for a class ofbackward stochastic differential equations , Discrete and Continuous Dynamical Systems-Series B, (2013), 2051–2067.[9] S Cr´epey, Financial Modeling , A Backward Stochastic Differential Equations Perspective, SpringerFinance, 2013.[10] Dan Crisan and Konstantinos Manolarakis,
Second order discretization of backward SDEs andsimulation with the cubature method , The Annals of Applied Probability (2014), no. 2, 652–678.[11] Lukasz Delong, Backward Stochastic Differential Equations with Jumps and Their Actuarial andFinancial Applications , BSDEs with Jumps, Springer Science & Business, June 2013.[12] Jim Douglas, Jin Ma, and Philip Protter,
Numerical methods for forward-backward stochasticdifferential equations , The Annals of Applied Probability (1996), no. 3, 940–968.[13] Nicole EL Karoui, Shige Peng, and M C Quenez, Backward stochastic differential equations infinance , Mathematical Finance. An International Journal of Mathematics, Statistics and FinancialEconomics (1997), no. 1, 1–71.[14] Anne Eyraud-Loisel, Backward stochastic differential equations with enlarged filtration: Optionhedging of an insider trader in a financial market with jumps , Stochastic Processes and theirApplications (2005), no. 11, 1745–1763.[15] Emmanuel Gobet, Jean-Philippe Lemor, and Xavier Warin,
A regression-based Monte Carlomethod to solve backward stochastic differential equations , The Annals of Applied Probability (2005), no. 3, 2172–2202.[16] Max D Gunzburger, Clayton G Webster, and Guannan Zhang, Stochastic finite element methodsfor partial differential equations with random input data , Acta Numerica (2014), 521–650.[17] Desmond J Higham and Peter E Kloeden, Numerical methods for nonlinear stochastic differentialequations with jumps , Numerische Mathematik (2005), no. 1, 101–119.[18] Antoine Lejay, Ernesto Mordecki, Soledad Torres,
Numerical approximation of Backward Stochas-tic Differential Equations with Jumps , (2007).[19] Ralf Metzler and Joseph Klafter,
The random walk’s guide to anomalous diffusion: a fractionaldynamics approach , Physics Reports (2000), no. 1, 1–77.[20] Bernt Øksendal, Agn`es Sulem, and Agn`es Sulem,
Maximum Principles for Optimal Control ofForward-Backward Stochastic Differential Equations with Jumps , SIAM Journal on control andoptimization (2010), no. 5, 2945–2976.[21] Etienne Pardoux and ShigePeng, Backward stochastic differential equations and quasilinearparabolic partial differential equations , Stochastic Partial Differential Equations and Their Ap-plications, Springer Berlin Heidelberg, Berlin/Heidelberg, January 1992, pp. 200–217.[22] Etienne Pardoux and Shige Peng,
Adapted solution of a backward stochastic differential equation ,Systems & Control Letters (1990), no. 1, 55–61.[23] Shige Peng, Backward stochastic differential equations and applications to optimal control , AppliedMathematics and Optimization (1993), no. 2, 125–144.[24] Shige Peng,, Nonlinear Expectations, Nonlinear Evaluations and Risk Measures , Stochastic Meth-ods in Finance, Springer Berlin Heidelberg, Berlin, Heidelberg, January 2004, pp. 165–253.[25] Eckhard Platen and Nicola Bruti-Liberati,
Numerical Solution of Stochastic Differential Equationswith Jumps in Finance , Stochastic Modelling and Applied Probability, vol. 64, Springer BerlinHeidelberg, Berlin, Heidelberg, 2010.[26] Marie-Claire Quenez and Agn`es Sulem,
BSDEs with jumps, optimization and applications to dy-namic risk measures , Stochastic Processes and their Applications (2013), no. 8, 3328–3357. [27] Situ Rong, On solutions of backward stochastic differential equations with jumps and applications ,Stochastic Processes and their Applications (1997), no. 2, 209–236.[28] Manuela Royer, Backward stochastic differential equations with jumps and related non-linear ex-pectations , Stochastic Processes and their Applications (2006), no. 10, 1358–1376.[29] Shanjian Tang and Xunjing Li,
Necessary Conditions for Optimal Control of Stochastic Systemswith Random Jumps , SIAM Journal on control and optimization (1994), no. 5, 1447–1475.[30] Guannan Zhang, Max Gunzburger, and Weidong Zhao, A sparse-grid method for multi-dimensional backward stochastic differential equations , Journal of Computational Mathematics (2013), no. 3, 221–248.[31] Guannan Zhang, Weidong Zhao, Clayton Webster, and Max Gunzburger, A numerical schemefor nonlocal diffusion equations via backward stochastic differential equations with jumps , ORNLtechnical report 2014/449 (2014).[32] Jianfeng Zhang,
A numerical scheme for BSDEs , The Annals of Applied Probability (2004),no. 1, 459–488.[33] Weidong Zhao, Lifeng Chen, and Shige Peng, A New Kind of Accurate Numerical Method forBackward Stochastic Differential Equations , SIAM Journal on Scientific Computing (2006),no. 4, 1563–1581.[34] Weidong Zhao, Yu Fu, and Tao Zhou, New Kinds of High-Order Multistep Schemes for CoupledForward Backward Stochastic Differential Equations , SIAM Journal on Scientific Computing (2014), no. 4, A1731–A1751.[35] Weidong Zhao, Guannan Zhang, and Lili Ju, A stable multistep scheme for solving backwardstochastic differential equations , SIAM Journal on Numerical Analysis (2010), no. 4, 1369–1394.[36] Weidong Zhao, Wei Zhang, and Lili Ju, A numerical method and its error estimates for thedecoupled forward-backward stochastic differential equations , Commun Comput Phys15