Weak Local Linear Discretizations for Stochastic Differential Equations: Convergence and Numerical Schemes F.Carbonell∗, J.C.Jimenez and R.Biscay

November 17, 2005

Abstract Weak Local Linear (WLL) discretizations are playing an increasing role in the construction of effective numerical integrators and inference methods for Stochastic Differential Equations (SDEs) with additive noise. However, due to limitations in the existing numerical implementations of WLL discretizations, the resulting integrators and inference methods have either been restricted to particular classes of autonomous SDEs or showed low computational efficiency. Another limitation is the absence of a systematic theoretical study of the rate of convergence of the WLL discretizations and numerical integratos. This task is the main purpose of the present paper. A second goal is introducing a new WLL scheme that overcomes the numerical limitations mentioned above. Additionally, a comparative analysis between the new WLL scheme and some conventional weak integrators is also presented.

AMC (2000): 60H10, 60H35, 65C30 Key words: Local Linearization, Weak Schemes, Weak Convergence, Stochastic Differential Equations

1

Introduction

In a number of problems in mathematical physics, biology, economics, finance and other fields the evaluation of Wiener functional space integrals and the estimation of diffusion processes play a prominent role. Examples are the computation of the Feynman-Kac formula in the analysis of wave scattering in random media [3], the computation of Lyapunov exponent of random dynamical systems [8], the estimation of the investment distribution that maximizes the expected utility function in an optimal portfolio problem [22], the estimation of continuous time stochastic volatilities models for stock prices [6], etc. In the solution of this kind of problems, weak numerical integrators for Stochastic Differential Equations (SDEs) have become an essential tool [31, 32, 15, 19]. During the last years, much progress has been made in the study of weak numerical integrators (see [21] for an updated review of these methods). Well-known are, for instance, the Euler, the Milstein, the Talay-Tubaro extrapolation, the Runge-Kutta and the Local Linearization (LL) methods. This paper focus on the class of weak LL schemes. The main idea of the LL method is obtaining a discretization of the d-dimensional nonlinear SDE dx(t) = f (t, x(t))dt + G(t)dw(t), t ∈ [t0 , T ]

(1)

x(t0 ) = x0 ,

by means of a local (piecewise) linear approximation of the drift coefficient f . Here, w is a m-dimensional standard Wiener process and G is a d × m matrix function. Two different ways of linearization have ∗ Instituto de Cibernética, Matemática y Física, Departamento de Sistemas Adaptativos, Calle 15, No. 551, e/ C y D, Vedado, La Habana 4, C.P. 10400, Cuba

1

been considered, namely, either based on truncated Taylor or Ito-Taylor expansions of the function f [18, 2, 24, 13]. They lead to two different weak LL discretizations that can be expressed by the recurrent relation (2) ytn+1 = ytn + φβ (tn , ytn ; tn+1 − tn ) + η(tn , ytn ; tn+1 − tn ), evaluated at an increasing sequence of discrete times tn ∈ [t0 , T ], where yt0 = x0 , φβ (t, y; δ) =

Rδ

e

∂f (t,y) (δ−s) ∂y

(f (t, y) + bβ (t, y)s)ds,

(3)

0

with bβ (t, y) = for all (t, y) ∈

R × Rd

⎧ ⎪ ⎨ ⎪ ⎩

∂f (t,y) ∂t , ∂f (t,y) ∂t

+

1 2

d P

2

⎫ ⎪ β=1 ⎬

(G(t)G| (t))k,l ∂∂yfk(t,y) , β=2 ⎪ ⎭ ∂yl

k,l=1

,

and δ > 0; and η(t, y; δ) is a zero mean stochastic process with variance matrix Σ(t, y; δ) =

Rδ

e

∂f (t,y) (δ−s) ∂y

G(t + s)G| (t + s)e(

∂f (t,y) | ) (δ−s) ∂y

ds.

(4)

0

The discretization with β = 1 corresponds to the truncated Taylor expansion, while other one corresponds to the truncated Ito-Taylor expansion. On the basis of these two discretizations, various weak LL schemes have been proposed [18, 24, 16], which differ in respect to the way of computing the integrals (3) and (4). In particular, the schemes introduced in [18], [24] have played a prominent role in the construction of effective inference methods for SDEs [23], [25], [5], [27] in the estimation of distribution functions in Monte Carlo Markov Chain methods [29, 20, 9] and the simulation of likelihood functions [17]. Moreover, the extensive simulation studies carried out in the papers just mentioned have showed that these LL schemes posses high numerical stability and remarkable computational efficiency. However, due to limitations in the numerical implementations of such LL schemes, they have been restricted to the class of autonomous SDEs with non singular Jacobian for the drift term. Moreover, the absence of a systematic theoretical study for these schemes has limited the study of the inference methods as well. Up to now, only the consistency of them has been studied in [28], while none estimate for the convergence rate has been given. On the other hand, an order two LL scheme recently introduced in [16] overcomes the restrictions of the previous schemes. However, this is based on quadrature formulas that, as will be shown in this paper, are inefficient and unnecessary. The main purpose of this paper is studying the convergence rate of both, the LL discretizations (2) and the first two LL schemes mentioned above. A second goal is introducing a new weak LL scheme that overcomes the drawbacks of the preceding LL schemes. Additionally, a simulation study is presented in order to illustrate the performance of the proposed integrator, which includes a comparison with conventional weak integrators as well.

2

Notations and Preliminaries Results

Let (Ω, F, P) be an underlying complete probability space and {Ft , t ≥ t0 } be an increasing right continuous family of complete sub σ-algebras of F associated to the m−dimensional standard Wiener process w defined in (1). Let M be the set of multi-indexes α = (j1 , ...jl(α) ), ji ∈ {0, 1, .., m}, i = 1, ..., l(α), where l(α) denotes the length of the multi-index α. Denote by −α and α− the multi-indexes in M that are obtained by deleting the first and the last component of α, respectively. The multi-index of length zero shall be denoted by ν. Let Γβ ⊂ M, β = 1, 2, be the hierarchical set Γβ = {α ∈ M : l(α) ≤ β} 2

and B(Γβ ) be the remainder set of Γβ , B(Γβ ) = {α ∈ M : l(α) = β + 1} . Denote by Hν , H(0) and H(1) the sets of adapted right continuous process h = {h(t), t ≥ 0} with left hand limits that satisfy Rt Rt |h(t)| < ∞, |h(s)| ds < ∞ and |h(s)|2 ds < ∞ w.p.1, 0

0

respectively. In addition define H(j) = H(1) for j = 2, ..., m, m ≥ 2. Then, for α = (j1 , ...jl(α) ), l(α) ≥ 2, define recursively the set Hα as the totality of adapted right continuous process h with left hand limits such that {Iα− [h(.)]ρ,t , t ≥ 0} ∈ Hjl(α) , where for h ∈ Hα the multiple Itô integral Iα [h(.)]ρ,τ is defined recursively by ⎫ ⎧ l(α) = 0 ⎬ ⎨ h(τ ), Rτ jl(α) Iα [h(.)]ρ,τ := . ⎩ Iα− [h(.)]ρ,s dWs , l(α) ≥ 1 ⎭ ρ

Let

L0 =

m d d P P ∂ ∂2 ∂ 1 P + fk k + Gk,i Gl,i k l ∂t k=1 ∂x 2 k,l=1 i=1 ∂x ∂x

be the diffusion operator of the SDE (1), and define Lj =

d P

k=1

Gk,j

∂ , j = 1, ..., m. ∂xk

Then, for each multi-index α = (j1 , ...jl(α) ) the Itô coefficient function λα is defined recursively by ¾ ½ Lj1 ...Ljl(α)−1 f , jl(α) = 0 . λα = Lj1 ...Ljl(α)−1 Gjl(α) , jl(α) 6= 0 The weak Itô-Taylor expansion of the process x solution of the SDE (1), corresponding to the hierarchical set Γβ , is given by P P Iα [λα (ρ, xρ )]ρ,τ + Iα [λα (., x.)]ρ,τ , (5) xτ = xρ + α∈Γβ /{ν}

α∈B(Γβ )

for any stopping times ρ, τ that satisfy 0 ≤ ρ ≤ τ ≤ T, w.p.1. Denote by CPl (Rd , R) the space of l time continuously differentiable functions g : Rd → R for which g and all its partial derivatives up to order l have polynomial growth. For l = 1, 2, ..., define Pl = {p ∈ {1, ..., d}l }, and for each p = (p1 , ..., pl ) ∈Pl define the function Fp : Rd → R as Fp (x) =

l Q

xpi .

i=1

We state below some results that will be referred to in section 4. Lemma 1 (Lemma 5.11.7 in [15]) Let β ∈ {1, 2, ...} be given and suppose that the components f k and Gk,j are such that 2(β+1) ([t0 , T ] × Rd , R), (6) f k , Gk,j ∈ CP

and satisfy Lipschitz conditions and linear growth bounds for all k = 1, ..., d . Then there exits constants K > 0 and r ∈ {1, 2, ...} such that ¯ ¯ ¯ ¯ P ¯ ¯ Iα [λα (t0 , x0 )]t0 ,t )ÁFt0 )¯ ≤ K(1 + |x0 |2r )(t − t0 )β+1 , ¯E(Fp (x(t) − x0 ) − Fp ( ¯ ¯ α∈Γβ /{ν}

for all l = 1, ..., 2β + 1 and p ∈ Pl .

3

The following theorem provides general conditions to assure that a discrete time approximation zδ converges weakly to x for some β. Theorem 2 (Theorem 14.5.2 in [15]) Let zδ be a time discrete approximation of the process x solution of (1) corresponding to a time discretization (t)δ . Suppose that E(|x0 |i ) < ∞, for i = 1, 2, .., ¯ ¯ ¯ ¯ ¯E(g(x0 )) − E(g(zδ0 ))¯ ≤ C0 δ β , 2(β+1)

for some β = 1, 2, .., a positive constant C0 and each g ∈ CP (Rd , R). Suppose that the drift coefficient is Lipschitz continuous with components that satisfy the condition (6) and the linear growth bound |f (t, x)| + |G(t)| ≤ K(1 + |x|), for all t ∈ [t0 , T ]. In addition suppose that for each q = 1, 2, ..., there exits constants K < ∞ and r ∈ N+ , which do not depend on δ, such that ¯ ¯2r ¯ ¯2q ¯ ¯ ¯ ¯ (7) E( max ¯zδtn ¯ ÁFt0 ) ≤ K(1 + ¯zδ0 ¯ ), 0≤n≤nT

¯ ¯2q ¯ ¯2r ¯ ¯ ¯ ¯ E(¯zδtn+1 − zδtn ¯ ÁFtn ) ≤ K(1 + max ¯zδtk ¯ )(tn+1 − tn )q , 0≤k≤n

(8)

and ¯ ¯ ¯ ¯ ¯ ¯2r P ¯ ¯ ¯ ¯ Iα [λα (tn , zδtn )]tn ,tn+1 )ÁFtn )¯ ≤ K(1 + max ¯zδtk ¯ )(tn+1 − tn )δ β , (9) ¯E(Fp (zδtn+1 − zδtn ) − Fp ( ¯ ¯ 0≤k≤n α∈Γβ /{ν} 2(β+1)

for all n = 0, 1, ..., nT − 1, l = 1, ..., 2β + 1and p ∈ Pl . Then for any g ∈ CP ¯ ¯ ¯ ¯ ¯E(g(xT )) − E(g(zδT ))¯ ≤ Cg δ β ,

(Rd , R),

(10)

for some positive constant Cg .

3

Weak Local Linear Discretizations

Let (t)δ = {t0 ≤ t1 ≤ ... ≤ tn < ... < ∞} be a time partition defined as a sequence of Ftn -measurable random times tn , n = 0, 1, .., that satisfy sup(δ n ) ≤ δ < 1, w.p.1, n

where δ n = tn+1 − tn . Define nt := max{n = 0, 1, 2, ..., : tn ≤ t} < ∞. A precise definition of LL discretization is the following. Definition 3 For a given time discretization (t)δ , the order β(= 1, 2) weak LL Discretization of the solution of (1) is defined by the recurrent relation (2), where the initial point y0 is a Ft0 -measurable random vector.

4

As was early mentioned in the introduction, depending on the way of computing the integrals (3) and (4) in (2) different LL schemes can be obtained. For instance, by integrating by part in (3) it is obtained

ytn+1

¶ µ ∂f (y ) ¶ tn ∂f (ytn ) −1 δ n = ytn + − Id f (ytn ) e ∂y ∂y ¶ µ ∂f (y ) ¶ µ tn ) ∂f (y ∂f (ytn ) −2 t δ n n δ n bβ (tn , ytn ) e ∂y − Id − + ∂y ∂y +η(tn , ytn ; δ n ), µ

(11)

∂G which for β = 1 and β = 2 corresponds to the LL schemes for autonomous equations (i.e. ∂f ∂t = ∂t ≡ 0) proposed in [18] and [23], [24] , respectively. For scalar SDEs, the variance of η(tn , ytn ; δ n ) is given by

¶ µ µ ¶ tn ) df (ytn ) −1 2 df (y Σ(tn , ytn ; δ n ) = 2 e dy δn − 1 G2 , dy which is obtained by integrating by part in (4) [18], [23]. For multidimensional SDEs the variance Σ of η(tn , ytn ; δ n ) is obtained as solution of the pencil equation [24] µ

∂f (ytn ) ∂y

¶

µ

∂f (ytn ) Σ+Σ ∂y

¶>

=e

∂f (ytn ) δn ∂y

|

 ∂f (y

GG e

tn ) ∂y

>

δn

− GG| .

(12)

Notice that the numerical scheme (11) is not always computational feasible since it can fails for SDEs where ∂f be either singular or ill-conditioned at some point ytn . Moreover, the equation the Jacobian matrix ∂y ∂f (12) might have no unique solution for some particular ∂y [24]. Another weak LL scheme for the nonautonomous case is the following [16]: ytn+1 = e

∂f (tn ,ytn ) δn ∂y

ytn + δ n e

∂f (tn ,ytn ) δn ∂y 2

(f (tn , ytn ) −

∂f (tn , ytn ) δn ytn + b2 (tn , ytn ) ) + η(tn , ytn ; δ n ), ∂y 2

where the variance Σ of η(tn , ytn ; δ n ) is approximated by Σ(tn , ytn ; δ n ) ≈ δ n e

∂f (tn ,ytn ) δ n ∂y 2

 ∂f (t

δn δn G(tn + )G(tn + )| e 2 2

n ,ytn ) ∂y

>

δn 2

.

This scheme can be obtained by making some algebraic manipulations in (2) and using quadrature formulas to compute the integrals (3) and (4) for the case β = 2. Note that this scheme overcomes the restrictions of the previous ones, but it requires an approximation to evaluate the integrals just mentioned. Next definition will be useful to obtain the convergence rate of the LL discretizations and numerical schemes in the next section. Definition 4 For a given time discretization (t)δ , the stochastic process yβδ = {yβδ (t), t ∈ [t0 , T ]} is called the order β(= 1, 2) weak LL Approximation of the solution of (1) if yβδ (t) = ytnt + φβ (tnt , ytnt ; t − tnt ) + η(tnt , ytnt ; t − tnt ),

(13)

where {ytn }, n = 0, 1, ..., is the LL Discretization (2). Note that, the LL Approximation (13) is a continuous time stochastic process that coincides with the LL Discretization (2) at each discretization time tn ∈ (τ )δ . 5

In addition, it is convenient to remark that the weak LL Approximation (13) coincides with the weak solution of the piecewise stochastic differential equation [11] dy(t) = pβ (t, y(t); tn , y(tn ))ds + G(t)dw(t), t ∈ [tn , tn+1 ], n = 0, 1, ..., nT − 1,

(14)

y(tn ) = ytn ,

where the function pβ is defined as pβ (s, v; r, u) = f (r, u) +

∂f (r, u) (v − u) + bβ (r, u)(s − r), for all v, u ∈ Rd , s, r ∈ R, s > r, ∂u

(15)

which for β = 1 and β = 2 is just the first order Taylor and Itô Taylor expansion of f , respectively.

4

Convergence of Weak Local Linear Discretizations

In this section, conditions for the weak convergence of the LL Approximation yβδ to the solution of the equation (1) are presented. The main result is stated in the following theorem. Theorem 5 Suppose that E(|x0 |i ) < ∞, for i = 1, 2, .., ¯ ¯ ¯ ¯ ¯E(g(x0 )) − E(g(yβδ (t0 )))¯ ≤ C0 δ β , 2(β+1)

for some C0 > 0 and all g ∈ CP constant such that

e be a positive (Rd , R). Assume also that condition (6) holds and let K

e + |x|), |f (t, x)| + |G(t)| ≤ K(1 ¯ ¯ ¯ ¯ ¯ ¯ 2 ¯ ∂f (t, x) ¯ ¯ ∂f (t, x) ¯ ¯ ∂ f (t, x) ¯ e ¯ ¯ ¯ ¯ ¯ ¯ ¯ ∂t ¯ + ¯ ∂x ¯ + ¯ ∂x2 ¯ ≤ K,

(16) (17)

for all t ∈ [t0 , T ] and x ∈Rd . Then there exits a positive constant Cg such that the LL Approximation yβδ satisfies ¯ ¯ ¯ ¯ δ E(g(x(T ))) − E(g(y (T ))) (18) ¯ ≤ Cg δ β . ¯ β

In order to proof the theorem above two additional lemmas shall be enunciated. The first lemma establishes uniform bounds for the moments of the LL Approximation (13), whereas the second one provides an Itô-Taylor expansion for the stochastic process that defines such LL Approximation. In all the proofs of this section, K and r will denote a positive constant and a natural number that do not depend on h, respectively. Whenever either K or r be used, they will only depend on the Lipschitz constants and growth bounds for f , G and its derivatives. For simplicity in the presentation, the values of them are not explicitly given from one occurrence to the other one. Lemma 6 Under the assumptions of Theorem 5 the LL Approximation yβδ , β = 1, 2 satisfies ¯ ¯2q ¯ ¯ E( sup ¯yβδ (t)¯ ÁFt0 ) ≤ K(1 + |y0 |2q ), t0 ≤t≤T

for each q = 1, 2, ..., where K is a positive constant.

6

Proof. The proof of this lemma is very similar to that of Theorem 7.1.2 in [1]. For simplicity in the exposition, the superscript δ shall be omitted in the subsequent. For q = 1 this is just Theorem 5 in [11]. For each q = 2, 3, ..., the Itô formula and equation (14) yield to the following piecewise SDE, |y(t)|2q = |y(tnt )|2q +

Rt

2q |y(u)|2q−2 y| (u)pβ (u, y(u); tnt , y(tnt ))du +

tnt

G(u)dw(u) +

Rt

tnt

q |y(u)|2q−2 |G(u)|2 du +

Rt

tnt

Rt

tnt

2q |y(u)|2q−2 y| (u)

2q(q − 1) |y(u)|2q−4 |y| (u)G(u)|2 du,

for all t ∈ [tnt , tnt +1 ]. By recursive application of this expression it is obtained, |y(t)|2q = |y0 |2q +

nP t −1 tn+1 R n=0

tn

2q |y(u)|2q−2 y| (u)pβ (u, y(u); tn , y(tn ))du +

pβ (u, y(u); tnt , y(tnt ))du +

Rt

t0

Rt

Rt

2q |y(u)|2q−2 y| (u)

tnt

2q |y(u)|2q−2 y| (u)G(u)dw(u) +

+ 2q(q − 1) |y(u)|2q−4 |y| (u)G(u)|2 du,

Rt

t0

q |y(u)|2q−2 |G(u)|2 du

t0

for all t ≥ t0 . This in turn yields to E( sup |y(s)|2q ÁFt0 ) ≤ |y0 |2q + t0 ≤s≤t

Rt

tnt

Rt

t0

Rt

t0

nP t −1 tn+1 R n=0

tn

2qE(|y(u)|2q−2 |y| (u)pβ (u, y(u); tn , y(tn ))|ÁFt0 )du +

2qE(|y(u)|2q−2 |y| (u)pβ (u, y(u); tnt , y(tnt ))|ÁFt0 )du +

2q(q − 1)E(|y(u)|2q−4 |y| (u)G(u)|2 ÁFt0 )du + qE(|y(u)|2q−2 |G(u)|2 ÁFt0 )du.

From conditions (16) and (17), and expression (15), the following inequality is derived |pβ (u, y(u); tnu , y(tnu ))| ≤

sup |pβ (s, y(s); tnu , y(tnu ))|

tnu ≤s≤u

e + e sup |y(s)| + K(1 ≤ 2K tnu ≤s≤u

where

Thus,

⎧ e ⎨ K, e Kβ = e ⎩ K(1 +

d2 2

eβ , sup |y(s)|) + K

tnu ≤s≤u

⎫ β=1 ⎬ sup |G(s)|2 ), β = 2 ⎭ ,

t0 ≤s≤T

|y| (u)pβ (u, y(u); tnu , y(tnu ))| ≤ K(1 +

7

sup |y(s)|2 ),

tnu ≤s≤u

for some positive constant K. This implies that E( sup |y(s)|2q ÁFt0 ) ≤ |y0 |2q + 2qK t0 ≤s≤t

Rt

t0

Rt

t0

Rt

t0

E( sup (1 + |y(s)|2 ) |y(s)|2q−2 ÁFt0 )du + t0 ≤s≤u

2q(q − 1)E( sup |y(s)|2q−4 |y| (s)G(s)|2 ÁFt0 )du + t0 ≤s≤u

qE( sup |y(s)|2q−2 |G(s)|2 ÁFt0 )du. t0 ≤s≤u

Moreover, from (16), E( sup |y(s)|2q ÁFt0 ) ≤ |y0 |2q + 2qK t0 ≤s≤t

Rt

t0

E( sup (1 + |y(s)|2 ) |y(s)|2q−2 ÁFt0 )du +

Rt

t0 ≤s≤u

e 2 q(2q − 1) E( sup (1 + |y(s)|2 ) |y(s)|2q−2 ÁFt )du, 2K 0 t0

and from inequality (1 + a2 )a2q−2 ≤ 1 + 2a2q ,

t0 ≤s≤u

e 2 (2q − 1))(t − t0 ) + E( sup |y(s)|2q ÁFt0 ) ≤ |y0 |2q + 2q(K + K t0 ≤s≤t

Rt e 2 (2q − 1)) E( sup |y(s)|2q ÁFt )du. 4q(K + K 0 t0

t0 ≤s≤u

The proof concludes by applying the Gronwall Lemma to this inequality.

Lemma 7 Let yβδ , β = 1, 2 be the LL Approximation (13) and zδβ = {zδβ (t), t ∈ [0, T ]} be the stochastic process defined by P P Iα [Λα (tnt , ytnt ; tnt , ytnt )]tnt ,t + Iα [Λα (., y.; tnt , ytnt )]tnt ,t , zδβ (t) = ytnt + α∈Γβ /{ν}

where Λα (s, v; r, u) =

α∈B(Γβ )

½

Lj1 ...Ljl(α)−1 pβ (s, v; r, u), Lj1 ...Ljl(α)−1 Gjl(α) ,

jl(α) = 0 jl(α) 6= 0

¾

and pβ (s, v; r, u) is defined as in (15). Then yβδ ≡ zδβ and Iα [Λα (tnt , ytnt ; tnt , ytnt )]tnt ,t = Iα [λα (tnt , ytnt )]tnt ,t ,

(19)

for all α ∈ Γβ /{ν} and t ∈ [t0 , T ]. Proof. Let Aγ be the hierarchical set 1 Aγ := {α ∈ M : l(α) + n(α) ≤ 2γ or l(α) = n(α) = γ + }, 2 where γ = β+1 2 , β = 1, 2, and let n(α) be the number of null components in the multi-index α. The thesis of the lemma is straighforwardly obtained from Lemma 7 in [11] by realizing that Γβ ⊂ Aγ , for β = 1, 2. 8

Proof of Theorem 5. By Lemma 6, E( max |ytn |2q ÁFt0 ) ≤ K(1 + |y0 |2r ),

(20)

0≤n≤nT

which is just the condition (7) of Theorem 2. The application of Theorem 4.5.4 in [15] to the piecewise linear equation (14) yields to ¯2q ¯ (21) E(¯ytn+1 − ytn ¯ ÁFtn ) ≤ K(1 + max |ytk |2r )(tn+1 − tn )q , 0≤k≤n

i.e., the condition (8). On the other hand, by applying Lemma 1 to the equation (14) it is obtained ¯ ¯ ¯ ¯ P ¯ ¯ Iα [Λα (tn , ytn ; tn , ytn )]tn ,tn+1 )ÁFtn )¯ ≤ K(1 + |ytn |2r )(tn+1 − tn )δ β , ¯E(Fp (ytn+1 − ytn ) − Fp ( ¯ ¯ α∈Γβ /{ν}

for all p ∈ Pl and l = 1, ..., 2β + 1 and some K > 0 and r ∈ {1, 2, ...} . Since by Lemma 7 the equality P P Iα [λα (tn , ytn )]tn ,tn+1 )ÁFtn ) = E(Fp ( Iα [Λα (tn , ytn ; tn , ytn )]tn ,tn+1 )ÁFtn ), E(Fp ( α∈Γβ /{ν}

α∈Γβ /{ν}

holds, then ¯ ¯ ¯ ¯ P ¯ ¯ Iα [λα (tn , ytn )]tn ,tn+1 )ÁFtn )¯ ≤ K(1 + |ytn |2r )(tn+1 − tn )δ β ¯E(Fp (ytn+1 − ytn ) − Fp ( ¯ ¯ α∈Γβ /{ν}

≤ K(1 + max |ytk |2r )(tn+1 − tn )δ β , 0≤k≤n

which is just the condition (9). The proof concludes by applying Theorem 2. Theorem (5) provides the global order of weak convergence for the LL Approximations at the time t = T . Thus, obviously, it also provides the order of convergence for the LL Discretization (2) and the scheme (11). Notice further that inequality (18) implies that the uniform bound ¯ ¯ ¯ ¯ sup ¯E(g(x(t))) − E(g(yβδ (t)))¯ ≤ Cg δ β t∈[t0 ,T ]

also holds for the LL Approximations yβδ , since the order condition (10) in Theorem 2 also holds uniformly in [t0 , T ] (see Theorem 14.5.1 and Exercise 14.5.3 in [15] for details).

5

Numerical Scheme

In this section a new weak LL scheme is presented and its rate of weak convergence is also studied. For simplicity, the cases of autonomous and nonautonomous G shall be treated separately.

5.1

Autonomous Case

To fix ideas, firstly consider the case of autonomous G, i.e. G(t) ≡ G. According to Lemma 9 in Appendix, the terms φβ and Σ in the LL Discretization (2) can be rewritten as φβ (tn , ytn ; δ n ) =

δRn

e

∂f (tn ,ytn ) (δ n −s) ∂y

f (tn , ytn )ds +

0

δRn Rs

e

∂f (tn ,ytn ) (δ n −s) ∂y

bβ (tn , ytn )ds1 ds

0 0

= B14 (tn , ytn ; δ n ), δRn ∂f (tn ,ytn ) ∂f (tn ,ytn ) | ∂f (tn ,ytn ) | s δn Σ(tn , ytn ; δ n ) = ( e ∂y (δn −s) GG| e− ∂y ds) e ∂y 0

= B12 (tn , ytn ; δ n )B|11 (tn , ytn ; δ n ), 9

where the block matrix B =(Blj ), l, j = 1, ..., 4 is defined as B(tn , ytn ;δ n ) = eA(tn ,ytn )δn , with ⎛ ∂f (t

⎜ ⎜ A(tn , ytn ) = ⎜ ⎝

GG|

n ,ytn ) ∂y

0 0 0

n ,ytn ) − ∂f (t∂y 0 0

Therefore, it is obtained the following LL scheme

|

⎞ bβ (tn , ytn ) f (tn , ytn ) ⎟ ⎟ 0 0 ⎟. ⎠ 0 1 0 0

yn+1 = yn + B14 (tn , ytn ; δ n ) + (B12 (tn , ytn ; δ n )B|11 (tn , ytn ; δ n ))1/2 ξ n ,

(22)

starting with y0 = x0 . Here, Σ1/2 denotes the square root matrix of Σ and {ξ n } is a sequence of ddimensional i.i.d normal random vectors. Obviously, the rate of convergence of the scheme (22) is provided by Theorem 5. It is well-known (see Section 5-12 in [15], for instance) that for purposes of weak convergence, the sequence {ξ n } of i.i.d normal random vectors in the scheme above can be replaced by other random variables that satisfy corresponding moment conditions. Thus, it is also obtained the simplified LL scheme ξn, yn+1 = yn + B14 (tn , ytn ; δ n ) + (B12 (tn , ytn ; δ n )B|11 (tn , ytn ; δ n ))1/2b

where {b ξ n } is a sequence of independent random variables with independent two-point distributed comk ponents {b ξ n } that satisfy k 1 P (b ξ n = ±1) = . 2 It should be noted that the LL scheme (22) is computationally feasible and its numerical implementation is reduced to the use of a suitable algorithm to compute matrix exponentials, e.g., those based on rational Padé approximations [7], the Schur decomposition [7] or Krylov subspace methods [26]. The selection of one of them will mainly depend on the size and structure of the matrix A. For instance, for many low dimensional systems of equations it is enough to use the algorithm developed in [33], which takes advantage of the special structure of the matrix A. Whereas, for large systems of equations, the Krylov subspace methods are strongly recommended. Notice that in contrast with the previous LL integrators (11), the LL scheme (22) has no restriction ∂f . Note further that, by using the Lemma 9, the integrals (3) and (4) were on the Jacobian matrix ∂y exactly computed in terms of a matrix exponential. This explains why the use of quadrature formulas to define weak LL schemes is unneeded. Moreover, they are not recommended because they introduce an unnecessary lack of precision.

5.2

Nonautonomous Case

For the case of nonautonomous G, consider the following approximation: G(tn + s) ≈ Gβ (tn + s), s ∈ [0, δ n ], where for each Gβ denotes the order β − 1 truncated Taylor expansion of G given by Gβ (tn + s) =

β−1 P i=0

This in turn implies that

di G(tn ) i s. dti

Σ(tn , ytn ; hn ) ≈ Σβ (tn , ytn ; δ n ), 10

where hRn

Σβ (tn , ytn ; δ n ) = (

e

∂f (tn ,ytn ) (δ n −s) ∂y

0 2β−2 P hRn Rs sR0

....

i=1 0 0 0

with

si−2 R

H0 (tn )e−

e

∂f (tn ,ytn ) | s ∂y

∂f (tn ,ytn ) (δ n −s) ∂y

ds +

Hi (tn )e−

∂f (tn ,ytn ) | s ∂y

dsi−1 ...ds0 ds)e

∂f (tn ,ytn ) | δn ∂y

,

0

P dl G(tn ) dj G(tn ) | Hi (tn ) = , i = 0, ..., 2β − 2. dtl dtj l+j=i

Hence, again by Lemma 9 it is obtained

φβ (tn , ytn ; δ n ) = B1,2β+2 (tn , ytn ; δ n ),

Σβ (tn , ytn ; δ n ) = B1,2β (tn , ytn ; δ n )B|11 (tn , ytn ; δ n ),

(23)

where eA(tn ,ytn )δn = B(tn , ytn ; δ n ) and ⎞ ⎛ ∂f (t ,y ) n tn H2β−2 (tn ) H2β−3 (tn ) · · · H0 (tn ) bβ (tn , ytn ) f (tn , ytn ) ∂y | ⎟ ⎜ n ,ytn ) ⎟ ⎜ 0 − ∂f (t∂y Id ··· 0 0 0 ⎟ ⎜ ⎟ ⎜ .. .. .. .. .. .. ∂f (tn ,ytn ) | ⎟ ⎜ . . . − ∂y . . . ⎟ ⎜ ⎟ ⎜ .. .. .. A(tn , ytn ) = ⎜ ⎟. . Id 0 0 0 . . ⎟ ⎜ ⎟ ⎜ .. .. .. .. ∂f (tn ,ytn ) | ⎟ ⎜ . 0 0 . . . − ∂y ⎟ ⎜ ⎟ ⎜ .. .. .. ⎠ ⎝ . . . ··· 0 0 1 0 0 0 ··· 0 0 0 This gives the following LL scheme

en+1 = y en + B1,2β+2 (tn , ytn ; δ n ) + (B1,2β (tn , ytn ; δ n )B|11 (tn , ytn ; δ n ))1/2 ξ n , y

(24)

and the simplified LL scheme

en + B1,2β+2 (tn , ytn ; δ n ) + (B1,2β (tn , ytn ; δ n )B|11 (tn , ytn ; δ n ))1/2b en+1 = y ξn, y

e0 = x0 . starting with y It is easily seen that the scheme (24) reduces to scheme (22) in the case G(t) ≡ G, whereas both weak schemes reduce to the deterministic LL scheme [12, 10] when G(t) ≡ G = 0. Note further, expression (23) provides a more efficient way to compute Σβ than that proposed in [14] based on Schur decomposition. In contrast with the autonomous case, Theorem 5 can not be applied to state the rate of convergence of the scheme (24). Therefore, the following theorem is in order. Theorem 8 Suppose that E(|x0 |i ) < ∞, for i = 1, 2, ..,

yβ (t0 )))| ≤ C0 δ β , |E(g(x0 )) − E(g(e 2(β+1)

(Rd , R). Assume also that the growth condition (16) holds. Then, there for some C0 > 0 and all g ∈ CP eβ satisfies exits a positive constant Cg such that the truncated LL Approximation y |E(g(x(T ))) − E(g(e yβ (T )))| ≤ Cg δ β . 11

Proof. Firstly, it should be noted that etn + φβ (tn , y etn ; δ n ) + Σβ (tn , y etn ; δ n )1/2 ξ n etn+1 = y y

(25)

is a weak solution at t = tn+1 of the piecewise linear SDE

e (t); tn , y e(tn ))dt + G(t)dw(t), t ∈ [tn , tn+1 ] de y(t) = pβ (t, y e(tn ) = y etn y

Hence, as in Theorem 5, it is straigthforward to show that

E( max |e ytn |2q ÁFt0 ) ≤ K(1 + |y0 |2r ), 0≤n≤nT

and

¯ ¯2q etn+1 − y etn ¯ ÁFtn ) ≤ K(1 + max |e E(¯y ytk |2r )(tn+1 − tn )q . 0≤k≤n

On the other hand, ¯ ¯ ¯ ¯ P ¯ ¯ etn ) − Fp ( etn )]tn ,tn+1 )ÁFtn )¯ ≤ e1 + e2 , ytn+1 − y Iα [λα (tn , y ¯E(Fp (e ¯ ¯ α∈Γβ /{ν} where

and

¯ ¯ ¯ ¯ P ¯ ¯ etn ) − Fp ( etn )]tn ,tn+1 )ÁFtn )¯ , e1 := ¯E(Fp (ztn+1 − y Iα [λα (tn , y ¯ ¯ α∈Γβ /{ν} ¯ ¯ etn ) − Fp (ztn+1 − y etn )ÁFtn )¯ e2 := ¯E(Fp (e ytn+1 − y etn + φβ (tn , y etn ; δ n ) + Σ(tn , y etn ; δ n )1/2 ξ n ztn+1 = y

(26)

denotes a weak solution at t = tn+1 of the SDE

etn )dt + G(t)dw(t), t ∈ [tn , tn+1 ] dz(t) = pβ (t, z(t); tn , y etn z(tn ) = y

Then, by applying Lemma 1 to the equation above it is obtained ¯ ¯ ¯ ¯ P ¯ ¯ etn ) − Fp ( etn )]tn ,tn+1 )ÁFtn )¯ ≤ K(1 + |e ytn |2r )δ β+1 , Iα [Λα (tn , y ¯E(Fp (ztn+1 − y ¯ ¯ α∈Γβ /{ν}

which by Lemma 7 is equivalent to

e1 ≤ K(1 + |e ytn |2r )δ β+1 .

From Lemma 10 follows that l(p)−1 ¯2 ¯4j ¯l(p)−1−j ¯ ¯ ¯ P etn+1 ¯ ÁFtn ))1/2 etn+1 ¯ e2 ≤ (E(¯ztn+1 − y (E(¯ztn+1 ¯ ÁFtn ))1/4 (E(¯y ÁFtn ))1/4 j=0

¯2 ¯ etn+1 ¯ ÁFtn ))1/2 . ≤ K(1 + |e ytn | )(E(¯ztn+1 − y 2r

12

(27)

By using the expressions (25) and (26) it is obtained ¯ ¯2 ¯2 ¯ ¯ ¯ etn+1 ¯ ÁFtn ))1/2 = (E(¯(Σ(tn , y etn ; δ n )1/2 − Σβ (tn , y etn ; δ n )1/2 )ξ n ¯ ÁFtn ))1/2 (E(¯ztn+1 − y

etn ; δ n )1/2 − Σβ (tn , y etn ; δ n )1/2 ]| [Σ(tn , y etn ; δ n )1/2 = (tr([Σ(tn , y etn ; δ n )1/2 ]E(ξ n ξ |n )))1/2 −Σβ (tn , y ¯ ¯ ¯ ¯ etn ; δ n )1/2 − Σβ (tn , y etn ; δ n )1/2 ¯ , = ¯Σ(tn , y

where tr() denotes the trace operator. On the other hand, by using the expressions for Σ and Σβ and some algebraic manipulations follows that ¯ ¯ ¯ ¯ | | etn ; δ n ) − Σβ (tn , y etn ; δ n )| ≤ Kδ n sup ¯G(tn + s)G (tn + s) − Gβ (tn + s)Gβ (tn + s)¯ |Σ(tn , y s∈[0,δ n ]

≤ Kδ n sup (|G(tn + s)| + |Gβ (tn + s)|) |G(tn + s) − Gβ (tn + s)| s∈[0,δ n ]

≤ Kδ

β+1

.

Now, the perturbation bounds for the Cholesky and SVD factorizations stated in [30] (Theorems 2.2.1 and 3.2.1 ) yield to ¯ ¯ ¯ 1/2 1/2 ¯ etn ; δ n ) − Σβ (tn , y etn ; δ n ) ¯ ≤ K |Σ(tn , y etn ; δ n ) − Σβ (tn , y etn ; δ n )| ¯Σ(tn , y ≤ Kδ β+1 ,

which gives ytn |2r )δ β+1 . e2 ≤ K(1 + |e Hence, this and (27) yield to ¯ ¯ ¯ ¯ P ¯ ¯ etn ) − Fp ( etn )]tn ,tn+1 )ÁFtn )¯ ≤ K(1 + |e ytn |2r )δ β+1 , ytn+1 − y Iα [λα (tn , y ¯E(Fp (e ¯ ¯ α∈Γβ /{ν} and the proof concludes by applying Theorem 2.

6

Numerical Experiments

In this section two examples are given in order to illustrate the performance of the introduced LL scheme. Example 1. Consider the 2-dimensional SDE µ ¶ µ ¶ µ ¶ 0 1 x1 = dw(t), t ∈ [0, 1] (28) d 1 2 dt + x2 x 0 2 1 µ ¶ 0 x(0) = , 0 which corresponds to Example 17.1.1 in [15]. It turns out (see Exercise 17.1.2 in [15]) that for the function g(x) = e−x2 , − 12

E(g(x(1))) = E(e r

R1

(w(s))2 ds

0

2e . 1 + e2

=

13

)

The quantity eb(δ) = |E(g(x(1))) − E(g(y(1)))|

was used to estimate the order β of weak convergence of the LL schemes (22). Here, y(t), t ∈ (t)δ denotes a simulated trajectory of equation (28) computed by (22) with step size δ. Then, the value E(g(y(1))) was estimated by the average of 10000 trajectories. A number of 2000 of such simulations were carried out. They we arranged into M = 20 batches with K = 100 simulations of E(g(y(1))) each, which were denoted by ebi,j (δ), i = 1, .., M ; j = 1, ..., K. Then, computing the sample mean error of the i-th batch and of all batches by K M 1 X 1 X ebi,j (δ), and eb(δ) = ebi (δ), ebi (δ) = K M j=1

i=1

respectively. The confidence interval for eb(δ) was computed as

[b e(δ) − ∆(δ), eb(δ) + ∆(δ)],

where

∆(δ) = t1−α/2,M −1

s

M

σ b2e (δ) 2 1 X , σ be (δ) = |b ei (δ) − eb(δ)|2 , M M −1 i=1

and t1−α/2,M −1 denotes the 1 − α/2 percentile of the Student’s t distribution with M − 1 degrees for the significance level 0 < α < 1. Table I shows the estimated values eb(δ) and their respective 90% confidence interval (i.e. the values ∆(δ) for α = 0.1). δ 2−4 2−5 2−6 2−7

Order1 LL (β = 1) eb(δ) ∆(δ) 0.00707225 0.00004896 0.00334453 0.00004251 0.00161903 0.00003934 0.00092813 0.00002245

Order 2 LL (β = 2) eb(δ) ∆(δ) 0.00176849 0.00002708 0.00035382 0.00000944 0.00008744 0.00000441 0.00002991 0.00000131

Table I. Global Discretization Errors and their respective 90% confidence intervals for the LL scheme (22) with β = 1 and β = 2. b of weak convergence was obtained from Table I as the slope For each β = 1, 2 the estimated order β e(δ i ))} , δ i = 2−(i+3) , i = 1, ..., 4. Figure 1 of the straight line fitted to the set of points {log2 (δ i ), log2 (b b = 0.9835 ± 0.0236 and β b shows the plot of such straight lines, where the values for the slopes are β = 1.9673 ± 0.0411, respectively. Notice that these results corroborate the theoretical estimate for β given in Theorem 5. Example 2. Consider the 1-dimensional SDE given by (Example 2 in [2]) 3

3e−t /3 dw(t), dx = −t x(t)dt + 2(t + 1) x(0) = 1, 2

whose exact solution satisfies

3 /3

E(x(t)) = e−t 14

.

t ∈ [0, 10]

(29)

For different numerical schemes, a number of N = 100 simulations yi (t), t ∈ (t)δ , δ = 2−4 , i = 1, ..., N N P yi (t), t ∈ (t)δ . Figure 2 shows the were used for estimating E(x(t)) on the basis of the average N1 i=1

comparison of this result for the order 2 explicit extrapolated weak Euler scheme (EE) [32], the order 3 explicit weak Taylor scheme (ET3) [15], the order 1 LL scheme (LL1) (given by (24) with β = 1) and the implicit weak Euler (IE) [15]. Notice that even higher order schemes such as the extrapolated Euler scheme and the order 3 weak Taylor scheme give explosive trajectories after time t = 9. In contrast, the LL scheme and implicit Euler scheme present better numerical stability. However, the implicit Euler method involves an additional computational effort due to the fact that it requires the solution of a nonlinear algebraic equation at each time step. This can be seen in Table II, which shows the relative CPU time for each numerical scheme, which is obtained as the ratio of the actual CPU time for each numerical scheme to the minimum of all CPU times. Notice that the CPU time for the LL method is around 33 times lower than that for the implicit Euler method. Therefore, the LL method shows not only numerical stability but also a low computational cost when compared with implicit methods that also behave numerically stable. Scheme Relative CPU time

ET3 1

EE 1.36

LL1 2.15

IE 71.29

Table II. Relative CPU times for different numerical schemes (order 3 Explicit Taylor (ET3), order 2 extrapolated Explicit Euler (EE), order 1 Local Linearization (LL1) and Implicit Euler (IE)).

7

Conclusions

The order of convergence of the weak Local Linear Discretization for the numerical integration of stochastic differential equations with additive noise was obtained. Additionally, a new weak LL scheme was introduced and its order of weak convergence was studied as well. Some advantages of the new scheme over previousreported LL schemes were pointed out. The simulation study carried out demonstrates that the proposed LL scheme has better accuracy and stability properties than some standard weak integrators.

8

Appendix

The following Lemma is a generalization of Theorem 1 in [33]. Lemma 9 (Theorem 1 in [4]) Let n, d1 , d2 , ..., dn be matrix defined by ⎛ A11 A12 ⎜ 0 A22 ⎜ A =⎜ ⎝ 0 0 0 0

positive integers and A a n × n block triangular ... ... .. . 0

⎞ A1n A2n ⎟ ⎟ .. ⎟ , . ⎠

Ann

where (Alj ) are dl × dj matrices, with l, j = 1, ..., n . Then for t > 0 ⎛ ⎞ B11 (t) B12 (t) ... B1n (t) ⎜ 0 B22 (t) ... B2n (t) ⎟ ⎜ ⎟ eAt = ⎜ .. ⎟ , .. ⎝ 0 . 0 . ⎠ 0 0 0 Bnn (t)

15

with Bll (t) = eAll t , l = 1, ..., n j−l−1 Rt (l,j) P Rt sR1 sRk Blj (t) = M (t, s1 )ds1 + ... k=1 0 0

0

P

l = 1, ..., n − 1, j = l + 1, ..., n,

where M(i1 ,...,ik ) (s1 , ..., sk ) =

M(l,i1 ,...,ik ,j) (t, s1 , ..., sk+1 )dsk+1 ...ds1 ,

0 l
µk−1 Q r=1

¶ eAir ir (sr −sr+1 ) Air ir+1 eAik ik sk .

Lemma 10 If p = (p1 , ..., pl ) ∈ Pl then |Fp (x) − Fp (y)| ≤ for all x, y ∈ Rd .

l−1 P

i=0

|x|i |y|l−1−i |x − y| ,

Proof. It is very easy to see that |Fp (x) − Fp (y)| ≤ |x| |F−p (x) − F−p (y)| + |x − y| |F−p (y)| ≤ |x| |F−p (x) − F−p (y)| + |x − y| |y|l−1

Thus, the results follows from applying recursively the inequality above.

References [1] Arnold, L., Stochastic Differential Equations: Theory and Applications, Wiley Interscience Publications, New York, 1974. [2] Biscay, R., Jimenez, J.C., Riera, J. and Valdes, P., Local linearization method for the numerical solution of stochastic differential equations, Annals Inst. Statist. Math., 1996, 48, 631-644. [3] Blankenship, G.I. and Baras, J.S., “Accurate evaluation of stochastic Wiener integrals with applications to scattering in random media and nonlinear filtering, SIAM J. Appl. Math., 1981, 41, 518-552. [4] Carbonell, F., Jimenez, J.C. and Pedroso, L., “Computing multiple integrals involving matrix exponentials, 2005, Technical Report # 2005-328, ICIMAF, Havana, Cuba. [5] Durham, G.B. and Gallant, A.R., Numerical techniques for maximum likelihood estimation of continuous-time diffusion processes, J. Buss. Econ. Stats., 2002, 20(3), 297-316. [6] Ghysels, E., Harvey, A.C. and Renault, E., Stochastic volatility,1996, In Maddala G.S. and Rao C.R. (eds.), Handbook of Statistics, Vol. 14, 19-191. [7] Golub, G. H. and Van Loan, C. F., Matrix Computations, The Johns Hopkins University Press, 2nd Edition, 1989. [8] Grorud, A. and Talay, D., “Approximation of Lyapunov exponents of nonlinear stochastic differential equations”, SIAM J. Appl. Math., 1996, 56, 627-650. 16

[9] Hansen, N.R., “Geometric ergodicity of discrete-time approximations to multivariate diffusions”, Bernoulli, 2003, 9, 725-743. [10] Jimenez, J. C., “A simple algebraic expression to evaluate the Local Linearization schemes for stochastic diferential equations”, Appl. Math. Lett., 2002, 15, 775-780. [11] Jimenez, J. C. and Biscay, R., “Approximation of continuous time stochastic processes by the local linearization method revisited ”, Stochastic Anal. Appl., 2002, 20, 105-121. [12] Jimenez, J.C., Biscay, R., Mora, C.M. and Rodriguez, L.M., “Dynamic properties of the Local Linearization method for initial-value problems”, Appl. Math. Comput., 2002, 131, 21-37. [13] Jimenez, J.C., Shoji, I. and Ozaki, T. “Simulation of Stochastic Differential Equations Through the Local Linearization Method. A Comparative Study”, J. Stat. Phys., 1999, 94, 587-602. [14] Jimenez, J. C., Valdes, P., Rodriguez, L. M., Riera, J. and Biscay, R., “Computing the noise covariance matrix of the Local Linearization scheme for the numerical solution of stochastic differential equations”, Appl. Math. Lett., 1998, 11, 19-23. [15] Kloeden, P. E., and Platen, E., Numerical solution of stochastic differential equations, 2nd Edition, Springer-Verlag, Berlin, 1995. [16] Mora, C. M., “Weak exponential schemes for stochastic differential equations with additive noise”, IMA J. Numer. Anal., 2005, 25, 486-506. [17] Nicolau, J., “A new technique for simulating the likelihood of stochastic differential equations”, Econom. J., 2002, 5, 91-103. [18] Ozaki, T., “A bridge between nonlinear time series models and nonlinear stochastic dynamical systems: a local linearization approach”, Statistica Sinica, 1992, 2, 113-135. [19] Prakasa-Rao, B.L.S., Statistical inference for diffussion type processes, Oxford University Press, 1999. [20] Roberts, G.O. and Stramer, O., “On inference for partially observed nonlinear diffusion models using the Metropolis-Hasting algorithm”, Biometrika, 2001, 88, 603-621. [21] Schurz, H., Numerical Analysis de Stochastic Differential Equations without tears, 2002, In Handbook of Stochastic Analysis and Applications, D. Kannan, V. Lakahmikamtham (Eds.) Marcell Dekker, 237-358. [22] Shiryaev A.N., Essential of stochastic finance: facts, models, theory, World Scientific, 1999. [23] Shoji, I. and Ozaki, T., “Comparative study of estimation methods for continuous time stochastic processes”, Journal of Time Series Analysis, 1997, 18, 485-506. [24] Shoji, I. and Ozaki, T., “A statistical method of estimation and simulation for systems of stochastic differential equations”, Biometrika, 1998, 85, 240-243. [25] Shoji, I. and Ozaki, T., “Estimation for nonlinear stochastic differential equations by a local linearization method”, Stochast. Anal. Appl., 1998, 16, 733-752. [26] Sidje, R. B., “EXPOKIT: software package for computing matrix exponentials”, AMC Trans. Math. Software, 1998, 24, 130-156.

17

[27] Singer, H., “Parameter estimation of nonlinear stochastic differential equations: Simulated maximum likelihood versus extended Kalman filter and Ito-Taylor expansion, J. Comput. Graph. Stats., 2002, 11(4), 972-995. [28] Stramer, O., “The local linearization scheme for nonlinear diffusion models with discontinuous coefficients”, Stat. Prob. Lett., 1999, 42, 249-256. [29] Stramer, O. and Tweedie R. L., “Langevin-type models I: diffussion with given stationary distributions and their discretizations”, Meth. Comput. Appl. Prob., 1999, 1, 283-306. [30] Sun, J., “Rounding-Error and Perturbation Bounds for the Cholesky and LDL| factorizations”, Linear Algebra and its Applications, 1992, 173,77-97. [31] Talay, D., “Second order discretizations schemes of stochastic differential systems for the computation of the invariant law ”, Stochastics and Stochastic Reports, 1990, 29, 13—36. [32] Talay, D. and Tubaro, L. "Expansion of the global error for numerical schemes solving stochastic differential equations", Stochast. Anal. Appl., 1990, 8(4), 94—120. [33] Van Loan, C. F., “Computing integrals involving the matirx exponential”, IEEE Transactions on Automatic Control, 1978, 23, 395-404.

18

Figures Caption: Figure 1: Estimated values for the errors e(δ) obtained by the order β = 1, 2 LL schemes applied to equation (28) with different step sizes δ. Figure 2: Comparison of different numerical schemes applied to equation (29) with step size δ = 2−4 .

19

-6

-8

log2(e(h))

-10

-12

-14

-16

-18 -8

-7.5

-7

-6.5

-6

-5.5

log2(h)

-5

-4.5

-4

-3.5

-3

Explicit Extrapolated Euler

Order 3 Explicit Taylor

2

2

1

1

0

0

-1

-1

-2

0

2

4

6

8

10

-2

0

2

Order 1 LL 2

1

1

0

0

-1

-1

0

2

4

6

6

8

10

8

10

Implicit Euler

2

-2

4

8

10

-2

0

2

4

6