On Numerical Approximations of Forward ... - Semantic Scholar

On Numerical Approximations of Forward-Backward Stochastic Differential Equations ∗ Jin Ma

Jie Shen

and

Yanhong Zhao

Department of Mathematics, Purdue University West Lafayette, IN 47907-1395, USA [email protected], [email protected], [email protected]

Abstract. A numerical method for a class of forward-backward stochastic differential equations (FBSDEs) is proposed and analyzed. The method is designed around the Four Step Scheme (Douglas-Ma-Protter, 1996) but with a Hermite-spectral method to approximate the solution to the decoupled quasilinear PDE on the whole space. A rigorous synthetic error analysis is carried out for a fully discretized scheme, namely a first-order scheme in time and a Hermite-spectral scheme in space, of the FBSDEs. Equally important, a systematical numerical comparison is made between several schemes for the resulting decoupled forward SDE, including a stochastic version of the Adams-Bashforth scheme. It is shown that the stochastic version of the Adams-Bashforth scheme coupled with the Hermite-spectral method leads to a convergence rate of 32 (in time) which is better than those in previously published work for the FBSDEs. Keywords: Forward-backward SDEs, four step scheme, spectral method, Hermite functions, stochastic Adam-Bashforth scheme. 2000 Mathematics Subject Classification. Primary: 60H10, Secondary: 34F05, 93E03

1

Introduction

Since the seminal work of Pardoux-Peng [25] in early 1990’s, the theory of backward and forwardbackward stochastic differential equations (BSDEs and FBSDEs, respectively, for short) has grown into an ubiquitous tool in the fields of stochastic optimizations and mathematical finance (see, for example, the books of El Karoui-Mazliak [12], Ma-Yong [21], and the survey of El Karoui-PengQuenez [13] for the detailed account for both theory and applications of such SDEs). In the mean time, after the earlier works of Bally [1] and Douglas-Ma-Protter [11], finding an efficient numerical scheme for both BSDEs and FBSDEs has also become an independent, but integral part of the theory. Tremendous efforts have been made during the past decade to circumvent the fundamental difficulties caused by the combination of the “backward” nature of the SDEs and the associated decoupling techniques for FBSDEs. In the “pure backward” (or “decoupled” forward-backward) case, various methods have been proposed. These include: PDE method in the Markovian case (cf. e.g., Chevance [7], Zhang-Zheng [28]), random walk approximations (cf. e.g., Briand-Delyon-Memin [6] and Ma-Protter-San Martin-Torres [19]), Malliavin calculus and Monte-Carlo method (cf. e.g., Zhang [27], Ma-Zhang [22], and Bouchard-Touzi [5]), and recently the quantization method (cf. e.g., ∗ This work is supported in parts by NSF grants DMS-0311915, DMS-0505427, DMS-0610646 and a Purdue Research Foundation grant.

1

Bally-Pag`es [2] and Bally-Pag´es-Printemps [3]). In the case of coupled FBSDEs, however, the results are very limited, due largely to the lack of solution method itself in such cases. It has been well-understood that in order to solve an FBSDE in an arbitrary duration, a (numerically) tractable method is to utilize the “decoupling PDE”, based on the so-called “Four Step Scheme” initiated in Ma-Protter-Yong [21]. Such an idea has led directly to most of the existing numerical results for FBSDEs: from the early work [11] to the recent improvements by Milstein-Tretyakov [24] and Delarue-Menozzi [10]. Extending such decoupling idea, and combining with an optimal control method as well as Monte-Carlo simulation, Cvitanic-Zhang [8] and Bender-Zhang [4] recently proposed a numerical scheme for FBSDEs without attacking the associated PDE directly. We should note that the Monte-Carlo method is most effective only when a single value of the solution is concerned at each computing cycle, which is quite different from the original problem where essentially the distribution of the solution at each time point was sought. There are two main technical obstacles in developing numerical schemes for FBSDEs: the dimensionality and the rate of convergence. The former is a natural consequence of the close relation between the FBSDEs and its decoupling quasilinear PDE, where the notorious “curse of dimensionality” is still a formidable difficulty for any numerical method. In fact, to our best knowledge, there is still no efficient numerical method for high dimensional PDEs that is directly applicable to our case. We note that almost all the existing numerical schemes for FBSDEs are constructed with α = 1 so their convergence rate is 12 , except in [24], where the Euler scheme for the forward SDE is replaced by Milstein’s first order scheme so the rate of convergence is improved to 1. We should mention that although higher order approximation is possible for forward SDEs (cf. e.g., [18]), it is by no means clear whether this can be extended to the coupled FBSDEs, even when all the coefficients are assumed to be smooth(!), given the intrinsic difficulties arising from the current numerical methods. This then raises an interesting question: is it possible to design a numerical scheme that has better than first-order convergence rate, and at the same time is applicable (in the practical sense) to high dimensional cases? This paper is can be considered as the first step towards this goal. We shall revisit the Four Step Scheme again, but replacing the usual finite difference method for the PDE by a Hermite-spectral method. Several main features of the Hermite-spectral method are worth noting: (i) the PDE on the whole space is approximated directly by using Hermite functions which form an orthonormal basis on the whole space, while the whole space is replaced by a boundary domain with ad-hoc artificial boundary conditions in previous approaches using finite differences; (ii) the rate of convergence is related explicitly to the regularity of the coefficients, to wit, the higher regularity implies the higher convergence rate. We should remark that this last feature marks the main difference between a spectral method and traditional finite difference and finite element methods. In fact, as we pointed out before, even given smooth coefficients, none of the existing schemes for FBSDE seem to be able to achieve higher than first-order convergence rate. Given the aforementioned spectral accuracy in space, it seems to be hopeful that one can produce arbitrarily higher order global error by applying higher order scheme in SDEs in, e.g., [18], at least when the coefficients are smooth. We shall carry out some numerical simulation to validate this point. To compare with the existing results, we shall use an example proposed in Milstein-Tretyakov [24], and we test several methods for the forward SDEs. These include: the standard Euler scheme, the Milstein scheme, the Platen-Wagner scheme, and the Stochastic Adam-Bashforth (SAB) scheme proposed in [14]. By the nature of these schemes, we expect that the Euler scheme should produce a 2

rate of convergence 1/2, the Milstein first-order scheme a rate of 1, and the Platen-Wagner and the SAB scheme with rate of 3/2. Our simulation results indicate that this is exactly the case. Also, this example shows that in the one-dimensional case our scheme has the best performance (in terms of computing time). Due to the apparent complexity of the error analysis and the length of the paper, in this paper we only give a rigorous proof of the rate of convergence with the Euler scheme for the forward SDE. It appears that a rigorous proof for higher order schemes will be quite similar, although conceivably much more tedious. We hope our numerical results are sufficiently convincing for this purpose. As a final remark, we would like to point out that the numerical method in this paper can be, albeit technical, extended to higher-dimensional cases with a tensor-product approach. However, such a method quickly becomes non-feasible for spatial dimensions higher than three, unless some non tensor-product method, such as those based on sparse grid or lattice rules, is introduced. We plan to address the higher dimensional issue in a forthcoming work. The rest of the paper is organized as follows. In Section 2 we give the necessary preliminaries; in Section 3 we introduce the Hermite-spectral method for the PDE and perform an error analysis. In Section 4 we study the synthetic error analysis of the full numerical scheme, and in Section 5 we carry out some numerical experiments.

2

Problem Formulation and Preliminaries

Throughout this paper we assume that (Ω, F, P ) is a complete probability space, on which is defined a d-dimensional Brownian motion B = {Bt : t ≥ 0}. We shall denote FB = {FtB : t ≥ 0} to be the natural filtration generated by B, with usual P -augmentation so that it is right continuous and contains all the P -null sets in F. We consider the following FBSDE: for t ∈ [0, T ],  Z t Z t   X(t) = x + b(s, Θ(s))ds + σ(s, X(s), Y (s))dW (s);   0 0 (2.1) Z T Z T     Y (t) = ϕ(X(T )) + g(s, Θ(s))ds − Z(s)dW (s), t

t

4

where Θ = (X, Y, Z). To simplify presentation, in what follows we shall assume that all processes involved are one dimensional. Moreover, we shall make use of the following Standing Assumptions in the sequel. (A1) The functions b, σ, g, and ϕ are continuously differentiable in all variables. Moreover, if we denote these functions by a generic one ψ, then there exists a constant K > 0, such that for any t ∈ [0, T ], ψ satisfies the uniform Lipschitz condition: |ψ(t, x, y, z) − ψ(t, x0 , y 0 , z 0 )| ≤ K(|x − x0 | + |y − y 0 | + |z − z 0 |), (A2) The function b, σ, g, and ϕ satisfying the following growth conditions: for some K > 0, |b(t, x, y, z)| + |g(t, x, y, z)| ≤ K(|y| + |z|); |σ(t, x, y)| ≤ K(1 + |y|); |ϕ(x)| ≤ K. 3

(2.2)

(A3) There exist constants 0 < c < C, such that ∀ (t, x, y) ∈ [0, T ] × R2 .

c ≤ σ(t, x, y) ≤ C,

(2.3)

Remark 2.1 We remark that the assumptions (A1)–(A3) are stronger than it is necessary for the well-posedness of FBSDE (2.1). In fact, in [9] it was shown that under (A1)–(A2) but without the differentiability assumption on the coefficients the FBSDE (2.1) already possesses a unique adapted solution over arbitrarily prescribed time duration. The extra smoothness condition is only needed for our numerical scheme, and therefore not essential, in principle. Note also, however, that even with the added differentiability conditions, our assumptions are still much weaker than that of [11]. Four Step Scheme. In [20] (see also [9]) it was shown that the unique adapted solution of FBSDE (2.1) can be obtained by following steps: Step 1. Define a function z [0, T ] × R3 → R by z(t, x, y, p) = −p σ(t, x, y),

∀ (t, x, y, p).

Step 2. Using the function z above, solve the quasilinear parabolic PDE :  1   ut + σ 2 (t, x, u)uxx + b(t, x, u, z(t, x, u, ux ))ux 2 +g(t, x, u, z(t, x, u, ux )) = 0, (t, x) ∈ (0, T ) × R,   u(T, x) = ϕ(x), x ∈ R. Step 3. Using the functions u and z, solve the forward SDE: Z t Z t ˜b(s, X(s))ds + X(t) = x + σ ˜ (s, X(s))dW (s), 0

(2.4)

(2.5)

0

where ˜b = b(t, x, u(t, x), z(t, x, u, ux )) and σ ˜ = σ(t, x, u(t, x)); and Step 4. Setting Y (t) = u(t, X(t)),

and

Z(t) = σ(t, X(t), u(t, X(t)))ux (t, X(t)),

Then (X, Y, Z) is the adapted solution to (2.1). It is readily seen that if a numerical scheme is designed along the lines of Four Step Scheme, then one essentially has to deal with two separate discretization schemes: one for the PDE (2.4), and one for the (forward) SDE (2.5). We should note that the PDE (2.4) takes the form of a Cauchy problem, an artificial “cut-off” procedure is usually needed (see, e.g., [11]). To avoid such an extra step, we shall adopt a Hermite-Galerkin method with numerical integration. To do this we found it convenient to rewrite PDE (2.4) in a divergence form, mainly for the sake of numerical stability:  ³1 ´   ut + ∂x σ 2 (t, x, u)ux + ˆb(t, x, u, ux )ux 2 (2.6) +g(t, x, u, σ(t, x, u)ux ) = 0, (t, x) ∈ (0, T ) × R,   u(T, x) = ϕ(x), x ∈ R. where ˆb(t, x, y, z) = b(t, x, y, σ(t, x, y)z) − σ(t, x, y)(σx (t, x, y) + σy (t, x, y)z). 4

Hermite polynomials, Hermite functions, and their properties Throughout this paper we shall denote L2 (R) to be the space of all square integrable functions u : R 7→ R, with the inner product Z hu, vi = u(x)v(x)dx, ∀u, v ∈ L2 (R), R 2 4

R

and the norm kuk = R |u(x)|2 dx. For simplicity, in what follows, without further specification we shall always denote k · k to be the L2 (R) norm. 4

Next, for each integer k ≥ 1, we denote the differential operator ∂ k = spaces H m (R) is defined by

dk . dxk

Then, the Sobolev

H m (R) = {u | ∂ k u ∈ L2 (R) ∀ 0 ≤ k ≤ m}, 4 4 Pm with the semi-norm |u|m = k∂ m uk, and the norm kuk2m = k=0 k∂ k uk2 . We note that the operator ∂ k can be easily extended to the higher dimensional case, namely as an partial differential operator, in an obvious way. We recall that (cf. e.g., [15]) the Hermite polynomials {Hn (x)}n≥0 is defined by

Hn (x) = (−1)n exp(x2 )

dn exp(−x2 ), dxn

n = 0, 1, · · ·

(2.7)

One can easily check that Hn is the solution to the following recursive equations:   Hn+1 (x) = 2xHn (x) − 2nHn−1 (x), H 0 (x) = 2nHn−1 (x),  n H0 (x) = 1, H1 (x) = 2x; and the following orthogonality condition holds: Z ∞ 2 Hm (x)Hn (x)e−x dx = γn δmn ,

m, n ≥ 1,

(2.8)

−∞

√ where γn = π2n n!, and δmn is the Kronecker delta. For each n ≥ 1, Hn (x) is called the “Hermite polynomial of degree n”, and the Hermite function of degree n is defined by b n (x) = √ 1 e−x2 /2 Hn (x). H 2n n! b n (x) enjoy the following recurrence and orthogIt follows immediately that the Hermite functions H onal relations: r r 2 b n b b Hn+1 (x) = x Hn (x) − Hn−1 (x), n ≥ 1 (2.9) n+1 n+1 r r d b nb n+1 b Hn (x) = Hn−1 (x) − Hn+1 (x) (2.10) dx 2 2 Z ∞ √ b m (x)H b n (x)dx = πδmn . (2.11) H −∞

Moreover, we derive from 4

Hn0

b n0 + xH bn = = 2nHn−1 that H

√

b n−1 . Hence, by introducing a new 2nH

operator Du = ∂x u + xu, we obtain the following recursive relation: √ b n = 2nH b n−1 . DH 5

(2.12)

Using the operator D we can define a Hilbert space similar to the Sobolev space H m (R), which will be important for our error analysis. To begin with, for any positive integer N , let PN be the space of polynomials of degree less or equal than N . We then define HN = e−x

2

/2

4

2

PN = {e−x

/2

p(x) : p(x) ∈ PN }.

(2.13)

Next, for each N ∈ N let us denote {xi }N i=0 to be the roots of the polynomial HN +1 (x); and define √ N π2 N ! wj = , 0 ≤ j ≤ N. (2.14) (N + 1)[HN (xj )]2 The pairs {(xi , ωi )}N i=0 are known as the Hermite-Gauss quadrature points and weights, respectively (cf. e.g., [26]). Using the Hermite-Gauss quadrature points and weights we can then define a discrete inner product in L2 (R) by 4

hf, giN =

N X

∀f, g ∈ L2 (R),

f (xj )g(xj )w ˆj ,

(2.15)

j=0

where w ˆj ’s are the normalized Hermite-Gauss quadrature weights defined by √ 2 π 4 , 0 ≤ j ≤ N. w ˆj = wj exj = b 2 (xj ) (N + 1)H N A direct consequence of the Hermite-Gauss quadrature is: hu, vi = hu, viN ,

∀u · v ∈ H2N +1 .

(2.16)

For any integer m ≥ 0, let us define 4

m (R) = {u | Dm u ∈ L2 (R)}, HD r equipped with norm kDm uk. Moreover, for any real number r ≥ 0, the space HD (R) and its norm can be defined by the usual space interpolation. In order to introduce our numerical approximation we need the following two operators: 1) the “orthogonal projection operator” PN : L2 (R) → HN , defined by

hu − PN u, vi = 0,

∀u ∈ L2 (R), v ∈ HN ;

2) the “Hermite-Gauss interpolation operator” IN : C(R) → HN , defined by IN f ∈ H N ,

such that IN f (xj ) = f (xj ),

j = 0, 1, · · · , N.

(2.17)

The following basic approximation results for PN can be found in, e.g., [16]: r Lemma 2.2 For any u ∈ HD (R), N ∈ N, and 0 ≤ s ≤ r, there exits a constant C > 0, independent of N and function u, such that

kDs (u − PN u)k ≤ CN

s−r 2

kDr uk,

(2.18)

1−r 2

kDr uk.

(2.19)

and k∂x (u − PN u)k ≤ CN 6

To end this section, let us give a backward discrete Gronwall inequality, which will be important in the proof of our main theorem. A similar, but simpler version of such an inequality can be found in [22]. Lemma 2.3 Suppose that the sequence {ξk , rk , ηk }k≥0 satisfies that ξk , rk , ηk ≥ 0, and that for some constants α < 1 and β > −1, (1 − α)ξk + ηk ≤ (1 + β)ξk+1 + rk ,

0 ≤ k ≤ M − 1.

(2.20)

Then for all 0 ≤ k ≤ M , it holds that ξk +

M −1 M −1 ³ 1 + β ´M −k (1 + β)M −k−1 X (1 + β)M −k−1 X ξ + η ≤ rj . M j (1 + α)M −k 1−α (1 + α)M −k j=k

(2.21)

j=k

The proof of this lemma is based on a straightforward backward induction, and we omit it here.

3

Numerical Schemes

In this section we present our numerical schemes, both for the PDE (2.6) and for the forward SDE (2.5). Although the scheme for the general quasilinear PDE (2.6) can be treated the same way with more complicated notation, to simplify presentation we shall assume that the coefficients b and g are both independent of variable z. We note that under such a simplification the PDE (2.6) becomes  ´ ³1   ut + ∂x σ 2 (t, x, u)ux + ˆb(t, x, u, ux )ux + g(t, x, u) = 0, 2 (3.1) (t, x) ∈ (0, T ) × R,   u(T, x) = ϕ(x), x ∈ R, where ˆb(t, x, y, z) = b(t, x, y)−σ(t, x, y)(σx (t, x, y)+σy (t, x, y)z), with a slight abuse of notations. We shall design a numerical scheme for (3.1) based on a semi-implicit discretization in time and Hermitecollocation method in space; and discuss three types of scheme for the corresponding forward SDE (2.5). The error analysis for the combined scheme will be carried out in the next two sections.

3.1

Numerical Scheme for the PDE (3.1)

T For any M ∈ N, let h = M be the length of time step and set tk = kh, k = 0, 1, . . . , M . We begin with a first-order semi-implicit time discretization: Let uM (x) = ϕ(x). For k = M − 1, . . . , 0, we compute uk by solving the the following discretized version of (3.1):

´ uk+1 − uk 1 ³ + ∂x σ 2 (tk , x, uk+1 )ukx + ˆb(tk , x, uk+1 , uk+1 )uk+1 + g(tk , x, uk+1 ) = 0. x x h 2

(3.2)

In other words, at each time step, the problem is reduced to solving the following elliptic equation: uk − where

´ h ³ 2 ∂x σ (tk , x, uk+1 )ukx = fhk+1 (x), 2

lim uk (x) = 0

x→±∞

n o k+1 k+1 fhk+1 (x) = uk+1 + h ˆb(tk , x, uk+1 , uk+1 )u + g(t , x, u ) . k x x

7

(3.3)

Since the problem is set on the whole domain, it is natural to consider a Hermite-spectral method. To be more precise, let {xj }N j=0 be the Hermite-Gauss quadrature points and HN be defined in (2.13). A Hermite-collocation method for (3.3) is to find ukN ∈ HN such that ukN (xj ) −

´ h ³ 2 k+1 k ∂x σ (tk , x, uk+1 (xj ), N )∂x uN (xj ) = fh 2

j = 0, 1, · · · , N.

(3.4)

Using the identity (2.16) and integration by parts, we can rewrite (3.4) in the following variational formulation: Find ukN ∈ HN such that hukN , vN iN +

h 2 k k hσ (tk , x, uk+1 N )∂x uN , ∂x vN iN = hfh,N , vN iN , 2

∀vN ∈ HN ,

(3.5)

k k where fh,N ∈ HN and fh,N (xj ) = f k (xj ), j = 0, 1, · · · , N. We remark here that while the first order (time discretization) scheme is simple to use and analyze, in practice one often prefers higher-order schemes. These schemes could be easily constructed and often implemented in a straightforward manner, but their error analysis becomes rather tedious. In fact, as we will see in the next sections, the error analysis for the first-order semi-implicit scheme presented in (3.5) is already unpleasantly lengthy. As an example, we present the following second-order BDF Adam-Bashforth scheme which will be used in our numerical experiments without theoretical error analysis.

The second order BDF Adam-Bashforth Scheme. M −1 Let uM (x) be computed from the equation (3.2). N (x) = ϕ(x) and let uN For k = M − 2, . . . , 1, 0, we compute uk from:

³ σ2 ´ −3uk+2 + 4uk+1 − ukN N N ˆ + ∂x ∂x ukN + (2∂x uk+1 − ∂x uk+2 N N )b + g = 0, 2h 2

(3.6)

where 4

σ 2 = σ 2 (tk , x, 2uk+1 − uk+2 N N ),

4

g = g(tk , x, 2uk+1 − uk+2 N N ),

4 ˆb = ˆb(tk , x, 2uk+1 − uk+2 , 2∂x uk+1 − ∂x uk+2 ). N N N N

We note that this scheme still leads to an elliptic equation of the form (3.3) for uk at each time step. In particular, this second-order scheme should be used together with the 32 th-order SDE scheme presented below.

3.2

Numerical schemes for the SDE (2.5)

We now turn our attention to the discretization of the forward SDE (2.5). We begin with the simplest one, known as the “forward Euler-scheme”. Assuming that for each N ∈ N we have obtained the approximating solution to the PDE (3.1), denoted by ukN , k = 0, 1, · · · , M , at each time tk = kh we define recursively the approximate solution to the SDE (2.5) by: ½ N X0 = x, (3.7) N = XkN + b(tk , XkN , ukN (XkN ))h + σ(tk , XkN , ukN (XkN ))∆k W, k = 0, 1, · · · M. Xk+1 4

where ∆k W = W (tk + h) − W (tk ). For notational simplicity in what follows we shall simply denote Xk = XkN , when context is clear. 8

It is well-understood that the Euler scheme is easy to implement and has the minimum requirements on the coefficients. Furthermore, if we denote the true solution to the PDE (3.1) by u, denote 4 ˆ by uk (x) = u(tk , x), and define an intermediate approximation X ½

ˆ 0 = x, X ˆ k+1 = X ˆ k + b(tk , X ˆ k , uk (X ˆ k ))h + σ(tk , X ˆ k , uk (X ˆ k ))∆k W, k = 0, 1, · · · M, X

(3.8)

then by the Fundamental Convergence Theorem (cf. e.g., [18]), the mean-square error of the Euler ˆ k |2 ∼ h, where X is the true solution of (2.5) (namely, the rate of converges scheme is E|X(tk ) − X 1 is 2 ). In order to obtain a higher order rate of convergence, one has to use a higher order scheme for the forward SDE. However, we should note that while in theory arbitrarily high-order schemes for (2.5) can be constructed with Taylor-Itˆo type expansions (see, e.g., [18]), the complexity of these schemes increases drastically. Consequently it often becomes too “expensive” in computational terms to implement. We will consider the following higher order schemes in our numerical experiments to test the numerical accuracy, and to compare our method with existing results. A. Milstein Scheme (cf. [18], [24]). This is a well-known scheme with first-order convergence rate. The recursive relation, adapted to our case, is given as follows.  X0 = x,    X k k  k+1 = Xk + b(tk , Xk , uN (Xk ))h + σ(tk , Xk , uN (Xk ))∆k W (3.9) 1  + σ(tk , Xk , ukN (Xk )){σx (tk , Xk , ukN (Xk ))+σu (tk , Xk , ukN (Xk ))∂x ukN (Xk )}   2  ×(∆2k W − h), k = 0, 1, . . . , M − 1; 4

where ∆2k W = (W (tk+1 ) − W (tk ))2 . The next two schemes requires higher order regularity of the coefficients. We shall assume that all such requirements, whenever needed, are fulfilled without further specification. B. Platen-Wagner Scheme (cf. [18]). Using the idea of Taylor-Itˆo expansion up to order 32 , and assuming that the coefficients are actually twice continuously differentiable, Platen and Wagner proposed the following scheme: for k = 0, 1, . . . , M − 1,  X0 = x,    1 ³ 0 1 2 00 ´ 2 1  bb + σ b h Xk+1 = Xk + bh + σ∆k W + σσ 0 {∆2k W − h} + b0 σ∆k Z + (3.10) 2 2 2  ³ ´ ³ ´  1 1 1   + bσ 0 + σ 2 σ 00 {h∆k W − ∆k W } + σ σσ 00 + (σ 0 )2 { ∆2k W − h}∆k W 2 2 3 R 4 4 t where ∆2k W = (W (tk+1 ) − W (tk ))2 , ∆k Z = tkk+1 [W (t) − W (tk )]dt, b = b(tk , Xk , ukN (Xk )), σ = σ(tk , Xk , ukN (Xk )), and for ϕ = b, σ, ϕ0 = ϕ + x + ϕu ∂x ukN ,

4

ϕ00 = ϕxx + 2ϕxu ∂x ukN + ϕuu (∂x ukN )2 + ϕu ∂x2 ukN .

(3.11)

Here the partial derivatives of b and σ should be evaluated at (tk , Xk , ukN (Xk )), and the approximate functions ukN ,

duk N dx

and

d2 uk N dx2

should be evaluated at point (Xk ). 9

Rt We should note that ∆k z = tkk+1 W (t) − W (tk )dt ∼ N (0, 13 h3 ) , and the covariance of ∆k W ³ ´ and ∆k W is E ∆k W ∆k Z = 12 h2 . The presence of the ∆k Z obviously complicates the analysis of the scheme. This will be even more so when the order of expansion increases. Treating these terms effectively will become more important. C. The Stochastic Adam-Bashforth (SAB) Scheme (cf. [14]). It is easily seen from (3.11) that Platen-Wagner’s 23 -order scheme requires that b and σ are both twice differentiable. Recently, Brian D. Ewald and Roger Temam [14] constructed a stochastic version of the Adam-Bashforth scheme which does not involve the second derivatives of coefficient function b but still achieves the 23 -order convergence rate. This Stochastic Adam-Bashforth (SAB) scheme, adopted for an one-dimensional diffusion of the form: dXt = b(t, Xt )dt + σ(t, Xt )dW (t), can be rewritten as follows: Xk+2 = Xk+1 +

∆t 3 [3b(tk+1 , Xk+1 ) − b(tk , Xk )] − ∆tAk (tk , Xk ) + Bk (tk , Xk ), 2 2

(3.12)

[σbx ](t, x)∆k W,

(3.13)

where Ak (t, x) Bk (t, x)

=

1 = σ(t, x)∆k W + [σt + bσx + σ 2 σxx ](t, x)I(0,1) + [σbx ](t, x)I(1,0) , 2 +[σσx ](t, x)I(1,1) + [σ((σx )2 + σσxx )](t, x)I(1,1,1) ;

with ∆k W = Wtk+1 − Wtk and h = tk+1 − tk as before, and the random coefficients I’s are defined by Z tk+2 I(0,1) = 2h∆k+1 W − [W (s) − W (tk+1 )]ds Z I(1,0)

tk+1 tk+2

= h∆k W +

[W (s) − W (tk+1 )]ds

(3.14)

tk+1

and I(1, 1), I(1, 1, 1) are the iterated Itˆo integrals Z I(1,1)

=

I(1,1,1)

=

tk+2

Z

Z

t

tk+1

Z

t

dW (s)dW (t) − tk Z tk+2

tk Z t

Z

dW (s)dW (t) tk

tk

s

Z

tk+1

Z tZ

s

dW (r)dW (s)dW (t)

dW (r)dW (s)dW (t) − tk

tk

tk

tk

(3.15)

tk

tk

To calculate the iterated Itˆo integrals I(1,...,1) , the useful formulae of Itˆo [17] are often been employed. More precisely, let us define a scaled Hermite polynomial hn by √ Hn (x) = 2n/2 hn ( 2x),

n = 0, 1, 2, · · · ,

then we have (cf. [26]) 1

2

hn (x) = (−1)n e 2 x

dn − 1 x2 (e 2 ); dxn 10

n = 0, 1, 2, · · · ,

and for and g ∈ L2 ([a, b]), b > a ≥ 0, it holds that Z b Z tn Z t2 θ n! ··· g(t1 )g(t2 ) · · · g(tn )dWt1 · · · dWtn = kgkn hn ( ), kgk a a a

(3.16)

where Z kgk = kgkL2 ([a,b])

and

θ=

b

g(t)dWt . a

Thus the coefficients I(1,1) and I(1,1,1) in (3.15) can be calculated explicitly as ´ 1³ 2 I(1,1) = ∆k+1 W + 2∆k+1 W ∆k W − h 2 ´ 1³ I(1,1,1) = (∆k+1 W + ∆k W )3 − (∆k W )3 − 6h∆k+1 W − 3h∆k W . 6 Using the above relations, we can rewrite the forward SDE (2.5) as X0

=

x,

Xk+2

=

Xk+1 +

3 2 -order

(3.17)

SAB scheme for the one dimensional

i hh 3b(tk+1 , Xk+1 , uk+1 (X )) − b , k+1 N 2 n o 3 1 − hσb0 ∆k W + σ∆k+1 W + σσ 0 (∆k+1 W )2 + 2(∆k+1 W )(∆k W ) − h 2 2 ³ 1 2 00 ´ 0 + bσ + σ σ {2h∆k+1 W − ∆k+1 Z} + σb0 {h∆k W + ∆k+1 Z} 2 ´n o 1 ³ 00 + σ σσ + (σ 0 )2 (∆k+1 W + ∆k W )3 − (∆k W )3 − 6h∆k+1 W − 3h∆k W 6 where k = 0, 1, . . . , M − 1 and ,similarly, b =

b(tk , Xk , ukN (Xk ))

σ

=

b0

=

σ0

=

σ 00

=

σ(tk , Xk , ukN (Xk )) ∂b ∂b dukN + ∂x ∂u dx ∂σ ∂σ dukN + ∂x ∂u dx ∂2σ ∂ 2 σ dukN ∂ 2 σ ³ dukN ´2 ∂σ d2 ukN +2 + + 2 ∂x ∂x∂u dx ∂u2 dx ∂u dx2

where the partial derivatives of b and σ are evaluated at (tk , Xk , ukN (Xk )), and the approximate duk

d2 uk

functions ukN , dxN and dx2N are evaluated at point (Xk ). Finally, the components Y and Z of the solution to (2.1) are approximated as Yk = ukN (Xk ),

4

Zk = σ(tk , Xk , Yk )

∂ukN (Xk ), ∂x

(3.18)

Error analysis for the PDE approximation

In this section we carry out an error analysis for the approximation of PDE (3.1). In order to simplify the presentation, we shall only consider the following Hermite-Galerkin scheme for (3.2): 11

Find ukN ∈ HN such that ³ ´ uk+1 − ukN 1 k N , vN i + h∂x σ 2 (tk , x, uk+1 N )(uN )x , vN i h 2 k+1 k+1 +hˆb(tk , x, uN , (uk+1 ) + g(tk , x, uk+1 N )x )∂x (u N ), vN i = 0,

h

(4.1) ∀vN ∈ HN .

We note that the only difference between this Hermite-Galerkin scheme and the Hermite-collocation scheme (3.5) is that in the latter the continuous inner product is replaced by the discrete inner product. It is well-known that for N sufficiently large, the difference between the two approaches are negligible. Hereafter, we shall use “A . B” to mean that there exists a constant C > 0 independent of N or h such that A ≤ CB. Our main result in this section is the following theorem. Theorem 4.1 Assume (A1)—(A3). Let u and ukN be the solutions of (3.1) and (4.1), respectively. m m Assume further that u, φ ∈ L∞ (0, T ; HD (R)) and ∂t u ∈ L2 (0, T ; HD (R)) with some m > 1. Then, there exists h0 > 0 such that for h ≤ h0 , the scheme (3.2) is stable, and the following error estimates hold: for each k = 0, 1, · · · , M , M ³ X ´ 12 ku(tk , ·) − ukN k + h k∂x (u(tk , ·) − ukN )k2 j=k

.

N

1−m 2

´ ³ m ) + kut kL2 (0,T ;H m ) + kϕk ∞ m+1 kukL∞ (0,T ;HD ) + h. L (0,T ;H D D

Proof. Let E k , k = M − 1, . . . , 1, 0, be the truncation error defined by ´ ³ ´ 1³ 4 E k (·) = u(tk+1 , ·) − u(tk , ·) + b(tk , ·, u(tk+1 , ·), ∂x u(tk+1 , ·)) ∂x u(tk+1 , ·) h ³1 ´ +∂x σ 2 (tk , ·, u(tk+1 , ·))∂x u(tk , ·) + g(tk , ·, u(tk+1 , ·)), 2

(4.2)

where u(tk , ·) is the exact solution to equation (3.1). In what follows we denote eˆkN = PN u(tk , ·)−ukN , e˜kN = u(tk , ·) − PN u(tk , ·), and thus ekN = u(tk , ·) − ukN = eˆkN + e˜kN . Next, multiplying vN ∈ HN on both sides of (4.2), using integration by parts, and then subtracting (4.1), we obtain that hE k , vN i =

1 k+1 1 he − ekN , vN i − h σ 2 (tk , ·, u(tk+1 , ·))∂x u(tk , ·) , ∂x vN i h N 2 1 2 k +h σ (tk , ·, uk+1 N )∂x uN , ∂x vN i 2 k+1 k+1 +hb(tk , ·, u(tk+1 , ·), ∂x u(tk+1 , ·))) ∂x u(tk+1 , ·) − b(tk , ·, uk+1 , vN i N , ∂x uN )∂x uN +hg(tk , ·, u(tk+1 , ·)) − g(tk , ·, uk+1 N ) , vN i.

Or equivalently, hE k , vN i =

1 k+1 1 k he − ekN , vN i − h σ 2 (tk , ·, uk+1 N )∂x eN , ∂x vN i h N 2 1 1 −h∂x u(tk , ·) , σ 2 (tk , ·, u(tk+1 , ·))∂x vN − σ 2 (tk , ·, uk+1 N )∂x vN i 2 2 k+1 k+1 +hb(tk , ·, u(tk+1 , ·), ∂x u(tk+1 , ·)) ∂x u(tk+1 , ·) − b(tk , ·, uk+1 , vN i N , ∂x uN )∂x uN +hg(tk , ·, u(tk+1 , ·)) − g(tk , ·, uk+1 N ) , vN i. 12

Since ekN = e˜kN + eˆkN , we have 1 k+1 hˆ e − eˆkN h N

, =

1 vN i − h σ 2 (tk , ·, uk+1 ˆkN , ∂x vN i N )∂x e 2 1 k+1 1 − h˜ eN − e˜kN , vN i + h σ 2 (tk , ·, uk+1 ˜kN , ∂x vN i (4.3) N )∂x e h 2 1 1 +h∂x u(tk , ·) , σ 2 (tk , ·, u(tk+1 , ·))∂x vN − σ 2 (tk , ·, uk+1 N )∂x vN i 2 2 k+1 k+1 , vN i −hb(tk , ·, u(tk+1 , ·), ∂x u(tk+1 , ·)) ∂x u(tk+1 , ·) − b(tk , ·, uk+1 N , ∂x uN )∂x uN k −hg(tk , ·, u(tk+1 , ·)) − g(tk , ·, uk+1 N ) , vN i + hE (·) , vN i.

We now choose vN = −hˆ ekN , and denote I1

= −hg(tk , ·, u(tk+1 , ·)) − g(tk , ·, uk+1 N ) , vN i;

I2

k+1 k+1 = −hb(tk , ·, u(tk+1 , ·), ∂x u(tk+1 , ·)) ∂x u(tk+1 , ·) − b(tk , ·, uk+1 , vN i; N , ∂x uN )∂x uN 1 1 = h∂x u(tk , ·) , σ 2 (tk , ·, u(tk+1 , ·))∂x vN − σ 2 (tk , ·, uk+1 N )∂x vN i; 2 2 1 = h σ 2 (tk , ·, uk+1 ˜kN , ∂x vN i; N )∂x e 2 1 k+1 = − h˜ e − e˜kN , vN i; h N = hE k (·) , vN i.

I3 I4 I5 I6

We shall estimate these terms separately under the assumptions (A1)-(A2). To begin with, noting the Lipschitz property of g and applying the Cauchy-Schwarz inequality, we have, for some constant C1 > 0, and any ε > 0, I1

k+1 C1 h|u(tk+1 , ·) − uk+1 N | , |vN |i = C1 h|eN | , |vN |i C1 h k+1 2 C1 ε C1 h k+1 2 C1 h k+1 2 C1 ε ke k + kvN k2 ≤ kˆ eN k + k˜ eN k + kvN k2 . 2ε N 2h ε ε 2h

≤ ≤

(4.4)

Similarly, we have C4 h k 2 C4 ε k˜ e k + k∂x vN k2 . (4.5) 2ε N 2h Second, note that under assumptions (A1) and (A2), the solution to (3.1) must have bounded first order derivative ∂x u (see, e.g., [9]). Thus, one has I4 ≤ C4 h|∂x e˜kN | , |∂x vN |i ≤

I2

=

k+1 −hb(tk , ·, u(tk+1 , ·), ∂x u(tk+1 , ·)) ∂x u(tk+1 , ·) − b(tk , ·, uk+1 N , ∂x uN ))∂x u(tk+1 , ·) , vN i k+1 k+1 −hb(tk , ·, uk+1 N , ∂x uN ))(∂x u(tk+1 , ·) − ∂x uN ) , vN i

≤

k+1 k+1 Kk∂x uk∞ h|ek+1 N |+|∂x eN |, |vN |i+kbk∞ h|∂x eN |, |vN |i

≤

K k∂x uk∞ h|ˆ ek+1 ˆk+1 ek+1 N | + |∂x e N | + |˜ N | k+1 k+1 +kbk∞ h|∂x eˆN | + |∂x e˜N | , |vN |i

≤

+

k+1 |∂x e˜N |

(4.6)

, |vN |i

C2 C2 ε C2 h k+1 2 2 2 2 (kˆ eN k + k˜ kvN k2 + C2 ε1 hk∂x eˆk+1 kvN k2 . ek+1 ˜k+1 N k + N k + k∂x e N k )+ ε h ε1 h

where C2 = 4 max{K k∂x uk∞ , kbk∞ } and ε1 is to be determined later. Similarly, we have C3 h k+1 2 C3 ε ke k∂x vN k2 k + 2ε N 2h C3 h k+1 2 C3 ε 2 (kˆ eN k + k˜ k∂x vN k2 . ek+1 N k )+ ε 2h

I3 ≤ L k∂x uk∞ h|ek+1 N | , |∂x vN |i ≤ ≤

13

(4.7)

k+1 ˜ k = − 1 (˜ where C3 = L k∂x uk∞ . Finally, if we denote E − e˜kN ), then we have N h eN

I5 ≤

h ˜k 2 ε kE k + kvN k2 , 2ε N 2h

and I6 ≤

h ε kE k (·)k2 + kvN k2 . 2ε 2h

(4.8)

To obtain the desired estimate, let us now look at the left hand side of (4.3) with vN = −h eˆkN . Note that 1 k+1 1 k 2 2 2 hˆ eN − eˆkN , −h eˆkN i = (kˆ e k − kˆ ekN − eˆk+1 ek+1 N k + kˆ N k ), h 2 N 1 −h∂x eˆkN , σ 2 (tk , ·, uk+1 ˆkN k2 , N )∂x vN i ≥ chk∂x e 2

(4.9) (4.10)

where c > 0 is a lower bound of σ 2 , with a slight abuse of notation (compare to the constant c > 0 in (A3)). Thanks again the Schwarz inequality, and denoting C = 5 max1≤i≤6 {Ci } + 1, we derive from (4.4)—(4.9) that kˆ ekN k2 + chk∂x eˆkN k2 ≤ (

Ch 1 k+1 2 Ch 1 k 2 2 + )kˆ e k + Cε1 hk∂x eˆk+1 + Cεh + εh + )kˆ e k N k +( ε 2 N ε1 2 N Ch k+1 2 2 +Cεhk∂x ekN k2 + (k˜ eN k + k∂x e˜k+1 (4.11) N k ) ε Ch k 2 h k h ˜k 2 + k˜ eN k + kE (·)k2 + kE k . ε ε ε N

Note that the truncation error E k (·) is of first order, i.e., kE k (·)k2 ≤ C8 h2 .

(4.12)

We see that (4.11) can now be rewritten as ³1 ³ C ´´ k 2 − h Cε + ε + kˆ eN k + (c − Cε)hk∂x eˆkN k2 2 ε1 ³ Ch 1 ´ 2 2 ≤ + kˆ ek+1 ˆk+1 N k + Cε1 hk∂x e N k ε 2 ´ h³ 2 2 2 ˜k 2 + Ck˜ ek+1 ekN k2 + Ck∂x e˜k+1 N k + Ck˜ N k + kEN k + C8 h . ε

(4.13)

2c−1 2C so that c − Cε = 1/2). 4 4 C9 = max{C + 1, C8 }, Cε =

Now, let us fix ε small enough so that c − Cε > 0 (for example, ε =

Multiplying both sides of (4.13) by 2, choosing ε1 = c−Cε 2C , setting C 2(Cε + ε + ε1 ), and finally assuming that the time step h satisfies the conditions hCε
0, and the ε ³ 2Ch ´(M −k−1) ε +1 is bounded from below by α2 > 0. We note that both constants α1 and sequence 1+h ˜ C α2 depend only on the constants c, C and T for fixed ε, and the bounds are all uniformly in k. We now take a closer look at all the errors on the right hand side above. First, by virtue of Lemma 2.2 with s = 0, 1, respectively, and r = m, we have 2 k˜ ek+1 ekN k2 ≤ C7 N −m kDm uk2 , N k + k˜

2 1−m k∂x e˜k+1 kDm uk2 . N k ≤ C7 N

(4.17)

˜ k k. Let e˜N (t, ·) = u(t, ·) − PN u(t, ·). Then, It remains to estimate kE N 1 k+1 1 k ˜N eN − e˜kN ) = − E = − (˜ h h

Z

tk+1

∂t e˜N (t, ·)dt. tk

Applying the Schwarz inequality and (2.18) again we obtain that k 2 ˜N kE k ≤

Consequently, recalling that hCε
0 to be a generic constant depending only on the constants in (A1)–(A3) and time duration T > 0, and is allowed to vary from line to line. Combining (5.8)–(5.11), noting (5.7), and assuming without loss of generality that h < 1, it follows easily from (5.6) that, for k = 1, 2, · · · , M , ³ ´ ˆ k+1 − Xk+1 |2 ≤ (1 + h)E|X ˆ k − Xk |2 + C(h2 + h) E|X ˆ k − Xk |2 + h E|X ³ ´ ³ ´ ˆ k − Xk |2 + h + C(h2 + h) ku(tk , ·) − ukN k2 + k∂x (u(tk , ·) − ukN )k2 + Ch E|X ³ ´ ˆ k − Xk |2 + Ch2 + Ch ku(tk , ·) − ukN k2 + k∂x (u(tk , ·) − ukN )k2 . ≤ (1 + Ch)E|X ˆ 0 = X0 = x, we obtain, for k = Apply the forward discrete Gronwall inequality , and notice X 1, 2, · · · , M , ˆ k − Xk |2 E|X

ˆ 0 − X0 |2 + C(1 + Ch)k−1 ≤ (1 + Ch)k E|X

k X

h2

j=1

+C(1 + Ch)k−1 h

k ³ X

ku(tj , ·) − ujN k2 + k∂x (u(tj , ·) − ujN )k2

´

j=1

h i ≤ C(1 + Ch)M T h + N 1−m Υ(u, φ) . In the above we used the fact that k ≤ M = T /h. Further, noting T = hM we also have (1 + Ch)M = (1 +

CT M ) → eCT , M

as M → ∞,

consequently we must have ˆ k − Xk |2 . h + N 1−m Υ(u, φ), E|X where C is a generic constant that is independent of the choice of the number of time steps M and the number of Hermite functions N . This, together with (5.4), yields that final mean square error estimate for X(tk ) − Xk : E|X(tk ) − Xk |2 . h + N 1−m Υ(u, φ).

(5.13)

The error estimates for E|Y (tk )−Yk |2 and E|Z(tk )−Zk |2 then follows easily from the relation (3.18), (5.1) and the estimate (5.13). The proof is now complete. Remark 5.2 From the proof of Theorem 5.1 we see that if we combine the higher order schemes for both PDE (2.6) and the resulting decoupled SDE (2.5), then we will obtain a higher order rate of convergence. The proof will be completely similar, but conceivably much more lengthy and tedious. We shall discuss this issue with extensive numerical examples in the next section.

6 6.1

Numerical Simulations and Conclusions Implementation details

Let {hj (x)}j=0,1,...,N be the Lagrange interpolation polynomials in PN associated with the Gauss 2 ˆ j (x) = hj (x) e−x2 /2 , then points {xj }j=0,1,...,N , which are zeros of HN +1 (x). We define h −x /2 j

e

ˆ j (xi ) = δij , h

ˆ j : j = 0, 1, . . . , N } HN = span{h 18

so any function g ∈ HN can be written as g(x) =

N X

ˆ j (x). g(xj )h

j=0 0

We derive easily from hj (xi ) that (

2

0

ˆ (xi ) = h j

e−xi /2 e

−x2j /2

0

(hj (xi ) − xi δij ) =

b N (xi ) H b N (xj ) , (xi −xj )H

0,

i 6= j i=j

ˆ (xi ) can be computed in a stable way for N large. Hence,h j Let us denote 0

u ¯ = (ukN (x0 ), ukN (x1 ), . . . , ukN (xN ))T , f¯ = (f k+1 (x0 ), f k+1 (x1 ), . . . , f k+1 (xN ))T 0

0

ˆ ,h ˆ )N = sij = (α(x)h j i

N X

0

0

ˆ (xl )h ˆ (xl )w α(xl )h ˆl + α(xi )w ˆi δij j i

l=0

with α(xi ) = system:

k+1 h 2 2 σ (tk , xi , uN (xi )).

Then, the equation (3.5) is reduced to the following linear (S + W )¯ u = W f¯,

where W = diag(w ˆ0 , w ˆ1 , . . . , w ˆN ), S = (sij )i,j=0,1,...,N

6.2

Numerical results and discussions

In order to make sensible comparisons to the existing results, we shall consider an example proposed by Milstein and Tretyakov [24], and use our spectral method to carry out the numerical approximation for u(t, x), and then use different stochastic numerical schemes for the resulting decoupled SDE. In particular, we shall use the standard Euler scheme, first-order Milstein scheme, 23 -order Platen-Wagner strong scheme, as well as 23 order stochastic Adam-Bashforth scheme (SAB) scheme, respectively. Example. (Milstein-Tretyakov [24]). Consider the following FBSDE: s  1 + X2 1 + 2Y 2 X(1 + X 2 )   dt + dX =  2 dW (t),  2 3 2 2  (2 + X ) 2+X 1 + Y + exp(− 2X t+1 )      dY = −g(t, X, Y )dt − f (t, X, Y )Zdt + ZdW (t),

X(0) = x,

X 2 (T ) ), Y (T ) = exp(− T +1

where g(t, x, u) = f (t, x, u) =

1 x2 h 4x2 (1 + x2 ) ³ 1 + x2 ´2 ³ 2x2 ´ x2 i exp(− ) + 1 − − , t+1 t+1 (2 + x2 )3 2 + x2 t+1 t+1 s 2x2 1 + u2 + exp(− t+1 ) x . (2 + x2 )2 1 + 2u2

19

(6.1)

Then the corresponding Cauchy problem has the form, for t < T, x ∈ R  ∂u 1 ³ 1 + x2 ´2 1 + 2u2 ∂ 2 u 2x(1 + x2 ) ∂u   0 = + +  2  2x  ∂t 2 2 + x2 1 + u2 + exp(− t+1 (2 + x2 )3 ∂x ) ∂x2    2x2 ´ 1 x2 h 4x2 (1 + x2 ) ³ 1 + x2 ´2 ³ x2 i  + exp(− ) + 1 − − ,   t+1 t+1 (2 + x2 )3 2 + x2 t+1 t+1      u(T, x) = exp(− x2 ).

(6.2)

T +1

It is easy to verify that the solution to problem (6.2) is u(t, x) = exp(−

x2 ). t+1

(6.3)

Since Y (t) = u(t, X(t)), plugging Y (t) into the forward equation in (6.1) we get dX =

X(1 + X 2 ) 1 + X2 dt + dW (t), 2 3 (2 + X ) 2 + X2

X(0) = x.

(6.4)

One can then easily check that the solution to this equation can be expressed as X(t) = Λ(x + arctan x + W (t)), where the function Λ(z) is defined by the equation Λ(z) + arctan Λ(z) = z.

(6.5)

Differentiating (6.5) with respect to z, we get Λ0 (z) =

1 + Λ2 (z) > 0 for all 2 + Λ2 (z)

z ∈ R, and Λ00 (z) =

2Λ(1 + Λ2 ) . (2 + Λ2 )3

(6.6)

which implies that Λ(z) is an one-to-one function. Hence X(0) = Λ(x + arctan x) = x. Furthermore, due to the Itˆo formula, we have dX

1 = Λ0 (x + arctan x + W (t))dW (t) + Λ00 (x + arctan x + W (t))dW (t) 2 X(1 + X 2 ) 1 + X2 = dt + dW (t), (2 + X 2 )3 2 + X2

Hence, the solution to (6.1) is X(t)

=

Λ(x + arctan x + W (t)), ³ X2 ´ , Y (t) = exp − t+1 ³ 2X(1 + X 2 ) X2 ´ Z(t) = − exp − . 2 (t + 1)(2 + X ) t+1

We first carry out numerical tests on the parabolic equation (6.2) using the second-order (in time) scheme (3.6). In Tables 1 to 3, we tabulate the errors max ku(tk , x) − ukN (x)kL2 (R) , max k ∂u ∂x (tk , x) − d k 2 dx uN (x)kL (R) ,

2

k

2

k

d k and max k ∂∂xu2 (tk , x) − dx 2 uN (x)kL2 (R) , respectively. We note that it is important to k

measure the errors of the approximation ukN (x) to the first and second derivatives of u(tk , x) since 20

Table 1: max ku(tk , x) − ukN (x)kL2 (R) k

h\N 0.5 0.2 0.05 0.02 0.005

64 0.0049 0.0010 1.5391e-004 1.5389e-004 1.5388e-004

70 0.0049 0.0010 9.0229e-005 8.3270e-005 8.3266e-005

100 0.0049 0.0010 7.8519e-005 1.3177e-005 3.9669e-006

128 0.0049 0.0010 7.8554e-005 1.3167e-005 8.4446e-007

150 0.0049 0.0010 7.8555e-005 1.3168e-005 8.4360e-007

RATE 1.7344 1.8351 1.9492 1.9822

the first-order Milstein scheme (3.9) and 32 -order Platen-Wagner scheme (3.10) need to use these derivatives. In the last column of these tables, the rates of convergence in time with N = 150 are reported. We observe that essentially second-order accuracy in time and exponential convergence in space are achieved for all three quantities.

Table 2: max k ∂u ∂x (tk , x) − k

h\N 0.5 0.2 0.05 0.02 0.005

64 0.0075 0.0018 0.0017 0.0017 0.0017

70 0.0074 0.0017 9.6178e-004 9.6179e-004 9.6178e-004

100 0.0074 0.0017 1.3741e-004 5.4862e-005 5.4859e-005

128 0.0074 0.0017 1.3685e-004 2.3165e-005 3.7980e-006

2

Table 3: max k ∂∂xu2 (tk , x) − k

h\N 0.5 0.2 0.05 0.02 0.005

64 0.0230 0.0188 0.0188 0.0188 0.0188

70 0.0226 0.0112 0.0112 0.0112 0.0112

100 0.0224 0.0055 7.7070e-004 7.7061e-004 7.7058e-004

128 0.0224 0.0055 4.5192e-004 7.8042e-005 6.0488e-005

d k 2 dx uN (x)kL (R)

150 0.0074 0.0017 1.3685e-004 2.3152e-005 1.4925e-006

RATE 1.6052 1.8174 1.9391 1.9777

d2 k 2 dx2 uN (x)kL (R)

150 0.0224 0.0055 4.5174e-004 7.7116e-005 8.0385e-006

RATE 1.5326 1.8029 1.9293 1.6310

Next, we examine the errors of the full simulation for the FBSDE with four different schemes presented in Section 3 for the forward SDE and with the Hermite-collocation scheme (3.5). For the sake of comparison with the results in [24], we also computed the same problem with a Monte Carlo simulation (with S = 1000 independent realizations of X(T ) and XN ) for the forward SDE. The results of the Euler scheme with Monte Carlo simulation are reported in Table 6.2. The averages presented in the table are computed as follows: r S ´2 DS 1 X ³ (k) (k) 2 X (T ) − XN ±2 E(X(T ) − XN ) = S S k=1

21

where DS =

S S ³ ´4 h 1 X ´2 i2 1 X ³ (k) (k) (k) X (T ) − XN X (k) (T ) − XN − S S k=1

k=1

Hence, h i1/2 E(X(T ) − XN )2

=

+

v u S ³ ´2 u1 X (k) t X (k) (T ) − XN S k=1 v v r u u S ³ S ´2 ´2 u1 X DS u 1 X ³ (k) (k) (k) t −t X (k) (T ) − XN X (T ) − XN ±2 S S S k=1

k=1

(k)

X (k) (T ) and XN are independent realizations of X(T ) and XN , respectively, k = 1, 2, . . . , S = 1000. Table 4: Euler Scheme with 1000 Monte Carlo simulation realizations with N = 150 h h 0.5 0.2 0.05 0.02 0.005

E(X(T ) − XN )2

i1/2

0.2453 ± 0.0121 0.1603 ± 0.0083 0.0745 ± 0.0037 0.0479 ± 0.0024 0.0259 ± 0.0013

h Rate 0.4643 0.5527 0.482 0.4435

E(Y (T ) − YN )2

i1/2

0.0325 ± 0.0022 0.0213 ± 0.0015 0.0101 ± 0.0007 0.0064 ± 0.0004 0.0036 ± 0.0002

h

E(Z(T ) − ZN )2

i1/2

0.0121 ± 7.9404e − 004 0.0079 ± 5.0520e − 004 0.0037 ± 2.2556e − 004 0.0023 ± 1.4899e − 004 0.0013 ± 7.6865e − 005

For conciseness, we now omit the Monte Carlo errors, which, as can be seen from Table 6.2, are significantly small than the approximation errors. In what follows, we denote [E(X(T ) − XN )2 ]1/2 by X, [E(Y (T ) − YN )2 ]1/2 by Y and [E(Z(T ) − ZN )2 ]1/2 by Z. In the following three tables, we report these errors by using the three schemes with 1000 independent realizations of X(T ) and XN .

Table 5: First Order SDE Scheme with N = 150 h .5 .2 .05 .02 .005

X 0.0937 0.0394 0.0095 0.0039 9.5023e-004

Rate 0.9455 1.0261 0.9717 1.0186

Y 0.0131 0.0052 0.0013 5.3451e-004 1.2996e-004

Rate 1.0084 1.000 0.9700 1.0201

Z 0.0048 0.0021 4.7994e-004 1.9852e-004 4.8413e-005

Rate 0.9022 1.0647 0.9634 1.0179

Observe that all three schemes produce expected convergence rates. The 32 -order SAB scheme is slightly less accurate than the 32 -order strong scheme but the convergence rates of these two schemes are essentially the same. 22

Table 6: h .5 .2 .05 .02 .005

X 0.0425 0.0115 0.0015 3.3583e-004 3.9962e-005

Rate 1.4266 1.4693 1.6333 1.5355

7

X 0.1301 0.0400 0.0048 0.0013 1.5749e-004

Rate 1.2872 1.5294 1.4256 1.5226

order strong scheme with N = 200

Y 0.0041 0.0015 1.7410e-004 3.6401e-005 7.6776e-006

Table 7: h .5 .2 .05 .02 .005

3 2

3 2

Rate 1.0974 1.5535 1.7080 1.1226

Z 0.0025 5.4051e-004 7.5482e-005 1.7780e-005 2.2442e-006

Rate 1.6714 1.4201 1.5779 1.4930

order SAB scheme with N = 200

Y 0.0177 0.0053 6.6087e-004 1.6757e-004 2.0907e-005

Rate 1.3160 1.5018 1.4975 1.5014

Z 0.0063 0.0020 2.3552e-004 6.5992e-005 7.9073e-006

Rate 1.2522 1.5430 1.3885 1.5305

Concluding remarks

We presented a numerical method for a class of forward-backward stochastic differential equations (FBSDEs). The method is based on the Four Step Scheme with a Hermite-spectral method to approximate the solution to the decoupled quasilinear PDE on the whole space. The use of Hermitespectral method not only avoids the use of artificial far field boundary conditions but also leads to spectrally accurate results in space. We carried out a rigorous error analysis for a fully discretized scheme for the FBSDEs with a first-order scheme in time and a Hermite-spectral scheme in space, and indicated that similar analysis can be extended to higher-order schemes in time. We presented detailed numerical comparisons between several schemes for the resulting decoupled forward SDE and showed that the stochastic version of the Adams-Bashforth scheme coupled with the Hermitespectral method leads to a convergence rate of 32 (in time). Although the analysis and computation is performed for the one-dimensional case, it is clear that similar analysis can be carried out for higher-dimensional cases, while computationally a direct tensor-product extension of the algorithm quickly becomes prohibitive as the dimension increases since . In a forthcoming work, we plan to introduce new elliptic solvers based on lattice rules and sparse grids which would allow us to handle FBSDEs with relatively large dimensional problems.

References [1] V. Bally. Approximation scheme for solutions of BSDE. Backward stochastic differential equations (Paris, 1995–1996), 364:177–191, 1997. [2] V. Bally and G. Pag`es. A quantization algorithm for solving multi-dimensional discrete-time optimal stopping problems. Bernoulli, 9(6):1003–1049, 2003. [3] V. Bally, G. Pag´es, and J. Printems. A quantization tree method for pricing and hedging multidimensional American options. Math. Finance, 15(1):119–168, 2005.

23

[4] C. Bender and J. Zhang. Preprint, 2006.

Time discretization and markovian iteration of coupled fbsdes.

[5] B. Bouchard and N. Touzi. Discrete-time approximation and Monte-Carlo simulation of backward stochastic differential equations. Stochastic Process. Appl., 111(2):175–206, 2004. [6] P. Briand, B. Delyon, and J. M´emin. Donsker-type theorem for BSDEs. Electron. Comm. Probab., 6:1–14 (electronic), 2001. [7] D. Chevance. Numerical methods for backward stochastic differential equations. In Numerical methods in finance, Publ. Newton Inst., pages 232–244. Cambridge Univ. Press, Cambridge, 1997. [8] J. Cvitani´c and J. Zhang. The steepest descent method for forward-backward SDEs. Electron. J. Probab., 10:1468–1495 (electronic), 2005. [9] F. Delarue. On the existence and uniqueness of solutions to FBSDEs in a non-degenerate case. Stochastic Process. Appl., 99(2):209–286, 2002. [10] F. Delarue and S. Menozzi. A forward-backward stochastic algorithm for quasi-linear PDEs. Ann. Appl. Probab., 16(1):140–184, 2006. [11] J. Douglas, Jr., J. Ma, and P. Protter. Numerical methods for forward-backward stochastic differential equations. Ann. Appl. Probab., 6(3):940–968, 1996. [12] N. El Karoui and L. Mazliak, editors. Backward stochastic differential equations, volume 364 of Pitman Research Notes in Mathematics Series. Longman, Harlow, 1997. Papers from the study group held at the University of Paris VI, Paris, 1995–1996. [13] N. El Karoui, S. Peng, and M. C. Quenez. Backward stochastic differential equations in finance. Math. Finance, 7(1):1–71, 1997. [14] B. D Ewald and R. Temam. Numerical analysis of stochastic schemes in geophysics. SIAM J. NUMER. ANAL., 42(6):2257–2276, 2005. [15] I.M. Gradshteyn and I.M. Ryzhik. Table of Integral Series adn Products, volume Sixth Edition. Academic Press, 2000. [16] Ben-yu Guo, Jie Shen, and Cheng-long Xu. Spectral and pseudospectral approximations using Hermite functions: application to the Dirac equation. Adv. Comput. Math., 19(1-3):35–55, 2003. Challenges in computational mathematics (Pohang, 2001). [17] K. Itˆo. Multiple wiener integral. J.Math. Soc. Japan, 3:157–169, 1951. [18] P.E. Kloeden and E. Platen. Numerical Solution of Stochastic Differential Equations. SpringerVerlag, Berlin, 1992. [19] J. Ma, P. Protter, J. San Martin, and S. Torres. Numerical method for backward stochastic differential equations. Ann. Appl. Probab., 12(1):302–316, 2002.

24

[20] J. Ma, P. Protter, and J. Yong. Solving forward-backward stochastic differential equations explicitly—a four step scheme. Probab. Theory Related Fields, 98(3):339–359, 1994. [21] J. Ma and J. Yong. Forward-backward stochastic differential equations and their applications, volume 1702 of Lecture Notes in Mathematics. Springer-Verlag, Berlin, 1999. [22] J. Ma and J. Zhang. Representations and regularities for solutions to BSDEs with reflections. Stochastic Process. Appl., 115(4):539–569, 2005. [23] G. N. Milstein and M. V. Tretyakov. Stochastic Numerics for Mathematical Physics. SpringerVerlag, Berlin, 2004. [24] G. N. Milstein and M. V. Tretyakov. Numerical algorithms for forward-backward stochastic differential equations. SIAM J. Sci. Comput., 28(2):561–582, 2006. ´ Pardoux and S. Peng. Backward stochastic differential equations and quasilinear parabolic [25] E. partial differential equations. In Stochastic partial differential equations and their applications (Charlotte, NC, 1991), volume 176 of Lecture Notes in Control and Inform. Sci., pages 200–217. Springer, Berlin, 1992. [26] G´abor Szeg˝o. Orthogonal polynomials. American Mathematical Society, Providence, R.I., fourth edition, 1975. American Mathematical Society, Colloquium Publications, Vol. XXIII. [27] J. Zhang. A numerical scheme for BSDEs. Ann. Appl. Probab., 14(1):459–488, 2004. [28] Y. Zhang and W. Zheng. Discretizing a backward stochastic differential equation. Int. J. Math. Math. Sci., 32(2):103–116, 2002.

25