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

Authors: Yousef Alnafisah
Journal: Abstract and Applied Analysis
 (2018)

---

## Abstract

Multiple stochastic integrals of higher multiplicity cannot always be expressed in terms of simpler stochastic integrals, especially when the Wiener process is multidimensional. In this paper we describe how the Fourier series expansion of Wiener process can be used to simulate a two-dimensional stochastic differential equation (SDE) using Matlab program. Our numerical experiments use Matlab to show how our truncation of Itô’-Taylor expansion at an appropriate point produces Milstein method for the SDE.


---

## Body

## 1. 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–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 O(h(2/3-ϵ)) for a one-dimensional diffusion process where h is the step-size and then they generalize the result to SDEs of any dimension with O(hlog⁡1/h) 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.


## 2. Stochastic Differential Equations (SDEs)

Definition.Let {W(t)}t≥0 be a d-dimensional standard Brownian motion on a probability space (Ω,F,P) equipped with a filtration F=(Ft)t≥0, a=a(t,x) a d-dimensional vector function (calleddrift coefficient), and b=b(t,x) a d×d-matrix function (calleddiffusion coefficient).The stochastic processX=X(t) considered in this paper can be described bystochastic differential equations(1)dXt=at,Xtdt+bt,XtdWt,t∈0,T.Let the initial condition X(0)=x be an F0-measurable random vector in Rd. An Ft-adapted stochastic process X=(X(t))t≥0 is called a solution of equation (1) if(2)Xt=X0+∫0tas,Xsds+∫0tbs,XsdWsholds a.s.The conditions that the integral processes,(3)∫0tas,Xsds,∫0tbs,XsdWs,are well defined are required for (2) to hold and for the functions a(s,X(s)) and b(s,X(s)) we have the following conditions that(4)E∫0tb2s,Xsds<∞,and almost surely for all t≥0,(5)∫0tas,Xsds<∞.And these conditions imply that (4) and (5) are well defined.One important property for the stochastic integral is that(6)∫0tWsdWs=12∫0tdW2s-12∫0tds=12W2t-t2,and for more details on stochastic integral see [1].


### 2.1. Existence and Uniqueness Theorems

The following theorem, which will be stated without proof, gives sufficient conditions for existence and uniqueness of a solution of a stochastic differential equation.(i)
Measurability: leta: [0,∞)×Rd→Rd and b: [0,∞)×Rd→Rd×d are jointly Borel measurable in [t0,T]×Rd.(ii)
Lipschitz condition: there is a constantA>0 such that a(t,x)-a(t,y)≤Ax-y, and b(t,x)-b(t,y)≤Ax-y, for all t∈[t0,T] and x,y∈R.(iii)
Growth condition: there is a constantK>0 such that a(t,x)2≤K2(1+x2), and b(t,x)2≤K2(1+x2), for all t∈[t0,T] and x,y∈R.Theorem 1.
Under these conditions ((i)–(iii)) the stochastic differential equation (1) has a unique solution X(t)∈[t0,T] with (7)supt0≤t≤TEXt2<∞.Proof.
See Kloeden and Platen [1], Theorem 4.5.3.


### 2.2. 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,T], F is the σ-algebra of Borel sets, and P is the Wiener measure. We consider an approximate solution xh of (1) which uses a subdivision of the interval [0,T] into a finite number N of subintervals which we assume to be of length h=T/N. Also we assume the approximate solutions xh are random variables on Ω. Now we say that the discrete time approximation xh with the step-size h converges strongly of order γ at time T=Nh to the solution X(t) if(8)Exh-XTp≤Chγp,h∈0,1,where the strong convergence will be in Lp space and X(T) is the solution to the stochastic differential equation. C is a positive constant and C is independent of h.


## 2.1. Existence and Uniqueness Theorems

The following theorem, which will be stated without proof, gives sufficient conditions for existence and uniqueness of a solution of a stochastic differential equation.(i)
Measurability: leta: [0,∞)×Rd→Rd and b: [0,∞)×Rd→Rd×d are jointly Borel measurable in [t0,T]×Rd.(ii)
Lipschitz condition: there is a constantA>0 such that a(t,x)-a(t,y)≤Ax-y, and b(t,x)-b(t,y)≤Ax-y, for all t∈[t0,T] and x,y∈R.(iii)
Growth condition: there is a constantK>0 such that a(t,x)2≤K2(1+x2), and b(t,x)2≤K2(1+x2), for all t∈[t0,T] and x,y∈R.Theorem 1.
Under these conditions ((i)–(iii)) the stochastic differential equation (1) has a unique solution X(t)∈[t0,T] with (7)supt0≤t≤TEXt2<∞.Proof.
See Kloeden and Platen [1], Theorem 4.5.3.


## 2.2. 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,T], F is the σ-algebra of Borel sets, and P is the Wiener measure. We consider an approximate solution xh of (1) which uses a subdivision of the interval [0,T] into a finite number N of subintervals which we assume to be of length h=T/N. Also we assume the approximate solutions xh are random variables on Ω. Now we say that the discrete time approximation xh with the step-size h converges strongly of order γ at time T=Nh to the solution X(t) if(8)Exh-XTp≤Chγp,h∈0,1,where the strong convergence will be in Lp space and X(T) is the solution to the stochastic differential equation. C is a positive constant and C is independent of h.


## 3. 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 order1/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(9)dXit=ait,Xtdt+∑k=1dbikt,XtdWkt,Xi0=Xi0,where i=1,…,d on an interval [0,T], for a d-dimensional vector X(t), with a d-dimensional Brownian path W(t).In order to approximate the solution, we assume[0,T] is divided into N equal intervals of length h=T/N.


### 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.(10)xij+1=xij+aijh,xjh+∑k=1dbikjh,xjΔWkj,where ΔWk(j)=Wk((j+1)h)-Wk(jh) and our numerical approximation to X(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 the Milstein scheme by adding the quadratic terms(11)∑k,l=1dρikljh,xjAkljto Euler scheme which gives the following scheme(12)xij+1=xij+aijh,xjh+∑k=1dbikjh,xjΔWkj+∑k,l=1dρikljh,xjAklj,where(13)ΔWkj=Wkj+1h-Wkjh,Aklj=∫jhj+1hWkt-WkjhdWlt,ρiklt,x=∑m=1qbmkt,x∂bil∂xmt,x.The implementation of the Euler scheme is easy to do as it only needs to generate the normal distribution for the standard Brownian motionΔWk(j) but it is not easy to generate the integral Akl(j) 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 Akl(j) using the Fourier method when we have two-dimensional SDEs.Before the implementation of Milstein scheme we need to mention some facts about the two-level approximation.


## 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.(10)xij+1=xij+aijh,xjh+∑k=1dbikjh,xjΔWkj,where ΔWk(j)=Wk((j+1)h)-Wk(jh) and our numerical approximation to X(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 the Milstein scheme by adding the quadratic terms(11)∑k,l=1dρikljh,xjAkljto Euler scheme which gives the following scheme(12)xij+1=xij+aijh,xjh+∑k=1dbikjh,xjΔWkj+∑k,l=1dρikljh,xjAklj,where(13)ΔWkj=Wkj+1h-Wkjh,Aklj=∫jhj+1hWkt-WkjhdWlt,ρiklt,x=∑m=1qbmkt,x∂bil∂xmt,x.The implementation of the Euler scheme is easy to do as it only needs to generate the normal distribution for the standard Brownian motionΔWk(j) but it is not easy to generate the integral Akl(j) 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 Akl(j) using the Fourier method when we have two-dimensional SDEs.Before the implementation of Milstein scheme we need to mention some facts about the two-level approximation.


## 4. Two-Level Approximation

We need to generate the incrementsΔWk(j) 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.That is,(14)ΔWkr,j=ΔWkr+1,2j+ΔWkr+1,2j+1,where r∈N and ΔWk(r,j)=Wk((j+1)h(r))-Wk(jh(r)) with h(r)=T/2r.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 r by firstly generating the increments in the LHS ΔWk(r,j) and then conditionally generating the increments in the RHS. We do the same process for each level r+2, r+3 and so on.We will see from the following section that the extension of Milstein tod≥2 is not easy to do. However we could implement special class of equations for Milstein scheme using only ΔWk(j). This could be done from the observation that Akl(j)+Alk(j)=2Bkl(j) where Bkl(j)=1/2ΔWk(j)ΔWl(j) if k≠l and Bkk(j)=1/2{(ΔWk(j))2-h}.


### 4.1. Empirical Estimation of the Error of a Numerical Method

Usually we do not know the solutions of the stochastic differential equation explicitly; therefore we could not directly estimate the mean errorEXT-xh which is the absolute value of the difference between the approximation solution xh and the solution X(T) of an SDE (1). Assume the approximate solution xh converges to the solution X(T) as we decrease the step-size and go to zero. Then we can estimate the order of convergence for a particular scheme by repeating R different independent simulations of sample paths. We will use the following estimator ϵ=1/RExr-x^r for different approximation solutions x(r) and x^(r) for different range value of h. So for any numerical method if we have a bound for the error E|xh-xh/2|≤C1hγ then Exh/2-xh/4≤C1(h/2)γ and then Exh/4-xh/8≤C1(h/22)γ and so on. Therefore we will get a geometric series; then we will obtain(15)EXT-xh≤∑h=0∞C1h2kγ=C1hγ1-2-γ.So from (15) we could estimate the convergence and the constant.


### 4.2. 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(16)dX1t=X2tdW1t+X1t+tdW2t,dX2t=e-X22tdW1t+X1t-X2tdW2t,for0≤t≤1,withX10=2,X20=0,where W1 and W2 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(17)xij+1=xij+∑k=12bikxjΔWkj+∑k=12ρiklxjAklj.But the major difficulty here is that the double stochastic integrals(18)Aklj=∫jhj+1hWlt-WljhdWktfor k≠l 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 ρikl(x)=∑m∂bik/∂xmbml(x) for the SDEs (16).We have(19)ρ111=e-X22,ρ112=X1-X2ρ121=X2ρ211=-2e-2X22ρ221=X2-e-X22ρ222=t+X2ρ122=X1+tρ212=-2X2X1-X2e-X22.Then, the Milstein approximation for (16) is(20)x1j+1=x1j+x2jΔW1j+x1j+jhΔW2j+e-x2j2A11+x1j-x2jA12+x2jA21+x1j-jhA22x2j+1=x2j+e-x2j2ΔW1j+x1j-x2jΔW2j-2x2je-2x2j2A11+x2j-e-x2j2A21+x2j-jhA22-2x2je-x2j2x1j-x2jA12.Here in this approximation we have the double Wiener integrals A11, A12, A21, and A22. The double Wiener integrals A11 and A22 in (20) are easily computed from the Wiener increments ΔW1(j) and ΔW2(j), respectively, so(21)A11=∫jhj+1hW1t-W1jhdW1t=12ΔW1j2-h,A22=∫jhj+1hW2t-W2jhdW2t=12ΔW2j2-h.On the other hand, the double Wiener integrals(22)A12=∫jhj+1hW2t-W2jhdW1t,A21=∫jhj+1hW1t-W1jhdW2tcould 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(23)Wkt-thWkhfor0≤t≤hhas the Fourier series(24)Wkt-thWkh=12ak,0+∑r=1∞ak,rcos⁡2rtπh+bk,rsin⁡2rtπh,where k=1,…,d.Here the coefficientsak,r and bk,r are independent random variables with N(0,h/2r2π2) distributed and we could derive them from the Fourier integrals,(25)ak,r=2h∫0hWks-shWkhcos⁡2rtπhds,bk,r=2h∫0hWks-shWkhsin⁡2rtπhds.For eachk=1,…,d and r=1,…,p, when we integrate (24) over the interval [0,h], we will obtain the approximation of multiple Stratonovich integrals(26)Jk1,k2p=12hξk1ξk2-12hak2,0ξk1-ak1,0ξk2+hAk1,k2pk1,k2=1,…,d.In formula (26), we have(27)Ak1,k2p=12π∑r=1p1rξk1,rηk2,r-ηk1,rξk2,r;(28)ak,0=-1π2h∑r=1p1rξk,r-2hρpμk,p;whereρp=112-12π2∑r=1p1r2.In addition, ξk,ξk,r,ηk,r and μk,p are independent standard Gaussian random variables.For the truncation of Fourier series we require a convergence rate of orderh for the global error for the Milstein scheme and we will use (26) to express the double integral Akl(j) for k≠l. 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 Ch3 for some positive constant C. The algorithm of Kloeden et al. [18] has an MSE of order h2/p and then we obtain that Ch3=h2/p which gives h=1/Cp. Hence we want the number of terms in the truncated sum p to be proportional to h-1.We know from the symmetry relation that for any double integral we haveA12(j)+A21(j)=2B12(j) where B12(j)=1/2ΔW1(j)ΔW2(j).


## 4.1. Empirical Estimation of the Error of a Numerical Method

Usually we do not know the solutions of the stochastic differential equation explicitly; therefore we could not directly estimate the mean errorEXT-xh which is the absolute value of the difference between the approximation solution xh and the solution X(T) of an SDE (1). Assume the approximate solution xh converges to the solution X(T) as we decrease the step-size and go to zero. Then we can estimate the order of convergence for a particular scheme by repeating R different independent simulations of sample paths. We will use the following estimator ϵ=1/RExr-x^r for different approximation solutions x(r) and x^(r) for different range value of h. So for any numerical method if we have a bound for the error E|xh-xh/2|≤C1hγ then Exh/2-xh/4≤C1(h/2)γ and then Exh/4-xh/8≤C1(h/22)γ and so on. Therefore we will get a geometric series; then we will obtain(15)EXT-xh≤∑h=0∞C1h2kγ=C1hγ1-2-γ.So from (15) we could estimate the convergence and the constant.


## 4.2. 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(16)dX1t=X2tdW1t+X1t+tdW2t,dX2t=e-X22tdW1t+X1t-X2tdW2t,for0≤t≤1,withX10=2,X20=0,where W1 and W2 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(17)xij+1=xij+∑k=12bikxjΔWkj+∑k=12ρiklxjAklj.But the major difficulty here is that the double stochastic integrals(18)Aklj=∫jhj+1hWlt-WljhdWktfor k≠l 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 ρikl(x)=∑m∂bik/∂xmbml(x) for the SDEs (16).We have(19)ρ111=e-X22,ρ112=X1-X2ρ121=X2ρ211=-2e-2X22ρ221=X2-e-X22ρ222=t+X2ρ122=X1+tρ212=-2X2X1-X2e-X22.Then, the Milstein approximation for (16) is(20)x1j+1=x1j+x2jΔW1j+x1j+jhΔW2j+e-x2j2A11+x1j-x2jA12+x2jA21+x1j-jhA22x2j+1=x2j+e-x2j2ΔW1j+x1j-x2jΔW2j-2x2je-2x2j2A11+x2j-e-x2j2A21+x2j-jhA22-2x2je-x2j2x1j-x2jA12.Here in this approximation we have the double Wiener integrals A11, A12, A21, and A22. The double Wiener integrals A11 and A22 in (20) are easily computed from the Wiener increments ΔW1(j) and ΔW2(j), respectively, so(21)A11=∫jhj+1hW1t-W1jhdW1t=12ΔW1j2-h,A22=∫jhj+1hW2t-W2jhdW2t=12ΔW2j2-h.On the other hand, the double Wiener integrals(22)A12=∫jhj+1hW2t-W2jhdW1t,A21=∫jhj+1hW1t-W1jhdW2tcould 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(23)Wkt-thWkhfor0≤t≤hhas the Fourier series(24)Wkt-thWkh=12ak,0+∑r=1∞ak,rcos⁡2rtπh+bk,rsin⁡2rtπh,where k=1,…,d.Here the coefficientsak,r and bk,r are independent random variables with N(0,h/2r2π2) distributed and we could derive them from the Fourier integrals,(25)ak,r=2h∫0hWks-shWkhcos⁡2rtπhds,bk,r=2h∫0hWks-shWkhsin⁡2rtπhds.For eachk=1,…,d and r=1,…,p, when we integrate (24) over the interval [0,h], we will obtain the approximation of multiple Stratonovich integrals(26)Jk1,k2p=12hξk1ξk2-12hak2,0ξk1-ak1,0ξk2+hAk1,k2pk1,k2=1,…,d.In formula (26), we have(27)Ak1,k2p=12π∑r=1p1rξk1,rηk2,r-ηk1,rξk2,r;(28)ak,0=-1π2h∑r=1p1rξk,r-2hρpμk,p;whereρp=112-12π2∑r=1p1r2.In addition, ξk,ξk,r,ηk,r and μk,p are independent standard Gaussian random variables.For the truncation of Fourier series we require a convergence rate of orderh for the global error for the Milstein scheme and we will use (26) to express the double integral Akl(j) for k≠l. 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 Ch3 for some positive constant C. The algorithm of Kloeden et al. [18] has an MSE of order h2/p and then we obtain that Ch3=h2/p which gives h=1/Cp. Hence we want the number of terms in the truncated sum p to be proportional to h-1.We know from the symmetry relation that for any double integral we haveA12(j)+A21(j)=2B12(j) where B12(j)=1/2ΔW1(j)ΔW2(j).


## 5. Numerical Example

In the M-file in Listing1, I will approximate the value of the double integrals A12 and A21 and some explanations are shown for the formulas ((26)-(27)).Listing 1:Code to approximate the double integrals A12 and A21.
function [w,v,J1,J2] = F_year_J(N,h) N=5; % Here N=P for the truncation of the series in (26) T=1; h=T/N; s=sqrt(h); N1=randn; N2=randn; N7=randn; w=s∗N1; v=s∗N2; f=0; g1=0; g2=0; z=0; N3=randn(1,N); N4=randn(1,N); N5=randn(1,N); N6=randn(1,N); for n=1:N; f=f+(n.^ (−2)); c=(1/12)−((2∗(pi). ^ 2). ^ (−1))∗f; g1=g1+(n.^ (−1))∗N3(n); g2=g2+(n.^ (−1)).∗N4(n); z=z+(n.^ (−1)).∗(N3(n).∗N6(n)−(N5(n).∗N4(n))); A=1./(2∗pi)∗z; % to calculate formula (27) end B=(1/2)∗(w∗v); d1=(−1/pi)∗sqrt(2∗h)∗g1−(2∗sqrt(h∗c)∗N7); % to calculate formula (28) d2=(−1/pi)∗sqrt(2∗h)∗g2−(2∗sqrt(h∗c)∗N7); % to calculate formula (28) J1=(1/2)∗h∗N1∗N2−(1/2)∗s∗(d2∗N1−(d1∗N2))+h∗A; % to calculate formula (26) J2=2∗B−J1; endNow, after we represent the approximation of the double integralsA11,A22,A12 and A21, we could substitute them in the Milstein approximation in (20). After that we need to estimate the error for the Milstein solution in two-dimensional case and test the convergence order.The Matlab code in Listing2 calculates the Milstein error over the interval [0,1], with step-size (200, 400, 800, 1600, 3200) with a number of simulation (R=20000).Listing 2:Code for estimating the error for two-dimensional Milstein scheme.
T=1; NN=10; h=T/NN; R=10000; q=0; for r=1:R, x=2; y=0; xx=2; yy=0; for m=1:NN, hh=h/2; N=10; [wL,vL,J1L,J2L] = F_year_J(N,hh); [ w r , v r , J 1 r , J 2 r ] = F_year_J(N,hh); w=wL+wr; v=vL+vr; J1=J1L+J1r+wr∗vL; J2=J2L+J2r+wL∗vr; u=x+y∗wL+(x+(m−1)∗h)∗vL+exp(−y. ^ 2)∗(1/2)∗(wL. ^ 2−hh) +(x−y)∗J1L+y∗J2L+(x+(m−1)∗h)∗(1/2)∗(vL. ^ 2−hh); y=y+exp(−y.^ 2)∗wL+(x−y)∗vL−2∗y∗exp(−2∗y. ^ 2)∗(1/2)∗(wL. ^ 2−hh) +(y−exp(−y. ^ 2))∗J2L+(y+(m−1)∗h)∗(1/2)∗(vL. ^ 2−hh) −2∗y∗exp(−y. ^ 2)∗(x−y)∗J1L; x=u; u=x+y∗wr+(x+(m−1/2)∗h)∗vr+exp(−y. ^ 2)∗(1/2)∗(wr. ^ 2−hh) +(x−y)∗J1r+y∗J2r+(x+(m−1/2)∗h)∗(1/2)∗(vr. ^ 2−hh); y=y+exp(−y.^ 2)∗wr+(x−y)∗vr−2∗y∗exp(−2∗y. ^ 2)∗(1/2)∗(wr. ^ 2−hh) +(y−exp(−y. ^ 2))∗J2r+(y+(m−1/2)∗h)∗(1/2)∗(vr. ^ 2−hh) −2∗y∗exp(−y. ^ 2)∗(x−y)∗J1r; x=u; u=xx+yy∗w+(xx+(m−1)∗h)∗v+exp(−yy. ^ 2)∗(1/2)∗(w. ^ 2−h) +(xx−yy)∗J1+yy∗J2+(xx+(m−1)∗h)∗(1/2)∗(v. ^ 2−h); yy=yy+exp(−yy.^ 2)∗w+(xx−yy)∗v−2∗yy∗exp(−2∗yy. ^ 2)∗(1/2)∗(w. ^ 2−h) +(yy−exp(−yy. ^ 2))∗J2+(yy+(m−1)∗h)∗(1/2)∗(v. ^ 2−h) −2∗yy∗exp(−yy. ^ 2)∗(xx−yy)∗J1; xx=u; end q=q+abs(x−xx)+abs(y−yy); end (q/R)It is obvious from Table1 and the plotting in Figure 1 that the convergence seems to occur when we decrease the step-size and we obtain O(h) convergence. By estimating a range of values of h we could get the estimation of the convergence and also the estimation of the constant by using (15), so(29)Exh-XT≤C1-1/2h∀h.Table 1
The error results for the Milstein scheme in 2d case.
step-size error (ϵ) 0.0050 0.1318 0.0025 0.0673 0.00125 0.0347 0.00062 0.0177 0.00031 0.0088Figure 1
Milstein method for the two dimension SDEs.



---

*Source: 3805042-2018-05-10.xml*