Abstract and Applied Analysis
Volume 2018 (2018), Article ID 3805042, 7 pages
https://doi.org/10.1155/2018/3805042

Research Article

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

Yousef Alnafisah

Mathematics Department, College of Science, Qassim University, P.O. Box 6644, Buraydah 51452, Saudi Arabia

Correspondence should be addressed to Yousef Alnafisah; uan_0261@hotmail.com

Received 23 January 2018; Revised 29 March 2018; Accepted 5 April 2018; Published 10 May 2018

Academic Editor: Jaume Giné

Copyright © 2018 Yousef Alnafisah. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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.

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 for a one-dimensional diffusion process where is the step-size and then they generalize the result to SDEs of any dimension with 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 be a -dimensional standard Brownian motion on a probability space equipped with a filtration , a -dimensional vector function (called drift coefficient), and a -matrix function (called diffusion coefficient).

The stochastic process considered in this paper can be described by stochastic differential equations

Let the initial condition be an -measurable random vector in . An -adapted stochastic process 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 ,

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].

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: let : and : are jointly Borel measurable in .(ii)Lipschitz condition: there is a constant such that , and , for all and .(iii)Growth condition: there is a constant such that , and , for all and .

Theorem 1. Under these conditions ((i)–(iii)) the stochastic differential equation (1) has a unique solution with

Proof. See Kloeden and Platen [1], Theorem 4.5.3.

2.2. Strong Convergence for SDEs

Let be a probability space satisfying the following: is the set of continuous functions with the supremum metric on the interval , is the -algebra of Borel sets, and is the Wiener measure. We consider an approximate solution of (1) which uses a subdivision of the interval 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 .

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 order 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 on an interval , for a -dimensional vector , with a -dimensional Brownian path .

In order to approximate the solution, we assume 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 .

Firstly, consider the Euler-Maruyama approximation scheme.

where and our numerical approximation to will be denoted .

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

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 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 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,

where and with .

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 , and so on.

We will see from the following section that the extension of Milstein to is not easy to do. However we could implement special class of equations for Milstein scheme using only . This could be done from the observation that where if and .

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 error which is the absolute value of the difference between the approximation solution and the solution of an SDE (1). Assume the approximate solution converges to the solution as we decrease the step-size and go to zero. Then we can estimate the order of convergence for a particular scheme by repeating different independent simulations of sample paths. We will use the following estimator for different approximation solutions and for different range value of . So for any numerical method if we have a bound for the error then and then 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.

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

where and 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 in this approximation we have the double Wiener integrals , , , and . The double Wiener integrals and in (20) are easily computed from the Wiener increments and , respectively, so

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 .

Here the coefficients and are independent random variables with distributed and we could derive them from the Fourier integrals,

For each and , when we integrate (24) over the interval , we will obtain the approximation of multiple Stratonovich integrals

In formula (26), we have

In addition, and are independent standard Gaussian random variables.

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 for some positive constant . The algorithm of Kloeden et al. [18] has an MSE of order and then we obtain that which gives . Hence we want the number of terms in the truncated sum to be proportional to .

We know from the symmetry relation that for any double integral we have where .

5. Numerical Example

In the M-file in Listing 1, I will approximate the value of the double integrals and and some explanations are shown for the formulas ((26)-(27)).

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=sN1; v=sN2;

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./(2pi)z;

% to calculate formula (27)

end

B=(1/2)(wv);

d1=(−1/pi)sqrt(2h)g1−(2sqrt(hc)N7);

% to calculate formula (28)

d2=(−1/pi)sqrt(2h)g2−(2sqrt(hc)N7);

% to calculate formula (28)

J1=(1/2)hN1N2−(1/2)s(d2N1−(d1N2))+hA;

% to calculate formula (26)

J2=2B−J1;

end

Listing 1: Code to approximate the double integrals and .

Now, after we represent the approximation of the double integrals and , 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 Listing 2 calculates the Milstein error over the interval , with step-size (200, 400, 800, 1600, 3200) with a number of simulation ().

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);

= F_year_J(N,hh);

w=wL+wr; v=vL+vr; J1=J1L+J1r+wrvL; J2=J2L+J2r+wLvr;

u=x+ywL+(x+(m−1)h)vL+exp(−y. 2)(1/2)(wL. 2−hh)

+(x−y)J1L+yJ2L+(x+(m−1)h)(1/2)(vL. 2−hh);

y=y+exp(−y. 2)wL+(x−y)vL−2yexp(−2y. 2)(1/2)(wL. 2−hh)

+(y−exp(−y. 2))J2L+(y+(m−1)h)(1/2)(vL. 2−hh)

−2yexp(−y. 2)(x−y)J1L; x=u;

u=x+ywr+(x+(m−1/2)h)vr+exp(−y. 2)(1/2)(wr. 2−hh)

+(x−y)J1r+yJ2r+(x+(m−1/2)h)(1/2)(vr. 2−hh);

y=y+exp(−y. 2)wr+(x−y)vr−2yexp(−2y. 2)(1/2)(wr. 2−hh)

+(y−exp(−y. 2))J2r+(y+(m−1/2)h)(1/2)(vr. 2−hh)

−2yexp(−y. 2)(x−y)J1r; x=u;

u=xx+yyw+(xx+(m−1)h)v+exp(−yy. 2)(1/2)(w. 2−h)

+(xx−yy)J1+yyJ2+(xx+(m−1)h)(1/2)(v. 2−h);

yy=yy+exp(−yy. 2)w+(xx−yy)v−2yyexp(−2yy. 2)(1/2)(w. 2−h)

+(yy−exp(−yy. 2))J2+(yy+(m−1)h)(1/2)(v. 2−h)

−2yyexp(−yy. 2)(xx−yy)J1; xx=u; end

q=q+abs(x−xx)+abs(y−yy); end

(q/R)

Listing 2: Code for estimating the error for two-dimensional Milstein scheme.

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

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

step-size
error ()

0.0050
0.1318

0.0025
0.0673

0.00125
0.0347

0.00062
0.0177

0.00031
0.0088

Figure 1: Milstein method for the two dimension SDEs.

Conflicts of Interest

The author declares that they have no conflicts of interest.

References

P. E. Kloeden and E. Platen, Numerical solution of stochastic differential equations, vol. 23 of Applications of Mathematics (New York), Springer-Verlag, Berlin, 1992.

T. Rydén and M. Wiktorsson, “On the simulation of iterated Itô integrals,” Stochastic Processes and Their Applications, vol. 91, no. 1, pp. 151–168, 2001.

M. Wiktorsson, “Joint characteristic function and simultaneous simulation of iterated Itô integrals for multiple independent Brownian motions,” The Annals of Applied Probability, vol. 11, no. 2, pp. 470–487, 2001.

A. M. Davie, “Pathwise approximation of stochastic differential equations using coupling,” http://www.maths.ed.ac.uk/~adavie/coum.pdf.

S. Kanagawa, “The rate of convergence for approximate solutions of stochastic differential equations,” Tokyo Journal of Mathematics, vol. 12, no. 1, pp. 33–48, 1989.

S. T. Rachev and L. Ruschendorff, Mass Transportation Problems, Volume 1, Theory; Volume 2, Applications, Springer-Verlag, 1998.

J. Komlós, P. Major, and G. Tusnády, “An approximation of partial sums of independent RV's and the sample DF,” Zeitschrift für Wahrscheinlichkeitstheorie und Verwandte Gebiete, vol. 32, no. 1-2, pp. 111–131, 1975.

N. Fournier, “Simulation and approximation of Lévy-driven stochastic differential equations,” ESAIM: Probability and Statistics, vol. 15, pp. 233–248, 2011.

E. Rio, “Upper bounds for minimal distances in the central limit theorem,” Annales de l'Institut Henri Poincaré Probabilités et Statistiques, vol. 45, no. 3, pp. 802–817, 2009.

E. Rio, “Asymptotic constants for minimal distance in the central limit theorem,” Electronic Communications in Probability, vol. 16, pp. 96–103, 2011.

A. Alfonsi, B. Jourdain, and A. Kohatsu-Higa, “Pathwise optimal transport bounds between a one-dimensional diffusion and its Euler scheme,” The Annals of Applied Probability, vol. 24, no. 3, pp. 1049–1080, 2014.

A. Alfonsi, B. Jourdain, and A. Kohatsu-Higa, “Optimal transport bounds between the time-marginals of a multidimensional diffusion and its euler scheme,” Electronic Journal of Probability, vol. 20, article no. A070, 2015.

A. B. Cruzeiro, P. Malliavin, and A. Thalmaier, “Geometrization of Monte-Carlo numerical analysis of an elliptic operator: strong approximation,” Comptes Rendus Mathematique. Academie des Sciences, Paris, vol. 338, no. 6, pp. 481–486, 2004.

B. Charbonneau, Y. Svyrydov, and P. F. Tupper, “Weak convergence in the Prokhorov metric of methods for stochastic differential equations,” IMA Journal of Numerical Analysis (IMAJNA), vol. 30, no. 2, pp. 579–594, 2010.

L. N. Vasershtein, “Markov processes over denumerable products of spaces describing large system of automata,” Problems of Information Transmission, vol. 5, no. 3, pp. 64–72, 1969.

I. Gyöngy and N. Krylov, “Existence of strong solutions for Itô's stochastic equations via approximations,” Probability Theory and Related Fields, vol. 105, no. 2, pp. 143–158, 1996.

A. Davie, “KMT theory applied to approximations of SDE,” in Stochastic Analysis and Applications, vol. 100, pp. 185–201, Springer, 2014.

P. E. Kloeden, E. Platen, and I. Wright, “The approximation of multiple stochastic integrals,” Stochastic Analysis and Applications, vol. 10, no. 4, pp. 431–441, 1992.