NUMERICAL PERFORMANCE OF SOME WONG-ZAKAI TYPE APPROXIMATIONS FOR STOCHASTIC DIFFERENTIAL EQUATIONS

We propose to use a Wong-Zakai type approximation method for the numerical evaluation of the solution of stochastic differential equations. The main feature of the method is that it takes advantage of well-developed techniques for solutions of ordinary differential equations and that its application to multidimensional problems is straightforward. We give focus to the evaluation of the numerical performance of the method. AMS Subject Classification: 65C30


Introduction
Ordinary differential equations (ODEs) are a common tool for modeling physical and social systems.Such models represent idealized versions of real systems, as they are purely deterministic.Stochastic differential equations (SDEs) are the instrument for building more realistic models, as they include random elements.Many areas of applications use SDEs including finance, economics, insurance, signal processing and filtering, several fields of biology and physics, population dynamics and genetics.
Since there are very few SDEs with exact analytical solutions, numerical techniques have to be used in the solution of SDEs.Unlike the solution of ODEs, SDEs can be approximated in two senses: strong and weak.We refer to strong approximations when we want to approximate the trajectory of the solution.Whenever we do not require the whole trajectory, but a functional from the solution, for instance, a moment, we talk about weak approximations.[2], [8], and [12] provide some introductions to numerical methods for SDEs.
A general Stratonovich SDE has the form where b : R n × [t 0 , T ] → R n is a d-dimensional vector valued function, σ : R n × [t 0 , T ] → R n×d is n × d-matrix-valued function that is twice continuously differentiable in the spatial variable and continuously differentiable in the time variable and W s is a d-dimensional Brownian motion.Commonly, b and σ are called the drift and the diffusion coefficient, respectively.It is the goal of this paper to illustrate how approximations of the Wong-Zakai type can are used as a numerical method for the approximations of a strong solution of SDEs.The main feature of the proposed approximations is that it uses the well-developed techniques of solution of ODEs.We point out how popular approximations techniques like the Euler approximations and the Milstein approximations can be seen as coarse versions of the Wong-Zakai approximations (see Section 3).
We also show that for some particular cases the Wong-Zakai approximation delivers the exact solution of the SDE; indeed we show that for those processes that can be written as functions of a Brownian Motion, the Wong-Zakai approximation delivers an exact solution (see Section 4).
The rate of convergence of these Wong Zakai approximations only have been studying for convergence in L p for p ≥ 2. (See [7] and [1]); with the help of Lyapunov inequality it is possible to obtain coarse rates of strong convergence in L 1 of order p for any p < 1/2 (see Theorem 1).However, in this paper we try to illustrate through numerical experiments the high-quality performance of these numerical approximations.We believe that the computer experiments presented make a strong case to study further properties of these Wong-Zakai approximations and more elaborated rates of convergence for the algorithm.
The remainder of this paper is organized as follows.In Section 2 we describe the proposed numerical technique.In Section 3 we show that the Milstein scheme can be seen as a coarse approximation of the proposed method.In Section 4 we derive a sufficient condition for the method to be exact.In Section 5 we prove that the rate of convergence of the algorithm is at least 1/4.In Section 6 some implementation issues are discussed.Finally, in Section 7 some numerical results are presented.

The Solution Method
We consider the class of Stratonovich stochastic differential equations where the diffusion coefficient is a function of the state variables.That is For this type of equations a Wong-Zakai type numerical solution method is proposed.[15] introduced Wong-Zakai type approximations.
Assume that we wish to approximate the solution to the above equations at points 0 = t 0 < t 1 < t 2 < • • • < t k = T of the interval [0, T ].Let Xj be the numerical approximation of X t j .For X0 = X 0 = x 0 , and for each sub-interval [t j , t j+1 ] ,j = 0, 1, . . ., k − 1, Xj+1 shall be calculated as the solution at time t j+1 to the following ordinary differential equation initial value problem: where ∆ j = t j+1 − t j and ∆W j = W t j+1 − W t j are the discrete time approximations of ds and dW s , respectively.Equation (3) may be rewritten in integral form:

Reduction to the Milstein Scheme
We shall now make some remarks on the 1-dimensional case n = d = 1.The Milstein scheme is a simple discrete time approximation derived from stochastic Taylor expansions.In the 1-dimensional case this scheme has the form where σ ′ (x, t) = ∂σ ∂x , and X M 0 = X 0 = x 0 .Under the assumption that b is once and σ twice continuously differentiable it can be shown that the Milstein scheme has the order of strong convergence γ = 1.0.A detailed explanation of the Milstein scheme is found in [8].We shall see that that scheme given in (4) becomes the Milstein scheme if some particular methods are used to approximate the Riemann-Stieltjes integrals where X(s), for t j ≤ s ≤ t j+1 is the solution of the ordinary differential equation (4) subject to the initial condition X(t j ) = X M j .The first integral may be approximated using the Euler method, which gives In the other hand, the second integral may be approximated using the 2nd order truncated Taylor method, which uses the approximation Replacing ( 6) and ( 7) in ( 4) we obtain which is the Milstein scheme for σ(x, t) ≡ σ(x)

A Particular Case where the Solution is Exact
We shall find a sufficient condition on b and σ, for which the proposed numerical scheme yields the exact solution to (8), given that we can solve analytically the ODE problem given by (3).We consider the class of scalar stochastic differential equation given by Suppose that Since Stratonovich integrals follow the classical chain rule formula, we conclude that Since X t is also the solution to (8), we see that and We obtain from ( 9) and ( 10) the identities and where Notice that for this equation to be satisfied the functions b and σ must be such that This implies that b(x) σ(x) = α, where α is an arbitrary constant.Thus We shall now see that ( 13) is a sufficient condition for the numerical method to produce the exact solution to (8).We shall then prove that the numerical method for the class of Stratonovich differential equations given by with initial value X 0 at t = 0, is exact.In [8] it is proved that the general solution to ( 14) is where h is given by For this particular type of stochastic differential equations, approximation (3) becomes We use the principle of separation of variables to solve this equation at t: Since X0 = X 0 we obtain which is precisely (15), the solution to (14).
We have proved that condition ( 13) is sufficient for the proposed numerical method to produce the exact solution to (8).Geometric Brownian motions are an important type of processes which satisfies this condition.

Rate of Convergence
It is possible under some conditions to obtain the rate of convergence of the Wong-Zakai algorithm for solutions of SDE's.In fact, the following result says that the strong rate of convergence of Wong-Zakai algorithm is at least of order 1/4.All the computer experiments made in this paper suggest a better rate of convergence.b denoted the class of bounded twice continuously differentiable functions and bounded one continuously differentiable functions respectively).Moreover, assume that ∂σ i,j /∂x i are bounded.Then for any where Xs the Wong-Zakai approximation described in section 2 for a time discretization of step size at most δ, and c p is a positive constant.
Proof.By the Lyapunov Inequality, and where M q is an appropriate constant.The result is a consequence of this last inequality.

Implementation Issues
We shall now focus on the computer implementation of the described numerical method.First, we rewrite equation ( 3) into a form that may be easier to implement.Recall that the increments W t − W s of a d-dimensional Brownian motion are N (0, (t − s)I) Gaussian distributed for 0 ≤ s < t, where I is the d-dimensional identity matrix.Thus, ∆W j is N (0, ∆ j I) Gaussian distributed.
Then in (3), we can replace ∆W j with ∆ j Z, to obtain where Z is a Gaussian distributed random variable with mean 0 and covariance the d-dimensional identity matrix.From equation ( 19), it becomes clear that the implementation of the numerical technique requires the generation of Gaussian distributed random variables and the solution of ODEs.Particularly, it is necessary to solve ODE initial value problems.
When solving ODE problems almost always, numerical techniques must be used since the available analytical techniques are not powerful enough to solve any ODE problem except the simplest.However, numerical methods for ODEs are rather more developed and available than numerical methods for SDEs.ODE problems are divided into two categories, stiff and non-stiff problems.[13] propose a definition of stiffness.A stiff problem is harder to solve than a non-stiff one.Therefore, the numerical methods for solving stiff equations are different from those for solving non-stiff equations.The methods used for solving non-stiff ODEs are based on the Runge-Kutta, Adams, or extrapolation methods.These methods are usually not suitable for solving stiff problems.The interested reader may refer to [6] for a survey on numerical methods for ODEs.
The hard work of the proposed method relies on the numerical method used to solve the ODE initial value problems stated in (19).Therefore, the accuracy of the described numerical method depends on the accuracy of the ODE initial value problem solver.For instance, if the Euler method is used as the ODE solver, the proposed method becomes the well known Euler-Maruyama method, and has thus strong order of convergence γ = 0.5.
For our numerical experiments, we used CVODE as ODE solver.CVODE is a solver, written in C, for stiff and non-stiff initial value problems for systems of ordinary differential equations.The underlying integration methods used in CVODE are variable-coefficient forms of the Adams and BDF (Backward Differentiation Formula) methods, for non-stiff and stiff problems, respectively.The interested reader may refer to [4] for a detailed presentation of CVODE and may download it from http://www.netlib.org/ode/cvode.tar.gz.

Numerical Performance
In this section, we will examine the numerical performance of the proposed numerical technique.Thus, we report numerical results for six test equations.We shall consider a time discretization of the interval [0, T ] given by where τ i = i * T /L, for i = 1, • • • , L and the step size h = T /L.
For the first three equations, we shall focus on the accuracy of the algorithm.To measure the accuracy of the algorithm we use the mean of the absolute error ǫ(T ).For this, we carried out a numerical simulation to get approximated for the mean of the absolute error ǫ(T ).For the computation of these quantities, we arrange the simulations into M batches of N simulations, and we estimate ∆ǫ(T ) in the following way.We denote Xi,j T the value produced by the Wong-Zakai algorithm in the k generated trajectory in the j batch at the time T , where the discretization of the time interval [0, T ] is the discretization discussed above, and we denote as X i,j T the corresponding value of the Itô process.We estimate the average error for each batch as Next, we estimate the mean of the batch averages as and therefore where t 1−α,M −1 is determined from the t Student t-distribution with M − 1 degrees of freedom.In this paper we shall use 20 batches, each with 100 trajectories, 1 − α = 0.9 (90% confidence intervals).Details of these computations can be found in [9][p.131].
On the other hand, for the last three test equations we shall mainly be concerned with the qualitative behavior of the simulated paths.
Test equation one: This test equation is a one-dimensional nonlinear problem, taken from [14].The nonlinear Stratonovich equation is given by, with α = −1 and β = 1.The exact solution of this equation is Notice that (20) is a particular case of (14) and thus the numerical scheme should produce the exact solution; the latter is in contrast with the error reported in [14], where Runge-Kutta type methods were used, where the mean absolute error is of order 0.01 for step size h = 2 −10 .

Test equation two:
The third test equation, taken from [10], is a linear two-dimensional system, whose Stratonovich form is with The exact solution of this equation is where and, Table 1 gives means of the absolute error for each of the two variables with step-size h = 2, σ ∈ [1, 10] and ρ = 1.The mean of the absolute errors is of a very small order, which suggests that the solution method returns the exact solution of the equation.The latter fact is in contrast with the numerical instabilities produced by the Euler approximation, as reported in [10].
Test equation 3: This test equation is a 2-dimensional linear SDE system, presented in [3], whose Stratonovich form is given by where U and V are the matrices defined by The exact solution of this equation is where p ± (t) = (−u − 1 2 v 2 ± u)t + vW t and P = P −1 .The stiffness of this system increases quadratically in terms of v. Table 2 gives the mean of the absolute error M of the numerical solution of (22) with u = 5, v ranging from 1 to 5 and step-size h = 1.As v increases, the error magnitude does not increase, and the stability of the solutions is maintained.The latter fact is in contrast with the poor stability properties of some Runge-Kutta type methods reported in [3] for v ≥ 3. The great accuracy obtained suggests that the method converges exactly for this test equation.

Test equation four:
This test equation corresponds to the interest rate model of Cox-Ingersoll-Ross (CIR) [5] for stochastic interest rates, whose Stratonovich form is where a ≥ 0, b ∈ R and σ > 0 are parameters.It is known that the solutions are non-negative.The proposed numerical scheme produces trajectories which take negative values (see Figure 1).A non-negative preserving numerical algorithm for the solution of ( 23) is presented in [11].Test equation five: The sixth example is the phased locked loop (PLL), from the filtering theory [2].The PLL is an important FM demodulator, whose model is The deterministic model has stable equilibrium at y 1 = 0, y 2 = 2πn(n = 0, ±1, . . .).Any trajectory that begins in the domain of attraction of a stable equilibrium will not leave this domain.However, in the stochastic model, small perturbations can cause a crossing of this domain.This phenomenon is known as slipping cycles.For small c slipping cycles are infrequent, but for big c they become more frequent.Figure 2 shows two trajectories produced by the pro-posed method, when the step-size h = 0.01.The numerical solution conserves the qualitative properties of the exact solution.
It is important to study the relationship between µ and b.For example if b ≫ µ the paths move back and forth between ± √ µt before settling down on one branch [2]. Figure 3 plots the numerical solution obtained with step-size h = 0.01, µ = 0.03 and b = 0.25, 0.75.The solutions follow the previously describe characteristics.

Conclusions
In this paper, we presented a Wong-Zakai type numerical method for the approximation of stochastic differential equations.The main feature of this method is that it uses the rather well developed and available numerical techniques for the simulation of ordinary differential equations.Moreover, contrary to standard approximation techniques, its application to multidimensional problems is straightforward.This allows one to implement it in any computer package that provides an ordinary differential equation solver.
We have also found a particular class of SDE for which the method produces the exact solution.In addition, we also showed that the method proposed delivers the exact solution when the solution is a function for of the Brownian Motion.We tested the performance of the method with several 1-dimensional and 2-dimensional equations.The results from the test problems suggested that the method produces very accurate approximations and that it often keeps the qualitative behavior of the solution paths.
We show that the Wong-Zakai algorithm has a coarse, strong converge rate of order p for any p < 1/2.It remains to be investigated a better strong convergence rate of the method and the extension of the method for general SDE where the diffusion term is a function of state and time.Other problems to address include adaptive versions of the Wong-Zakai Algorithm, studies of properties that are preserved when it is used this algorithm to approximate the solution of SDE's, and the study of the algorithm for explosive processes.

2 b and C 1 b
18) be a Stratonovich Stochastic differential equation, where we assume that b = (b i ) and σ = (σ i,j ) are bounded Lipschitz continuous functions, of class C (where C 2 b and C 1

Figure 2 :
Figure 2: Phase plane of the PLL model with c = 0.001(left) and c = 10 (right)

Table 1 :
Mean of the absolute error of (21) with ρ = 1 and T = 2, step size=2.M 1 is for Y 1 t and M 2 for Y 2

Table 2 :
Error and convergence rate of (22) with u = 5 and T = 1, step-size h = 1, M 1 is for X 1 t and M 2 for X 2