A regression-based Monte Carlo method to solve two-dimensional forward backward stochastic differential equations

The purpose of this paper is to investigate the numerical solutions to two-dimensional forward backward stochastic differential equations(FBSDEs). Based on the Fourier cos-cos transform, the approximations of conditional expectations and their errors are studied with conditional characteristic functions. A new numerical scheme is proposed by using the least-squares regression-based Monte Carlo method to solve the initial value of FBSDEs. Finally, a numerical experiment in European option pricing is implemented to test the efficiency and stability of this scheme.


Introduction
In this paper, we consider the numerical solutions to the two-dimensional decoupled forward backward stochastic differential equations (FBSDEs): where X t = (X 1 is how to discrete conditional expectation. Up to now, there have been lots of methods to solve this problem. Zhang et al. [2] constructed a kind of sparse-grid Gauss-Hermite quadrature rule and hierarchical sparse-grid interpolation to approximate this conditional expectation. Fu et al. [3] gave a method of spectral sparse grid approximations to deal with high-dimensional conditional expectation. It is known that the Fourier transform is an important tool in option pricing, not only in the SDE framework but also in the ODE framework. With the Fourier cos transform, one can convert conditional expectation to a series form. By constructing a function basis from truncating series, one can approximate the conditional expectation. More efficient methods and fast algorithms follow this idea. For details, one can refer to [4][5][6]. These schemes are applied to many fields, such as option pricing [7][8][9][10], portfolio optimization [11,12], and so on. Ruijter and Oosterlee [13] extended the Fourier cos method to two-dimensional FBSDE, named Fourier cos-cos method, and gave a numerical scheme to pricing European option and Bermudan option under GBM model and Heston stochastic volatility. Recently, Meng and Ding [14] investigated a Fourier sin-sin method named modified Fourier sin-sin method to price rainbow options within two-dimensional BSDE. Numerical experiments showed that its convergence and efficiency were expected. Inspired by the literature, we extend the idea to solving two-dimensional FBSDEs by using the Fourier cos-cos transform and the leastsquare Monte Carlo regression to obtain the numerical solution to FBSDEs (1) and (2). It is a supplement to our previous work in [15], and it can be extended to high-dimensional FBSDEs. The paper is organized as follows. In Sect. 2, some assumptions about FBSDEs are given to ensure the existence of solution. In the discretization scheme of forward equation (1), we use the classical Euler scheme which was used by Zhao et al. [16]. For backward equation (2), we use the theta scheme. In Sect. 3, we give the approximations and their error analysis of conditional expectations from the discretization of backward equation (2). In Sect. 4, we present a numerical scheme based on the least-squares Monte Carlo regression and provide an example in option pricing for a numerical experiment. In Sect. 5, we conclude our investigation.

Discretization of FBSDEs
In this section, we denote by L 2 T (R 2 ) the set of F T -measurable random variables X : → R 2 which are square integrable, and by H 2 T (R) the set of predictable processes η : where | · | is the standard Euclidean norm in the Euclidean space R. The terminal condition Y T in equation (2) is F T -measurable and square integrable. We give some assumptions: (A1) The function g(x) is uniformly global functional Lipschitz continuous. (A2) The functions μ(x) and σ (x) are uniformly Lipschitz continuous and satisfy a linear growth condition. (A3) The generator f (t, y, z) satisfies the following continuity condition: for any (t 2 , y 2 , z 2 ), (t 1 , y 1 , Assumptions (A1), (A2), and (A3) can guarantee the existence and uniqueness of solution (X t , Y t , Z t ) to FBSDEs (1)- (2). Now we are in the position to discretize FBSDEs (1) and (2) by using the Euler scheme. Given a partition : In the time interval [0, T], we rewrite BSDE (2) to the following form: Considering Y t to be an (F t )-adapted process, we take conditional expectations on both sides of equation (3) with respect to filtration F t m-1 , and then we have an iteration backward equation where Multiplying W * m and taking conditional expectations on both sides of (4), we have Applying the theta discretization method to (4),(5), we obtain a discrete solution (Y m-1 , Z m-1 ) to approximate the solution (Y t-1 , Z t-1 ) to BSDE (2): Here, θ 1 and θ 2 are two parameters in the theta discretization scheme. As a consequence of the Feynman-Kac theorem, the terminal values Y M and Z M are both deterministic func- where ∇ is a normal gradient operator with respect to the augment. Now, in combination with equations (6) and (7), we know that the solution (Y m-1 , Z m-1 ) is represented by the kinds of conditional expectations for some function υ(x). In these expectations, the first conditional expectation is onedimensional and the second is two-dimensional. Motivated by successful use of the Fourier cos-cos method in two-dimensional BSDEs, we use the Fourier transform to obtain the approximation expressions of the above conditional expectations.

Approximation of conditional expectation and error analysis
In this section, we give the approximation of conditional expectations U(x), V (x) and their error analysis. First, we give the approximation of U(x). Let p(y|x) denote the conditional density function of X m given by X m-1 = x. The symbol in theorems below means that the first term in the summation is weighted by one-half, and Re{·} denotes the real part of a complex number. 2 , the conditional expectation U(x) has the following expansion: where Proof For a truncated finite integration region D, we have By using the Fourier cos-cos transform to p(y 1 , y 2 |x 1 , x 2 ) in D, we have where A k 1 ,k 2 (x 1 , x 2 ) is the Fourier cosine coefficient of p(y 1 , y 2 |x 1 , x 2 ): Then we have With the cos formula 2 cos α cos β = cos(α + β) + cos(αβ), the integral in equation (9) will be changed to the following form: Substituting the above equation into (10), we can obtain the form of equation (8).
Next, we consider the two-dimensional conditional expectation The difficulty of V (x) is to deal with the Brownian motion W m . Note that the components of W m are independent, we can handle them separately. Denote and assume that the given condition is (X 1, m-1 , X 2, m-1 ) = (x 1 , x 2 ). Then, from the forward scheme for equation (1), we can revise it to another form defined by We find that υ(y 1 , y 2 )ρ x 2 (y 2 )p(y 1 , y 2 |x 1 , x 2 ) dy 1 dy 2 .
These integrals are similar to U(x) and can be calculated by using the method in Theorem 3.1. Next we directly give the expansion of V 1 (x) and V 2 (x).
Remark 1 There are some results to approximate conditional expectation, such as using polynomial basis functions, Malliavin approach, and Monte Carlo sequence convergence (see [17][18][19][20][21][22]). Most of them consider the time-spatial approximation. In fact, it needs much more time to implement, especially in dealing with high-dimensional conditional expectation. The results in Theorem 3.1 and Theorem 3.2 show that they can contain much information and have many advantages to deal with high-dimensional FBSDEs.
Theorem 3.1 and Theorem 3.2 give us some idea to approximate conditional expectations. For suitable integers N 1 , N 2 , the conditional expectations U(x), V (x) can be approximated by the truncation terms with the error with the error υ(y 1 , y 2 )ρ x j (y j )p(y 1 , y 2 |x 1 , x 2 ) dy 1 dy 2 for j = 1, 2. Now, we give the error analysis of the approximations. Ruijter and Oosterlee [13] pointed out that the coefficients A k 1 ,k 2 (x 1 , x 2 ) usually decay faster than B k 1 ,k 2 . Thus, we find that the error 2 (x) converges exponentially in N 1 and N 2 for density functions in the class C ∞ ([a 1 , b 1 ] × [a 2 , b 2 ]), i.e., for some positive constants P 2 , β, N , where β ≥ N, N = min{N 1 , N 2 }. On the other hand, according to [23], B k 1 ,k 2 exhibits at least algebraic convergence and gives us information of algebraic convergence of Fourier series, i.e., for suitable positive constants N, n, P, Q, we have +∞ After interchanging the summation and integration, we rewrite 3 (x) in another form: It then follows that From (11)-(15), with a properly chosen truncation of the integration range, the overall error (x) converges. With the same method, we can also prove that the overall error Therefore, if we choose a suitable region [a 1 , b 1 ] × [a 2 , b 2 ], then the errors of the approximation can be well controlled. We can use approximations (11) and (12) as a substitution to conditional expectations. The key is to choose basis functions. In the next section, we state our basis functions and employ the least-squares Monte Carlo regression method to numerical FBSDEs (1) and (2).

Numerical experiment
In this section, we give the numerical scheme to FBSDEs (1) and (2) based on the leastsquare Monte Carlo regression and perform a numerical experiment in pricing European option. In the following, we give the basis functions and corresponding coefficients α j,m at time t m . For approximations Y m , Z m , we use truncation functions to represent.
First, we state the numerical scheme. Under the condition of value (Y m , Z m ), we implement the following least-square regression by using finite-dimensional basis functions, respectively, to approximate Y m-1 , Z m-1 at each time t m : Notice that if θ = 1 in (16), the scheme will deduce to the representation in Gobet et al. [24]. Many numerical experiments show that the theta scheme is of second-order convergence when θ = 0.5. Following this idea, we also consider θ = 0.5. On the choice of basis functions, Gobet et al. use hypercubes and global polynomials as basis functions to test the effectiveness and stability under the assumption of assets following geometric Brownian motions (GBMs). Unfortunately, they only give one-dimensional numerical experiments to test the efficiency. We want to know the stability and efficiency in a high-dimensional space. From Sect. 3, we can use conditional characteristic functions to express the basis functions. Next, we simplify the basis functions by following Theorems 3.1 and 3.2. We assume that the underly assets follow GBM, i.e., μ j (x j ) = μ j x j and σ j (x j ) = σ j x j for j = 1, 2. Then the basis functions with respect to U(x) are given by m,k (x) = cos θ 1 · cos θ 2 · exp(β m,k ), and for V j (x) (j = 1, 2), the function bases are given j m,k (x) = cos θ 1 · cos θ 2 · exp(β m,k )/σ j x. Here, β m,k = -1 2 and ξ i j denotes the jth cumulant of the stochastic variable X i T . Denote by e Y = Y 0 -Y 0 the error of the numerical solution of Y . In the experiment, we present an application of our scheme to financial problems, i.e., pricing European option and hedging strategy. We consider option pricing of a basket call option in the Black-Scholes model. Someone has two kinds of assets. Denote by p t and X t = (X 1 t , X 2 t ) the bond price and the prices of two independent stocks, respectively, that satisfy dp t = rp t dt,