The Implementation of Milstein Scheme in Two-Dimensional SDEs Using the Fourier Method

and Applied Analysis 3 called Euler-Maruyama scheme) which utilizes only the first two terms of the Taylor expansion and it attains the strong convergence γ = 1/2. Firstly, consider the Euler-Maruyama approximation scheme. x(j+1) i = x(j) i + ai (jh, x(j)) h + d ∑ k=1 bik (jh, x(j)) ΔW(j) k , (10) where ΔW(j) k = Wk((j + 1)h) − Wk(jh) and our numerical approximation toX(jh) will be denoted x(j). 3.2. The Milstein Scheme. We shall now introduce the Milstein scheme which gives an order one strong Taylor scheme. We could obtain theMilstein scheme by adding the quadratic terms


Introduction
Numerical analysis for stochastic differential equation (SDE) has seen a considerable development in recent years.There are many numerical methods for solving SDEs.Kloeden and Platen [1] described a method based on the stochastic Taylor series expansion but the major difficulty with this approach is that the double stochastic integrals cannot be easily expressed in terms of simpler stochastic integrals when the Wiener process is multidimensional.In the multidimensional case, the Fourier series expansion of Wiener process has been used to represent the double integrals by [1][2][3] but it needs to generate many random variables at each time.Therefore, it takes a lot of time to compute and also it is hard to extend to higher order.
There have been many studies for the numerical resolution of Stochastic differential equations.Davie [4] uses coupling and gives order one for the strong convergence for stochastic differential equations (SDEs).Kanagawa [5] investigates the rate of convergence in terms of two probability metrics between approximate solutions with i.i.d random variables.Rachev and Ruschendorff [6] developed Kanagawa's method by using the Komlós et al. theorem in [7].Fournier [8] uses the quadratic Vaserstein distance for the approximation of the Euler scheme and the results of Rio [9] which gives a very precise rate of convergence for the central limit theorem in Vaserstein distance.Also, Rio [10] provided precise bound estimates.Under uniform ellipticity, Alfonsi et al. [11,12] studied the Vaserstein bound for Euler method and they proved an (ℎ (2/3−) ) for a one-dimensional diffusion process where ℎ is the step-size and then they generalize the result to SDEs of any dimension with (ℎ√log(1/ℎ)) bound when the coefficients are time-homogeneous.Cruzeiro et al. [13] get an order one method and under the nondegeneracy they construct a modified Milstein scheme which obtains an order one for the strong approximation.Charbonneau et al. [14] investigate the Vaserstein bound [15] by using the weak convergence and Strassen-Dudley theorem.Convergence of an approximation to a strong solution on a given probability space was established by Gyöngy and Krylov in [16] using coupling.Davie in [17] applied the Vaserstein bound to solutions of vector SDEs and uses the Komlós et al. theorem to get order one approximation under a nondegeneracy assumption.The rest of this paper is organized as follows.Section 2 reviews some results concerning SDE.Section 3 presents some existing schemes for numerical resolution of SDE.In Section 4 we show the theoretical and implementation of Milstein scheme using the Fourier method.In the last section we give numerical example to the show the convergence behaviour.
The stochastic process  = () considered in this paper can be described by stochastic differential equations Let the initial condition (0) =  be an F 0 -measurable random vector in R  .An F  -adapted stochastic process  = (()) ≥0 is called a solution of equation (1) if holds a.s.The conditions that the integral processes, are well defined are required for (2) to hold and for the functions (, ()) and (, ()) we have the following conditions that and almost surely for all  ≥ 0, And these conditions imply that (4) and ( 5) are well defined.One important property for the stochastic integral is that and for more details on stochastic integral see [1].

Strong Convergence for SDEs.
Let (Ω, F, P) be a probability space satisfying the following: Ω is the set of continuous functions with the supremum metric on the interval [0, ], F is the -algebra of Borel sets, and P is the Wiener measure.
We consider an approximate solution  ℎ of (1) which uses a subdivision of the interval [0, ] into a finite number  of subintervals which we assume to be of length ℎ = /.Also we assume the approximate solutions  ℎ are random variables on Ω.Now we say that the discrete time approximation  ℎ with the step-size ℎ converges strongly of order  at time  = ℎ to the solution () if where the strong convergence will be in   space and () is the solution to the stochastic differential equation. is a positive constant and  is independent of ℎ.

Numerical Method for Approximating the SDEs
There are many numerical methods for solving stochastic differential equation; here we will mention two important schemes.The first one is the Euler-Maruyama scheme which will give strong order 1/2 and the second one is the Milstein scheme which has an order one for the strong convergence.
We will illustrate by a numerical example their convergence behaviour of Milstein scheme.Suppose we have the stochastic differential equation where  = 1, . . .,  on an interval [0, ], for a -dimensional vector (), with a -dimensional Brownian path ().
In order to approximate the solution, we assume [0, ] is divided into  equal intervals of length ℎ = /.

3.1.
Euler-Maruyama Scheme.The simplest numerical method for approximating the solution of stochastic differential equations is the stochastic Euler scheme (also called Euler-Maruyama scheme) which utilizes only the first two terms of the Taylor expansion and it attains the strong convergence  = 1/2.
Firstly, consider the Euler-Maruyama approximation scheme.

The Milstein Scheme.
We shall now introduce the Milstein scheme which gives an order one strong Taylor scheme.We could obtain the Milstein scheme by adding the quadratic terms (ℎ,  () )  ()  (11) to Euler scheme which gives the following scheme where The implementation of the Euler scheme is easy to do as it only needs to generate the normal distribution for the standard Brownian motion Δ ()  but it is not easy to generate the integral  ()  for the Milstein scheme when we have two or more dimensional SDEs.We will show by a numerical example in the next section how we could generate the integral  ()  using the Fourier method when we have twodimensional SDEs.
Before the implementation of Milstein scheme we need to mention some facts about the two-level approximation.

Two-Level Approximation
We need to generate the increments Δ ()  when we approximate the solution to (1) by using Milstein or other schemes; therefore Levy's construction of the Brownian motion will be used to simulate a sequence of approximations which converge to the solution.
We will call the two-level approximation in ( 14) the trivial coupling.We could generate the normal distribution in (14) for the increments for a given level  by firstly generating the increments in the LHS Δ (,)  and then conditionally generating the increments in the RHS.We do the same process for each level  + 2,  + 3 and so on.
We will see from the following section that the extension of Milstein to  ≥ 2 is not easy to do.However we could implement special class of equations for Milstein scheme using only Δ  )  and so on.Therefore we will get a geometric series; then we will obtain So from (15) we could estimate the convergence and the constant.

Two-Dimensional Stochastic Differential Equation.
In this section, we consider the two-dimensional stochastic differential equations and we need to test the convergence by using Milstein scheme.The SDEs that we will choose to implement our methods on are where  1 and  2 are independent standard Brownian motions.
For the two-dimensional SDEs ( 16), we could simply implement the Euler method by only generating some normal distributions.Now, we need to apply Milstein method to (16) and show the convergence between the final solutions of these methods.We need to find an approximation for the Milstein scheme for two-dimensional SDE.
For the SDEs ( 16), we have the Milstein scheme But the major difficulty here is that the double stochastic integrals for  ̸ =  cannot be so easily expressed in terms of simpler stochastic integral when the Wiener process is multidimensional.Therefore we will use the Fourier series expansion of Wiener process to represent the double integrals.
Before explaining the Fourier method let us start by applying the Milstein scheme ( 17) to ( 16) and then explain which terms that the Fourier method will be represented to.Before writing the Milstein approximation, we need to find the derivative terms   () = ∑  (  /  )  () for the SDEs (16).
We have Then, the Milstein approximation for ( 16) is Here On the other hand, the double Wiener integrals could not be expressed in terms of simpler stochastic integrals when the Wiener process is multidimensional.Therefore, for these integrals the Fourier series expansion will be used to approximate them.Now we will explain the idea of Fourier method as described in Kloeden, Platen [1,18].The Brownian bridge process has the Fourier series where  = 1, . . ., .
For the truncation of Fourier series we require a convergence rate of order ℎ for the global error for the Milstein scheme and we will use (26) to express the double integral  ()  for  ̸ = .So in order to have this convergence rate we need to compare the mean square error (MSE) of the approximation of the iterated Itô integrals to the discretization error of the Milstein scheme.As described in Kloeden and Platen [1], Corollary 10.6.5, and equation 10.6.16 we require an MSE bounded by ℎ 3 for some positive constant .The algorithm of Kloeden et al. [18] has an MSE of order ℎ 2 / and then we obtain that ℎ 3 = ℎ 2 / which gives ℎ = 1/.Hence we want the number of terms in the truncated sum  to be proportional to ℎ −1 .
We know from the symmetry relation that for any double integral we have

Numerical Example
In the M-file in Listing 1, I will approximate the value of the double integrals  12 and  21 and some explanations are shown for the formulas ((26)-( 27)).Now, after we represent the approximation of the double integrals  11 ,  22 ,  12 and  21 , we could substitute them in the Milstein approximation in (20).After that we need to estimate the error for the Milstein solution in twodimensional case and test the convergence order.
It is obvious from Table 1 and the plotting in Figure 1 that the convergence seems to occur when we decrease the step-size and we obtain (ℎ) convergence.By estimating a range of values of ℎ we could get the estimation of the convergence and also the estimation of the constant by using (15), so

Figure 1 :
Figure 1: Milstein method for the two dimension SDEs.

Table 1 :
The error results for the Milstein scheme in 2 case.