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 equationsLet the initial condition be an -measurable random vector in . An -adapted stochastic process is called a solution of equation (1) ifholds 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 thatand almost surely for all ,And these conditions imply that (4) and (5) are well defined.
One important property for the stochastic integral is thatand 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 ifwhere 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 equationwhere 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 termsto Euler scheme which gives the following schemewhere
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 arewhere 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 integralsfor 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) isHere 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, soOn the other hand, the double Wiener integralscould 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 processhas the Fourier serieswhere .
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 integralsIn formula (26), we haveIn 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)).
| ||||||||||||||||||||||||||||||||||||||||||
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 ().
| ||||||||||||||||||||||||||||||||||||||||||||||
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
Conflicts of Interest
The author declares that they have no conflicts of interest.