{'fraction_non_alphanumeric': 0.061191358281646756, 'fraction_numerical': 0.026016733963285184, 'mean_word_length': 3.7041315841002547, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Second-order schemes for solving decoupled forward backward stochastic differential equations

In this paper, by using trapezoidal rule and the integration-by-parts formula of Malliavin calculus, we propose three new numerical schemes for solving decoupled forward-backward stochastic differential equations. We theoretically prove that the schemes have second-order convergence rate. To demonstrate the effectiveness and the second-order convergence rate, numerical tests are given.


Introduction
Let (Ω, F , F, P ) be a filtered complete probability space, where F = (F t ) 0 t T is the natural filtration of the standard d-dimensional Brownian motion W = (W t ) 0 t T on (Ω, F , F, P ), and T is a fixed finite horizon. Let L 2 = L 2 F (0, T ) be the set of all F t -adapted and mean-square-integrable vector/matrix processes.
In this paper, on (Ω, F , F, P ), we consider numerical solutions to the decoupled forward-backward stochastic differential equations (FBSDEs) where X 0 is the initial condition of the forward stochastic differential equation (SDE), b is the drift coefficient valued in R d , σ is the diffusion matrix valued in R d×d , ϕ and f valued in R m are the terminal and generator functions of the backward stochastic differential equation (BSDE), respectively, and W t is the standard d-dimensional Brownian motion. Note that the two integrals with respect to W s in (1.1) are the Itô-type integrals.
A triple (X t , Y t , Z t ) : [0, T ] × Ω → R d × R m × R m×d is called an L 2 -adapted solution to (1.1) if it is {F t }-adapted, L 2 -integrable, and satisfies (1.1). In 1990, under some standard conditions on the coefficients of (1.1), Pardoux and Peng [13] originally proved the existence and uniqueness of the solution to nonlinear backward stochastic differential equations (BSDEs). The solution (Y t , Z t ) to (1.1) can be represented as (see [4,14]) (1.2) where u(t, x) is the smooth solution to the following parabolic partial differential equation (PDE): The representations in (1.2) are called the Feynman-Kac formula for the nonlinear PDE. Since the work [13], FBSDEs have been extensively studied and have been found important applications in many fields, such as finance, risk measure, stochastic control, and so on. Thus, methods for solving the solution (X t , Y t , Z t ) to FBSDEs, theoretically or numerically, are needed. It is often difficult to get the solution in closed form. So it is more important to study numerical methods for solving FBSDEs. There are many numerical schemes for solving FBSDEs. Some of them are Euler-type schemes for solving BSDEs or FBSDEs (see [1][2][3]5,10,11,16]); the accuracy of these Euler-type schemes are often low with half-order strong convergence rate. Some of them are second-order (or higher) schemes for solving BSDEs (see [9,17,18,21]).
In this paper, we first extend the second-order scheme proposed for BSDEs in [17] to the decoupled FBSDEs. Then by using the trapezoidal rule and the integration-by-parts formula of Malliavin calculus, we propose two other new numerical schemes for solving decoupled FBSDEs. By Malliavin calculus theory we rigorously prove that the two new schemes have second-order convergence rate. Results of our numerical tests given in this paper show that the proposed schemes are second-order schemes, which is consistent with our theoretical results.
The paper is organized as follows. In Section 2, after introducing some needed preliminaries, we will propose new numerical schemes for solving the decoupled FBSDEs (1.1). In Section 3, we will give general error estimate results, which imply the stability of the schemes. By using the error estimate results obtained in Section 3, we obtain second-order error estimate results of the new schemes for the decoupled FBSDEs (1.1) in Section 4. In Section 5, two numerical examples are given to verify the theoretical results. Some conclusions are given in Section 6. Now we introduce the following notations: • (·) : The transpose operator for a matrix or vector; (·) −1 : the inverse operator for a matrix.
• C l,k,k,k b : The set of continuously differential functions ψ : [0, T ]×R d ×R m ×R m×d → R with uniformly bounded partial derivatives ∂ l1 t ψ and ∂ k1 X ∂ k2 Y ∂ k3 Z ψ for l 1 l and k 1 + k 2 + k 3 k. The notations C l,k b and C l,k,k b are similarly defined. • C k,α b (α ∈ (0, 1)): φ ∈ C k,α b means that it has uniformly bounded partial derivatives up to order k, and ∂ k x φ(x) is Hölder continuous with index α w.r.t. x. • F t,x s (t s T ): σ-field generated by the diffusion process {X r , t r s, X t = x} starting from the time-space point (t, x). When s = T , we use F t,x to denote F t,x T .
The conditional mathematical expectation of the random variable η under the σ-field

Second-order schemes
Suppose that H is a real separable Hilbert space with scalar product denoted by ·, · H . The norm of an element h ∈ H will be denoted by h H . Let B = {B(h), h ∈ H} denote an isonormal Gaussian process associated with the Hilbert space H on (Ω, F , F, P ). Let the operator D k [12] be the Malliavin derivative of order k. For any integer p 1, D k,p is the domain of D k (k ∈ N) in L p (Ω), i.e., D k,p is the closure of the class of smooth random variables F with respect to the norm For the Malliavin derivative operator D, we introduce the following two lemmas.
Lemma 2.1 (See [12]). For F ∈ D 1,2 and u ∈ L 2 (Ω; H), we have and the following Integration-by-parts formula: To propose our numerical schemes, we introduce a regular time partition: 0 = t 0 < · · · < t N = T for the time interval [0, T ]. Let Δt n = t n+1 − t n and Δt = max N −1 i=0 Δt n . We denote the increment W tn+1 − W tn of the Brownian motion W t by ΔW tn+1 .

Numerical schemes for solving SDEs
Since the FBSDEs (1.1) are decoupled, the SDE and the BSDE in (1.1) can be solved separately. There have been many numerical schemes for solving SDEs in literature, such as the Euler scheme, the Milstein scheme and the Itô-Taylor type schemes [7,15].
We assume that the numerical scheme used for solving the SDE in (1.1) is in the following form: where X 0 = X 0 , φ is a given R d -valued function, and ξ n+1 n is a random vector related to ΔW tn+1 = W tn+1 −W tn , and under some conditions on the coefficients b and σ the function φ has the estimates (2.1) and (2.2) in the following two Hypothesis.
and the Malliavin derivative D t X n+1 of X n+1 defined by (2.1) has the approximate property Hypothesis 2.2. The approximation solution X n+1 to (2.1), n = 0, 1, 2, . . . , N −1, have the properties where g ∈ C 2β+2 b , β and γ are nonnegative numbers, C g is a positive number which does not depend on time partition.
Numerical schemes for one-dimensional and multi-dimensional SDEs have been extensively studied [7], and many of them can be written in the form of (2.1), where the φ satisfies the above two hypothesis, such as the Euler scheme, the Milstein scheme and the order-2 weak Itô-Taylor type schemes. In the following, we take the three schemes for one-dimensional SDEs as examples to illustrate Hypothesis 2.1.
1. The Euler scheme is in the form of (2.1) with the φ defined by 2. The Milstein scheme is in the form of (2.1) with the φ defined by The Malliavin derivative 3. The following scheme is an order-2 weak Itô-Taylor type scheme: Obviously, the scheme (2.4) is in the form of (2.1). And the Malliavin derivative Here, for function ψ = ψ(t, x), we use ψ n to denote ψ(t n , X n ).

Remark 2.1.
It is well known in [7] that under proper conditions on b and σ, the following conclusions hold: 1. The Euler scheme is a strong-order 0.5 and weak order 1.0 scheme, and Hypothesis 2.1 and 2.2 hold with β = γ = 1; 2. The Milstein scheme is a strong-order 1.0 and weak order 1.0 scheme, and Hypothesis 2.1 and 2.2 hold with β = γ = 1; 3. The order-2 weak Itô-Taylor scheme is a strong-order 1.0 and weak order 2.0 scheme, and Hypothesis 2.1 and 2.2 hold with β = γ = 2.

Reference equations for solving BSDEs
Let the triple (X tn,X n t , Y tn,X n t , Z tn,X n t ) t∈ [tn,T ] be the solution to the following FBSDEs: ]. (2.9) Now, we multiply BSDE (2.7) by ΔW tn+1 and then take E X n tn [·] on both sides of the derived equation to obtain which by trapezoidal rule again yields Based on the two reference equations (2.8) and (2.10), we will propose our second order schemes in the following subsections.

Remark 2.2.
Note that except for the initial condition for the forward SDEs, (1.1) and (2.6) have the same coefficients b, σ, f and ϕ. Thus, according to the nonlinear Feynman-Kac formula, their solution can be represented as respectively, where u is the unique smooth solution to the PDE (1.3).

Scheme I: Crank-Nicolson-type scheme
Based on (2.8) and (2.10), Zhao et al. [17] proposed a kind of θ-scheme for solving BSDEs. The θscheme becomes a Crank-Nicolson-type Scheme by letting the parameters equal 1 2 in the θ-scheme. Now by directly withdrawing the truncation terms in (2.8) and (2.10), we extend the Crank-Nicolson-type Scheme in [17] to the decoupled FBSDEs (1.1) as follows.

Scheme 2.1.
Given Y N = ϕ(X N ), Z N = ϕ x (X N )σ(t N , X N ) for n = N − 2, . . . , 1, 0, solve random variables X n+1 , Z n and Y n by Different from getting the above Crank-Nicolson-type scheme, in the following two subsections, we will propose two new second-order schemes for solving FBSDEs (1.1) by using the integration-by-parts formula in Lemma 2.1 to the second term on the right-hand side of (2.10).

Scheme II
By Lemma 2.1, we deduce Then using trapezoidal rule to approximate the above integral gives where I d is the d× d identity matrix, σ j (·) is the j-th column of the matrix σ(·), ∂ x σ j is the Jacobi-matrix of σ j (·) w.r.t x, and ∇ x X tn,X n t is the variation of X tn,X n t w.r.t X n . From (2.13) and (2.15), we obtain which by using the Euler method to approximate ∇ x X tn,X n tn+1 defined by (2.16) yields Then by (2.10) and (2.17), we deduce where R n Z , R n Z1 and R n Z2 are defined by (2.11), (2.14) and (2.18), respectively, and Based on (2.1) and the two reference equations (2.8) and (2.19), we propose our second scheme in the following:

Scheme III
In this subsection, without any discretization of E X n tn [Y tn,X n tn+1 ΔW tn+1 ], we rewrite (2.10) in the following equivalent form: Just as what we did to (2.12), by Lemmas 2.1 and 2.2, we get Then by (2.21) and (2.23), we deduce (2.24) From (2.1), (2.21) and (2.24), we propose our third scheme for solving the decoupled FBSDEs (1.1) as follows: . . , 1, 0, solve random variables X n+1 , Z n and Y n by and for the order-2 weak Itô-Taylor type scheme (2.4),

Remark 2.4.
Since the terminal condition for BSDEs is only given to Y T . In Schemes 2.1-2.3, the approximation Z N of Z T is needed. From the theory of PDEs (see [4]), under certain conditions on data b, σ, f and ϕ, it holds that u x (T, x) = ϕ x (x). Then according to (1.2), we have Z T = ϕ x (X T )σ(T, X T ). In this case, we let Z N = ϕ x (X N )σ(T, X N ). Otherwise, if the data are not good enough, we can first obtain Y N −1 and Z N −1 by the Euler scheme and then use these schemes backwardly to get Y n and Z n for n = N − 2, . . . , 1, 0. Since the local truncation error of (2.25) is O(Δt 2 ), finite steps of computations by (2.25) do not influence the accuracy of the three schemes.

Remark 2.5.
Schemes 2.1-2.3 are time discrete schemes. To use these schemes to solve decoupled FBSDEs, the conditional mathematical expectations should be approximated. In this paper, we will obtain the error estimate results for the time discrete Schemes 2.2 and 2.3 in Sections 3 and 4, and use the Guass-Hermite quadrature rule to approximate the conditional expectations in our numerical experiments. The details of the approximations will be introduced in Section 5.

General error estimates
In this section, we only give general error estimate results for Schemes 2.2 and 2.3. It should be noted that we do not know how to derive the error estimate for Scheme 2.1.
We assume that the time partition step Δt n has the regularity constraint: where c 0 is a positive constant. In the sequel, C is a generic constant which does not depend on the time partition and may be different from line to line. In order to get the error estimates of Schemes 2.2 and 2.3, we also make the following Hypothesis.
Then by (2.8) and Scheme 2.2, we obtain Δt n e n f + where L is the Lipschitz constant of f . It follows from taking square on both sides of (3.1) and applying the inequality (a + b) 2 By (2.19) and Scheme 2.2, we deduce By taking square on both sides of the equation (3.3), the inequality (a 1 + a 2 + · · · + a 7 ) 2 7(a 2 1 + a 2 2 + · · · + a 2 7 ) and Hölder's inequality, we get Thanks to Cauchy-Swartz inequality, we obtain Under the conditions of the theorem, we have where L is the Lipschitz constant of σ, we deduce By the inequalities (3.4), (3.5) and (3.7), we obtain Dividing both sides of the inequality (3.8) by 7(Δtn) 2 Δt leads to By choosing γ 1 = 14Cc 0 in (3.2), multiplying (3.2) by c 0 and plugging the derived inequality into (3.9), we deduce which by combining the following inequality: The proof is completed. 0, 1, . . . , N) be the solutions to the FBSDEs (2.6) and Scheme 2.3, respectively. Then under the conditions of Theorem 3.1, for sufficiently small time step Δt, we have . . , N. By (2.24) and Scheme 2.3, we deduce Then under the conditions of the theorem, by the techniques used in the proof of Theorem 3.1, we can get the estimate (3.10) similarly. The proof is completed.

Error estimates
In this section, we first give error estimates for Schemes 2.2 and 2.3.
whereȲ tn andZ tn are the values of Y t and Z t at the time-space point (t n , X n ).
Remark 4.1. Theorem 4.1 implies that the convergence rate is 2.0 for solving Y t , and 1.5 for solving Z t . But our numerical tests in Section 5 show that the convergence rate is 2.0 for solving both Y t and Z t . So the error estimate in this theorem is not optimal for Z t . Up to now, by our analysis we can only get the estimate (4.1). The optimal error estimate for these schemes should be studied in our future work.
To prove the above theorem, we introduce Lemmas 4.1-4.4. and ϕ ∈ C 4,α b (α ∈ (0, 1]), we have the estimates For simple presentation, in Lemmas 4.2 and 4.3, we only consider FBSDEs with m = d = 1. The error estimate results in these two lemmas are true for general multidimensional FBSDEs.

Numerical experiments
In this section, we will perform two numerical experiments to demonstrate the effectiveness and accuracy of the three schemes proposed in this paper. In order to numerically solve the FBSDEs by Schemes 2.1-2.3, the space partition is also needed. We partition the d-dimensional Euclidean space R d by where R k h (k = 1, . . . , d) are the partition of (−∞, ∞), To solve the values of Y n and Z n at every space grid point x i for n = N −1, . . . , 1, 0 backwardly, we need to compute the conditional mathematical expectations By the definitions of the schemes, these random variables are functions of X n+1 and ΔW tn+1 , and then functions of t n , x i and ξ n+1 n , where ξ n+1 n is a Gaussian distributed random vector. By the Gram-Schmidt orthogonalization process, this random vector can be represented as ξ n+1 n = ξ n+1 n (U ) = BU in law, where B is a deterministic matrix and U = (U 1 , U 2 , . . . , U k ) with U i (i = 1, 2, . . . , k) independent and identity distributed standard normal random variables. Thus, in the sense of distribution, G can be represented as a function of t n , x i and ξ n+1 n . We denote this function by H(t n , x i , U), i.e., G = H(t n , x i , U). Now, we turn the approximation of E x i tn [G] to the approximation of E[H(t n , x i , U)], the expectation of functions of standard normal random vector.
It is well known that the Gauss-Hermite Quadrature is a very efficient and accurate method to approximate the expectations of functions of Gaussian type random variable. In our numerical experiments, we choose this method to approximate the conditional expectations in our schemes.
The Gauss-Hermite quadrature rule is an approximation for the values of integrals in the following form by a weighted sum of function values at specified points, where ω i (i = 1, 2, . . . , L) are weights defined by ω i = 2 L+1 L!  R(g, L) is the truncation error expressed as where ζ is some real number. For a k-dimensional function g(x 1 , x 2 , . . . , x k ), the Gauss-Hermite quadrature formula is According to this formula, the conditional mathematical expectations appearing in the proposed schemes can be approximated as ) on time level t n+1 are needed. Generally, X n+1 does not belong to R h . In our numerical experiments, we use the Lagrange interpolation method to approximate them.
By the theory of Lagrange interpolation, when the Lagrange polynomials of order p are used, the errors of the interpolation is O(h p+1 ). Since the aim of our experiments is to check the convergence rate with respect to time, the discrete space step h should be chosen in such a way that the errors resulted from the discretization of FBSDEs and from the space interpolation are balanced. If the theoretical convergence rate of the time discrete FBSDE scheme is α, then h should be selected such that the local errors caused by the FBSDE scheme and the errors caused by the Lagrange interpolation are of the same, i.e., We note that in our numerical experiments, the number of the Gauss-Hermite quadrature points L is also set to be big enough (we take 8 points) such that the errors contributed by the use of Gauss-Hermite quadrature are small and affect the convergence rate very little. The three schemes proposed in this paper have common structures. Besides the terminal conditions, each scheme consists of three parts. The first part is a numerical scheme for solving the forward SDE, the second part is an explicit scheme for Z n (x i ) and the last part is an implicit scheme for Y n (x i ). The detailed solving procedure is given in the following algorithm.
where H stands for the necessary functions of t n , x i and ξ n+1 n whose conditional expectations will be used in the second and third parts in the three schemes, respectively.
(b) Compute the conditional expectations E x i tn [H] by the Gauss-Hermite quadrature rule. (c) Compute Z n (x i ) by the second part of each scheme for each i ∈ Z d , which is explicit for Z n (x i ). (d) Compute Y n (x i ) according to the third part of each scheme by iteration methods, which is implicit for Y n (x i ).
Step 3. If n = 0, then stop, otherwise go to Step 2. Remark 5.1. In many practical problems, such as option pricing problems, we are only interested in the values of (Y 0 (x i ), Z 0 (x i )) for x i in compact domain. Thus, under some accuracy, we just compute (Y n (x i ), Z n (x i )) for x i in a bounded domain instead of all the x i ∈ Z d .
According to Theorem 4.1, the convergence rate of Schemes 2.2 and 2.3 are also restricted by the weak convergence rate of the SDE scheme. Therefore, in our numerical tests, three different schemes for SDEs, the Euler scheme, the Milstein scheme and the order-2.0 weak Itô-Taylor scheme given in Subsection 2.1, are used to solve SDEs.
In the following examples, we take uniform partition in the time interval [0, T ], and set the terminal time T = 1. The time step sizes ranges from T 2 5 to T 2 9 . In the following tables, the notations |Y 0 − Y 0 | and |Z 0 − Z 0 | represent the errors between the exact solution (Y t , Z t ) of (1.1) at time t = 0 and the solution (Y n , Z n ) to the three schemes at n = 0. CR stands for the convergence rate with respect to time step Δt, and is obtained by using linear least square fitting to the errors.

Example 5.1 (Linear case). The considered FBSDEs is
with X 0 = 1.0 and Y T = exp(T ) ln(sin(X T ) + 3.0). The analytical solutions Y t and Z t are Y t = exp(t) ln(sin(X t ) + 3.0) and Z t = exp(t) cos 2 (X t ).
Then the exact solution (Y t , Z t ) at the initial time is (Y 0 , Z 0 ) = (1.044, 0.292). Choose the number of points used in the Lagrange interpolation to be 4. The errors |Y 0 − Y 0 |, |Z 0 − Z 0 |, and the convergence rates for different time partitions and different schemes are listed in Tables 1-3.
The run time of the three schemes for Example 5.1 are listed in Table 4, from which we can see that when the Euler scheme is applied, the run time of the three FBSDE schemes are shorter than the other two SDE schemes. And for the same SDE scheme, Scheme 2.3 is the most efficient.
The terminal condition of Y t is Then the analytic solutions of (5.3) are and Z t = (exp(t + X t )) 2 (1 + exp(t + X t )) 3 , respectively. Assume     In this example, the BSDE is nonlinear, which is different from Example 5.1, in which BSDE is linear. Table 8 shows the run time of Schemes 2.1-2.3 when the three different SDE schemes are applied, respectively.
From the errors and convergence rates listed in Tables 1-7, we conclude that: 1. The convergence rates of the three schemes for solving the decoupled FBSDEs, linear or nonlinear, depend only on the weak order of schemes used for solving the forward SDE. The convergence rates of the three schemes are of order 1 if the Euler scheme or the Milstein scheme is used for solving the SDE, and are of order 2 if the order-2 weak Itô-Taylor scheme is used for solving the SDE.
2. The convergence rates of Schemes 2.2 and 2.3 listed in Tables 2, 3, 6 and 7 are consistent with the theoretical results in Theorem 4.1. Tables 1 and 5 tell us that Scheme 2.1 is also a second-order scheme if order-2 weak schemes for solving SDEs are used.
3. Tables 2, 3, 6 and 7 also show that among the three schemes Schemes 2.1 and 2.3 are more accurate than Scheme 2.2. For Scheme 2.1, we have done many numerical tests, and all the tests showed that it was a second-order scheme for solving decoupled FBSDEs. Unfortunately, for Scheme 2.1, up to now, we still do not know how to theoretically get the error estimate results as stated in Theorem 4.1, though in [20], we have theoretically gotten second-order error estimate results when Scheme 2.1 is used to solve BSDEs.   8 show that the computational cost increased dramatically when we use the Milstein scheme and the weak order-2 scheme to solve the SDEs instead of the Euler scheme. The increasing is mainly caused by the fast increasing of the computation times of space interpolations when h decreases. 5. The three schemes are more efficient when second-order (weak) scheme is used to solve SDEs. For example, in Table 8 of Example 5.2, the errors for solving Y 0 by Scheme 2.3 are similar when the Euler scheme is used with N = 512 and the weak order-2.0 scheme is used with N = 32, but the run times are 43.643 seconds and 0.806 seconds, respectively. The ratio is approximately 54.

Conclusions
We proposed three new high-order numerical schemes for solving the decoupled FBSDEs (1.1) by trapezoidal rule and the integration-by-parts formula of Malliavin calculus. By using Malliavin calculus and stochastic analysis theory, under some reasonable regularity conditions, we rigorously obtained general error estimate results for Schemes 2.2 and 2.3 for solving decoupled FBSDEs. Then we theoretically proved that the accuracy of Schemes 2.2 and 2.3 depends on the accuracy of numerical methods for solving SDEs in weak sense. Numerical experiments showed that the three numerical schemes are first-order ones if the weak-order 1 Euler and Milstein schemes are used for solving SDEs, and are second-order ones if weak-order 2 Itô-Taylor-type schemes are used. These convergence rate results obtained in the experiments are consistent with our theoretical results. Our theoretical and numerical experiment results show that the proposed schemes in this paper are stable, effective and accurate for solving decoupled FBSDEs. These schemes may be used in finance, stochastic control and other related problems. We will do some study on our schemes for these problems in our future research.