- Email: [email protected]
A note on convergence rate of a linearization method for the discretization of stochastic differential equations
Commun Nonlinear Sci Numer Simulat 16 (2011) 2667–2671
Contents lists available at ScienceDirect
Commun Nonlinear Sci Numer Simulat journal homepage: www.elsevier.com/locate/cnsns
Short communication
A note on convergence rate of a linearization method for the discretization of stochastic differential equations Isao Shoji Graduate School of Systems and Information Engineering, University of Tsukuba, Tsukuba Ibaraki 305-8573, Japan
a r t i c l e
i n f o
Article history: Received 12 May 2010 Received in revised form 7 September 2010 Accepted 8 September 2010 Available online 17 September 2010 Keywords: Stochastic differential equations Local linearization method Euler method Rate of convergence
a b s t r a c t This note discusses convergence rate of a linearization method for the discretization of stochastic differential equations with multiplicative noise. The method is to approximate the drift coefficient by the local linearization method and the diffusion coefficient by the Euler method. The mixed method guarantees the approximated process converges to the original one with the rate of convergence Dt, where Dt is the time interval of discretization. Ó 2010 Elsevier B.V. All rights reserved.
1. Introduction Diffusion processes are widely used for modeling time evolution of dynamic phenomena in engineering, physics, economics, and neuroscience. Since the dynamics of a diffusion process is formulated by its stochastic differential equation (SDE), how to implement SDEs is essential for those applications. Apart from simple cases, most SDEs used for real applications are nonlinear in drift and thus the manipulation is quite intractable because those SDEs can not be solved explicitly. Therefore, approximate methods of discretization are often used to get sample paths of a diffusion, or another method, called the Exact Algorithm, is also used recently; detail discussion about numerical comparison can be seen in Pelcuhetti and Roberts [10] for example. Among them, the local linearization (LL) method is easily manipulated and shows better performance than the Euler method, as reported by Hurn et al. [5] and Iacus [6] for example. Considering the LL method is derived basically from the first order expansion differently from the Euler method, the LL method is thought to be more efficient in approximation. The basic idea of the LL method is to convert a target nonlinear SDE into an approximate linear SDE just like the linearization around a critical point of a nonlinear dynamical system. Linear SDEs have their explicit solutions by which numerical experiments can be conducted easily. Therefore the LL method is used not only for generating sample paths of nonlinear SDEs but also for samplers to get their intractable probability densities numerically; for example, when making a Bayesian inference on a diffusion or filtering the unobservable components of a diffusion, its transition density, which is almost always computationally intractable, is needed and thus such a sampler is used as shown by Stramer and Tweedie [15,16]. The LL method is basically developed for SDEs with the constant diffusion coefficient and so is not directly applied to SDEs with multiplicative noise. Indeed, those SDEs can be transformed into the SDE with the constant diffusion by the change of variables, but quite a few SDEs are still difficult to transform. In that case, while we could apply the LL method to the diffusion coefficient, unexpectedly however, the resulting SDE becomes numerically intractable. Instead, as in Roberts and
E-mail address: [email protected] 1007-5704/$ - see front matter Ó 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.cnsns.2010.09.008
2668
I. Shoji / Commun Nonlinear Sci Numer Simulat 16 (2011) 2667–2671
Stramer [11], we can apply the Euler method to the diffusion coefficient while applying the LL method to the drift coefficient. This combination of the two methods, simply called the LL/Euler method, produces a linear SDE locally, thereby we can get its exact discretization as easily as the LL method. Although the LL/Euler method has a merit of numerical tractability, its convergence rate is still unknown. So, it may be interesting to ask at what speed this mixed method converges to the true process when time interval for discretization going to zero. The convergence of the LL method is discussed in Carbonell et al. [2], Jimenez and Biscay [8], and Shoji [12], where the diffusion coefficient is assumed to be constant or a deterministic function of time. But, considering applications to real problems, many interesting models, which are handled by Chan et al. [3], Heston [4] and Jones [7] in financial economics for example, are formulated by SDEs with multiplicative noise. Thus, the LL/Euler method is thought as one of practical methods for those applications. The aim of this note is to show the convergence rate of the LL/Euler method. 2. Convergence rate of the LL/Euler method We assume a process Xt is bounded by M (0 < M < 1) and satisfies the following SDE starting at X0 = x,
dX t ¼ lðX t Þdt þ rðX t ÞdBt ;
ð1Þ
where Bt is the standard Brownian motion, l is twice continuously differentiable on I = [M, M] (compact set) and r is continuous on I with satisfying the Lipschitz condition. In general, a diffusion process is not always bounded, but we can take the stopped process defined by Xt^T for the stopping time T = inf{t P 0kXtj P M}. The stopped process is bounded and is thought to have almost no problem for practical applications when taking a large M. Hence, we focus on bounded diffusions in this paper. Let time s be arbitrarily fixed and divide it equally into n-subintervals by 0 = t0 < t1 < < tn = s. For each time interval [tk, tk+1), we approximate the original SDE by applying the LL method to the drift coefficient and alternatively the Euler method to the diffusion coefficient as follows:
et ¼ X et þ X kþ1 k
Z
t kþ1
l~ ð Xe u Þdu þ
tk
Z
t kþ1
r~ ð Xe u ÞdBu ;
ð2Þ
tk
l~ ð Xe u Þ ¼ lð Xe tk Þ þ l0 ð Xe tk Þð Xe u Xe tk Þ þ
r~ ð Xe u Þ2 2
l00 ð Xe tk Þðu tk Þ ðLL methodÞ;
r~ ð Xe u Þ ¼ rð Xe tk Þ ðEuler methodÞ e t ¼ X t ¼ x. Here l0 and l00 mean the first and second derivatives of l. See Shoji and Ozaki [13] for example. Another with X 0 0 types of linearization are also presented by Biscay et al. [1], where the second derivative in (2) is dropped off, and Stramer [14], where SDEs with discontinuous coefficients can be handled. Throughout this paper, however, we consider the linearization given by (2). e u and r ~ of (2) is locally linear in X ~ is locally constant, we get the following exact discretization, Since l
n o n o e 2 e t þ l0 ð X e t Þ1 expðl0 ð X e t ÞDtÞ 1 lð X e t Þ þ l0 ð X e t Þ2 expðl0 ð X e t ÞDtÞ 1 l0 ð X e t ÞDt rð X tk Þ l00 ð X et Þ et ¼ X X kþ1 k k k k k k k k 2 Z tk þDt et Þ e t Þðt k þ Dt uÞ dBu ; þ rð X exp l0 ð X k k
ð3Þ
tk
where Dt = tk+1 tk = s/n. Thanks to the property of the Brownian motion Bu, the last stochastic integration follows the normal distribution with mean zero and variance,
rð Xe tk Þ2
e t ÞDtÞ 1 expð2l0 ð X k : et Þ 0 2l ð X
ð4Þ
k
e t as above, we eventually arrive at X e t , or X e s . If necessary, replace X e t with M X e t =j X e t j to make X e Recursively updating X n k k k k e bounded at discrete times whenever j X tk j > M. Under these settings, we can get the following claim. Proposition 1. Suppose l and r in (1) belong to C2-class and C1-class on I, respectively. Let time s be arbitrarily fixed and s = nDt. e given by (3), For the approximate process X
e s X s j2 ¼ 0; lim Ej X
Dt!0
where X means the original process satisfying (1). Proof. In the following, t and s is used as tk+1 and tk for simplicity. Define eu on u 2 [s, t) by
eu ¼ Xe u X u :
ð5Þ
2669
I. Shoji / Commun Nonlinear Sci Numer Simulat 16 (2011) 2667–2671
Clearly,
e u dX u ¼ ðl e u Þ lðX u ÞÞdu þ ðr e u Þ rðX u ÞÞdBu : ~ ðX ~ ðX deu ¼ d X Applying the Itô’s formula to
ð6Þ
2 u,
e we get,
~ u lu Þdu þ 2eu ðr ~ u ru ÞdBu þ ðr ~ u ru Þ2 du; de ¼ 2eu deu þ dheiu ¼ 2eu ðl 2 u
e u Þ and lu = l(Xu) and similarly r ~u ¼ l ~ ðX ~ u and ru. Therefore, the expectation conditional of Fs , simply denoted by meaning l Es[], is,
Es
e2t e2s ¼ 2Es
Z s
t
eu ðl~ u lu Þdu þ Es
Z
t
~ u ru Þ2 du : ðr
ð7Þ
s
In the first place, we evaluate the first term of the right-hand side of (7) and then the second term. Application of the second order Taylor expansion of l around Xs produces,
1 2
lðX u Þ ¼ lðX s Þ þ l0 ðX s ÞðX u X s Þ þ l00 ðX h ÞðX u X s Þ2 ; where Xh = (1 h)Xs + hXu for some h 2 [0, 1] which of course depends on Xu as well as Xs. Hence,
l~ u lu ¼ lð Xe s Þ lðX s Þþ l0 ð Xe s Þð Xe u Xe s Þþ
rð Xe s Þ2 2
1 2 rð Xe s Þ2
l00 ð Xe s ÞðusÞ l0 ðX s ÞðX u X s Þ l00 ðX h ÞðX u X s Þ2
e s Þeu þ lð X e s Þ lðX s Þ l0 ð X e s Þð X e s X s Þþðl0 ð X e s Þ l0 ðX s ÞÞðX u X s Þþ ¼ l0 ð X
2
1 2
l00 ð Xe s ÞðusÞ l00 ðX h ÞðX u X s Þ2
e s Þeu þAu ; ¼ l0 ð X
ð8Þ 2
2
where Au means the latter part of the right-hand side of the second equation. Since 2ab 6 a + b ,
e s Þe2 þ 2eu Au 6 2jl0 ð X e s Þje2 þ e2 þ A2 ¼ ð1 þ 2jl0 ð X e s ÞjÞe2 þ A2 : ~ u lu Þ ¼ 2l0 ð X 2eu ðl u u u u u u Rt Next, we try to evaluate Es ½ s Au du. Using the Schwartz inequality, ( ) 4 2 e s Þ lðX s ÞÞ2 þ l0 ð X e s Þ2 ð X e s X s Þ2 þ ðl0 ð X e s Þ l0 ðX s ÞÞ2 ðX u X s Þ2 þ rðX s Þ l00 ð X e s Þ2 ðu sÞ2 þ 1 l00 ðX h Þ2 ðX u X s Þ4 Au 6 5 ðlð X 4 4 ( ) 4 e s Þ lðX s ÞÞ2 þ l0 ð X e s Þ2 ð X e s X s Þ2 þ 1 ðl0 ð X e s Þ l0 ðX s ÞÞ4 þ 1 ðX u X s Þ4 þ rðX s Þ l00 ð X e s Þ2 ðu sÞ2 þ 1 l00 ðX h Þ2 ðX u X s Þ4 6 5 ðlð X 2 2 4 4 e s Þ lðX s ÞÞ2 þ l0 ð X e s Þ2 ð X e s X s Þ2 þ 1 ðjl0 ð X e s Þj þ jl0 ðX s ÞjÞ2 ðl0 ð X e s Þ l0 ðX s ÞÞ2 þ 1 ðX u X s Þ4 6 5 ðlð X 2 2 ) 4 rðX s Þ 00 e 2 1 þ l ð X s Þ ðu sÞ2 þ l00 ðX h Þ2 ðX u X s Þ4 : 4 4 From the first to the second inequality, 2ab 6 a2 + b2 is used. Since l and l0 , which are defined on the compact set I, hold the Lipschitz condition,
e s Þ lðX s ÞÞ2 6 K 1 ð X e s X s Þ2 ; ð lð X e s Þ l0 ðX s ÞÞ2 6 K 1 ð X e s X s Þ2 ; ðl0 ð X e s are bounded, for some positive number K1. Additionally, since l0 , l00 and r are continuous and Xs, Xu and X l0 ð Xe s Þ; l0 ðX s Þ; rðX s Þ; l00 ð Xe s Þ, and l00 (Xh) are all bounded. Thus, without loss of generality, we can assume,
e s X s Þ2 þ ðu sÞ2 þ ðX u X s Þ4 g: A2u 6 K 1 fð X It is well known Es[(Xu Xs)2m] = O((u s)m) for any natural number m, and thus for some positive number K2,
e s X s Þ2 þ ðu sÞ2 þ Es ½ðX u X s Þ4 g 6 K 2 fe2 þ 2ðu sÞ2 g: Es ½A2u 6 K 1 fð X s Hence,
Es
Z
t s
n o A2u du 6 K 2 Dte2s þ ð2=3ÞðDtÞ3 :
ð9Þ
Taking some positive number K, we can assume,
2Es
Z s
t
eu ðl~ u lu Þdu 6 KEs
Z s
t
e2u du þ K Dte2s þ KðDtÞ3 :
ð10Þ
2670
I. Shoji / Commun Nonlinear Sci Numer Simulat 16 (2011) 2667–2671
On the other hand, using the Lipschitz condition of r,
~ u ru Þ 2 ¼ ðr
n
rð Xe s Þ rðX s Þ ðrðX u Þ rðX s ÞÞ
o2
2 e s Þ rðX s Þ þ ðrðX u Þ rðX s ÞÞ2 6 2 rð X
n o 6 L e2s þ ðX u X s Þ2 ;
ð11Þ
for some positive number L. Therefore, using the same argument on (Xu Xs)2 as above and rearranging L, we get,
Es
Z
t
s
~ u ru Þ2 du 6 LDt e2s þ LðDtÞ2 : ðr
ð12Þ
Using (10) and (12) and applying the unconditional expectation,
E½e2t e2s 6 bDtE½e2s þ bðDtÞ2 þ bE
Z
t s
e2u du ;
for some positive number b. Here define /(t) by /ðtÞ ¼ E½e2t and then,
/ðtÞ 6 aðtÞ þ b
Z
t
ð13Þ
/ðuÞdu: s
where a(t) = (1 + b(t s))/(s) + b(t s)2. From the Gronwall inequality,
/ðtÞ 6 aðtÞ þ b
Z
t
expðbðt sÞÞaðuÞdu
s
Z ¼ 1 þ bDt þ b
t
Z expðbðt uÞÞð1 þ bðu sÞÞdu /ðsÞ þ bðDtÞ2 þ b
s
¼
Z 1 þ bDt þ b expðbDtÞ
t
expðbðt uÞÞbðu sÞ2 du
s
Dt
Z expðbuÞð1 þ buÞÞdu /ðsÞ þ bðDtÞ2 þ b expðbDtÞ
0
Dt
expðbuÞbu2 du:
0
Here note,
Z
Dt
expðbuÞð1 þ buÞÞdu ¼ OðDtÞ;
0
Z
Dt
expðbuÞbu2 du ¼ OððDtÞ2 Þ;
0
since application of the l’Hôpital’s rule produces,
R Dt lim
Dt!0
lim
Dt!0
0
R Dt 0
expð2uÞð1 þ buÞÞdu ¼ 1; Dt expð2uÞbu2 du ¼ 0: ðDtÞ2
Without loss of generality, we can assume,
/ðtÞ 6 ð1 þ bDtÞ/ðsÞ þ bðDtÞ2 :
ð14Þ
According to lemma 1.3 in Milstein [9], we get,
/ðt n Þ 6 ðexpðbsÞ 1ÞDt:
ð15Þ
e t ¼ X t ¼ x by definition and thus /(t0) = 0. Therefore, Here, note that X 0 0
lim /ðt n Þ ¼ lim /ðsÞ ¼ 0:
Dt!0
Dt!0
This completes the proof. h Theorem 1. The rate of convergence in Proposition 1 is Dt. Proof. Immediately from (15). Zhao et al. [17] discusses the convergence rate of the other method when the diffusion coefficient is a time-homogeneous function of a diffusion process and it shows the rate is the same as the LL/Euler method. It can be seen in the above proof that the determinant of the convergence rate is the rate of the last term of the right-hand side of (14). This rate comes from the lower rate of K(D t)3 in (10) or L(Dt)2 in (12). In other words, the convergence rate becomes worse because of the
I. Shoji / Commun Nonlinear Sci Numer Simulat 16 (2011) 2667–2671
2671
approximation of the diffusion coefficient by the Euler method, that is L(D t)2 in (12). So, it may be possible to improve the convergence rate if a more efficient approximation method of the diffusion coefficient is used. For example, we could apply the similar linear approximation to r and then, as shown in (8), the difference in (11) can be formulated as,
~ u ru Þ 2 ¼ ðr
r0 ð Xe s Þeu þ C u
2
62
r0 ð Xe s Þ2 e2u þ C 2u ;
where,
Es
Z
t s
n o C 2u du 6 L Dte2s þ ð2=3ÞðDtÞ3 ;
for some positive number L. See the inequality (9). Therefore, the convergence rate would be improved to (Dt)2. This is the same rate as in the original setting of the LL method that the diffusion coefficient is constant. h Furthermore, we could apply higher order Taylor expansions to l and r, but the convergence rate would not be improved because the term jXu Xsj4 would be left, thereby the rate remains O((Dt)2). In addition, from a practical point of view, applying a higher order method to r as well as l may not be promising because such a method would produce a more complicated approximation of the diffusion coefficient so that the discretization is not so easy as the LL/Euler method and thus be quite intractable for manipulation. 3. Concluding remarks This note shows the mixed local linearization method guarantees the approximated process converges to the true one as time interval for discretization Dt going to zero with the rate of convergence Dt in the mean squared error. It might be possible to improve the convergence rate by using a more efficient approximation of the diffusion coefficient which probably would produce a more complicated approximation so that the discretization might not be so easy as the LL/Euler method although the higher rate, or (Dt)2, can be achieved in the original setting of the LL method. This note handles one-dimensional diffusions, but the same argument can be applied to multi-dimensional diffusions with somewhat cumbersome notations. Acknowledgements We are grateful to the editor and to the anonymous referees for their helpful comments and suggestions. References [1] Biscay RJ, Jimenez JC, Riera J, Valdes P. Local linearization method for the numerical solution of stochastic differential equations. Ann Inst Statist Math 1996;48:631–44. [2] Carbonell F, Jimenez JC, Biscay RJ. Weak local linear discretizations for stochastic differential equations: convergence and numerical schemes. J Comput Appl Math 2006;197:578–96. [3] Chan KC, Karolyi GA, Longstaff FA, Sanders AS. An empirical comparison of alternative models of the short-term interest rate. J Financ 1992;47:1209–27. [4] Heston SL. A closed-form solution for options with stochastic volatility with applications to bonds and currency options. Rev Financ Stud 1993;6:327–43. [5] Hurn AS, Jeisman JI, Lindsay KA. Seeing the wood for the trees: a critical evaluation of methods to estimate the parameters of stochastic differential equations. J Financ Econom 2007;5:390–455. [6] Iacus SM. Simulation and inference for stochastic differential equations: with r examples. Springer Series in Statistics. Springer; 2008. [7] Jones CS. The dynamics of stochastic volatility: evidence from underlying and options markets. J Econom 2003;116:181–224. [8] Jimenez JC, Biscay RJ. Approximation of continuous time stochastic processes by the local linearization method revisited. Stoch Anal Appl 2002;20:105–21. [9] Milstein GN. Numerical Integration of Stochastic Differential Equations. Kluwer Academic Publishers; 1988. [10] Pelcuhetti S, Roberts GO. An empirical study of the efficiency of EA for diffusion simulation, Working Paper No.08-14, Center for Research in Statistical Methodology (CRiSM), The University of Warwick, 2008. [11] Roberts GO, Stramer O. Langevin diffusions and metropolis-hastings algorithms. Meth Comput Appl Probab 2002;4:337–57. [12] Shoji I. Approximation of continuous time stochastic processes by a local linearization method. Math Comput 1998;67:287–98. [13] Shoji I, Ozaki T. A statistical method of estimation and simulation for systems of stochastic differential equations. Biometrika 1998;85:240–3. [14] Stramer O. The local linearization scheme for nonlinear diffusion models with discontinuous coefficients. Statist Probab Lett 1999;42:249–56. [15] Stramer O, Tweedie RL. Langevin-type models I: diffusions with given stationary distributions, and their discretizations. Meth Comput Appl Probab 1999;1:283–306. [16] Stramer O, Tweedie RL. Langevin-type models II: self-targeting candidates for MCMC algorithms. Meth Comput Appl Probab 1999;1:307–28. [17] Zhao W, Tian L, Ju L. Convergence analysis of a splitting method for stochastic differential equations. Int J Numer Anal Model 2008;5:673–92.
Contents lists available at ScienceDirect
Commun Nonlinear Sci Numer Simulat journal homepage: www.elsevier.com/locate/cnsns
Short communication
A note on convergence rate of a linearization method for the discretization of stochastic differential equations Isao Shoji Graduate School of Systems and Information Engineering, University of Tsukuba, Tsukuba Ibaraki 305-8573, Japan
a r t i c l e
i n f o
Article history: Received 12 May 2010 Received in revised form 7 September 2010 Accepted 8 September 2010 Available online 17 September 2010 Keywords: Stochastic differential equations Local linearization method Euler method Rate of convergence
a b s t r a c t This note discusses convergence rate of a linearization method for the discretization of stochastic differential equations with multiplicative noise. The method is to approximate the drift coefficient by the local linearization method and the diffusion coefficient by the Euler method. The mixed method guarantees the approximated process converges to the original one with the rate of convergence Dt, where Dt is the time interval of discretization. Ó 2010 Elsevier B.V. All rights reserved.
1. Introduction Diffusion processes are widely used for modeling time evolution of dynamic phenomena in engineering, physics, economics, and neuroscience. Since the dynamics of a diffusion process is formulated by its stochastic differential equation (SDE), how to implement SDEs is essential for those applications. Apart from simple cases, most SDEs used for real applications are nonlinear in drift and thus the manipulation is quite intractable because those SDEs can not be solved explicitly. Therefore, approximate methods of discretization are often used to get sample paths of a diffusion, or another method, called the Exact Algorithm, is also used recently; detail discussion about numerical comparison can be seen in Pelcuhetti and Roberts [10] for example. Among them, the local linearization (LL) method is easily manipulated and shows better performance than the Euler method, as reported by Hurn et al. [5] and Iacus [6] for example. Considering the LL method is derived basically from the first order expansion differently from the Euler method, the LL method is thought to be more efficient in approximation. The basic idea of the LL method is to convert a target nonlinear SDE into an approximate linear SDE just like the linearization around a critical point of a nonlinear dynamical system. Linear SDEs have their explicit solutions by which numerical experiments can be conducted easily. Therefore the LL method is used not only for generating sample paths of nonlinear SDEs but also for samplers to get their intractable probability densities numerically; for example, when making a Bayesian inference on a diffusion or filtering the unobservable components of a diffusion, its transition density, which is almost always computationally intractable, is needed and thus such a sampler is used as shown by Stramer and Tweedie [15,16]. The LL method is basically developed for SDEs with the constant diffusion coefficient and so is not directly applied to SDEs with multiplicative noise. Indeed, those SDEs can be transformed into the SDE with the constant diffusion by the change of variables, but quite a few SDEs are still difficult to transform. In that case, while we could apply the LL method to the diffusion coefficient, unexpectedly however, the resulting SDE becomes numerically intractable. Instead, as in Roberts and
E-mail address: [email protected] 1007-5704/$ - see front matter Ó 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.cnsns.2010.09.008
2668
I. Shoji / Commun Nonlinear Sci Numer Simulat 16 (2011) 2667–2671
Stramer [11], we can apply the Euler method to the diffusion coefficient while applying the LL method to the drift coefficient. This combination of the two methods, simply called the LL/Euler method, produces a linear SDE locally, thereby we can get its exact discretization as easily as the LL method. Although the LL/Euler method has a merit of numerical tractability, its convergence rate is still unknown. So, it may be interesting to ask at what speed this mixed method converges to the true process when time interval for discretization going to zero. The convergence of the LL method is discussed in Carbonell et al. [2], Jimenez and Biscay [8], and Shoji [12], where the diffusion coefficient is assumed to be constant or a deterministic function of time. But, considering applications to real problems, many interesting models, which are handled by Chan et al. [3], Heston [4] and Jones [7] in financial economics for example, are formulated by SDEs with multiplicative noise. Thus, the LL/Euler method is thought as one of practical methods for those applications. The aim of this note is to show the convergence rate of the LL/Euler method. 2. Convergence rate of the LL/Euler method We assume a process Xt is bounded by M (0 < M < 1) and satisfies the following SDE starting at X0 = x,
dX t ¼ lðX t Þdt þ rðX t ÞdBt ;
ð1Þ
where Bt is the standard Brownian motion, l is twice continuously differentiable on I = [M, M] (compact set) and r is continuous on I with satisfying the Lipschitz condition. In general, a diffusion process is not always bounded, but we can take the stopped process defined by Xt^T for the stopping time T = inf{t P 0kXtj P M}. The stopped process is bounded and is thought to have almost no problem for practical applications when taking a large M. Hence, we focus on bounded diffusions in this paper. Let time s be arbitrarily fixed and divide it equally into n-subintervals by 0 = t0 < t1 < < tn = s. For each time interval [tk, tk+1), we approximate the original SDE by applying the LL method to the drift coefficient and alternatively the Euler method to the diffusion coefficient as follows:
et ¼ X et þ X kþ1 k
Z
t kþ1
l~ ð Xe u Þdu þ
tk
Z
t kþ1
r~ ð Xe u ÞdBu ;
ð2Þ
tk
l~ ð Xe u Þ ¼ lð Xe tk Þ þ l0 ð Xe tk Þð Xe u Xe tk Þ þ
r~ ð Xe u Þ2 2
l00 ð Xe tk Þðu tk Þ ðLL methodÞ;
r~ ð Xe u Þ ¼ rð Xe tk Þ ðEuler methodÞ e t ¼ X t ¼ x. Here l0 and l00 mean the first and second derivatives of l. See Shoji and Ozaki [13] for example. Another with X 0 0 types of linearization are also presented by Biscay et al. [1], where the second derivative in (2) is dropped off, and Stramer [14], where SDEs with discontinuous coefficients can be handled. Throughout this paper, however, we consider the linearization given by (2). e u and r ~ of (2) is locally linear in X ~ is locally constant, we get the following exact discretization, Since l
n o n o e 2 e t þ l0 ð X e t Þ1 expðl0 ð X e t ÞDtÞ 1 lð X e t Þ þ l0 ð X e t Þ2 expðl0 ð X e t ÞDtÞ 1 l0 ð X e t ÞDt rð X tk Þ l00 ð X et Þ et ¼ X X kþ1 k k k k k k k k 2 Z tk þDt et Þ e t Þðt k þ Dt uÞ dBu ; þ rð X exp l0 ð X k k
ð3Þ
tk
where Dt = tk+1 tk = s/n. Thanks to the property of the Brownian motion Bu, the last stochastic integration follows the normal distribution with mean zero and variance,
rð Xe tk Þ2
e t ÞDtÞ 1 expð2l0 ð X k : et Þ 0 2l ð X
ð4Þ
k
e t as above, we eventually arrive at X e t , or X e s . If necessary, replace X e t with M X e t =j X e t j to make X e Recursively updating X n k k k k e bounded at discrete times whenever j X tk j > M. Under these settings, we can get the following claim. Proposition 1. Suppose l and r in (1) belong to C2-class and C1-class on I, respectively. Let time s be arbitrarily fixed and s = nDt. e given by (3), For the approximate process X
e s X s j2 ¼ 0; lim Ej X
Dt!0
where X means the original process satisfying (1). Proof. In the following, t and s is used as tk+1 and tk for simplicity. Define eu on u 2 [s, t) by
eu ¼ Xe u X u :
ð5Þ
2669
I. Shoji / Commun Nonlinear Sci Numer Simulat 16 (2011) 2667–2671
Clearly,
e u dX u ¼ ðl e u Þ lðX u ÞÞdu þ ðr e u Þ rðX u ÞÞdBu : ~ ðX ~ ðX deu ¼ d X Applying the Itô’s formula to
ð6Þ
2 u,
e we get,
~ u lu Þdu þ 2eu ðr ~ u ru ÞdBu þ ðr ~ u ru Þ2 du; de ¼ 2eu deu þ dheiu ¼ 2eu ðl 2 u
e u Þ and lu = l(Xu) and similarly r ~u ¼ l ~ ðX ~ u and ru. Therefore, the expectation conditional of Fs , simply denoted by meaning l Es[], is,
Es
e2t e2s ¼ 2Es
Z s
t
eu ðl~ u lu Þdu þ Es
Z
t
~ u ru Þ2 du : ðr
ð7Þ
s
In the first place, we evaluate the first term of the right-hand side of (7) and then the second term. Application of the second order Taylor expansion of l around Xs produces,
1 2
lðX u Þ ¼ lðX s Þ þ l0 ðX s ÞðX u X s Þ þ l00 ðX h ÞðX u X s Þ2 ; where Xh = (1 h)Xs + hXu for some h 2 [0, 1] which of course depends on Xu as well as Xs. Hence,
l~ u lu ¼ lð Xe s Þ lðX s Þþ l0 ð Xe s Þð Xe u Xe s Þþ
rð Xe s Þ2 2
1 2 rð Xe s Þ2
l00 ð Xe s ÞðusÞ l0 ðX s ÞðX u X s Þ l00 ðX h ÞðX u X s Þ2
e s Þeu þ lð X e s Þ lðX s Þ l0 ð X e s Þð X e s X s Þþðl0 ð X e s Þ l0 ðX s ÞÞðX u X s Þþ ¼ l0 ð X
2
1 2
l00 ð Xe s ÞðusÞ l00 ðX h ÞðX u X s Þ2
e s Þeu þAu ; ¼ l0 ð X
ð8Þ 2
2
where Au means the latter part of the right-hand side of the second equation. Since 2ab 6 a + b ,
e s Þe2 þ 2eu Au 6 2jl0 ð X e s Þje2 þ e2 þ A2 ¼ ð1 þ 2jl0 ð X e s ÞjÞe2 þ A2 : ~ u lu Þ ¼ 2l0 ð X 2eu ðl u u u u u u Rt Next, we try to evaluate Es ½ s Au du. Using the Schwartz inequality, ( ) 4 2 e s Þ lðX s ÞÞ2 þ l0 ð X e s Þ2 ð X e s X s Þ2 þ ðl0 ð X e s Þ l0 ðX s ÞÞ2 ðX u X s Þ2 þ rðX s Þ l00 ð X e s Þ2 ðu sÞ2 þ 1 l00 ðX h Þ2 ðX u X s Þ4 Au 6 5 ðlð X 4 4 ( ) 4 e s Þ lðX s ÞÞ2 þ l0 ð X e s Þ2 ð X e s X s Þ2 þ 1 ðl0 ð X e s Þ l0 ðX s ÞÞ4 þ 1 ðX u X s Þ4 þ rðX s Þ l00 ð X e s Þ2 ðu sÞ2 þ 1 l00 ðX h Þ2 ðX u X s Þ4 6 5 ðlð X 2 2 4 4 e s Þ lðX s ÞÞ2 þ l0 ð X e s Þ2 ð X e s X s Þ2 þ 1 ðjl0 ð X e s Þj þ jl0 ðX s ÞjÞ2 ðl0 ð X e s Þ l0 ðX s ÞÞ2 þ 1 ðX u X s Þ4 6 5 ðlð X 2 2 ) 4 rðX s Þ 00 e 2 1 þ l ð X s Þ ðu sÞ2 þ l00 ðX h Þ2 ðX u X s Þ4 : 4 4 From the first to the second inequality, 2ab 6 a2 + b2 is used. Since l and l0 , which are defined on the compact set I, hold the Lipschitz condition,
e s Þ lðX s ÞÞ2 6 K 1 ð X e s X s Þ2 ; ð lð X e s Þ l0 ðX s ÞÞ2 6 K 1 ð X e s X s Þ2 ; ðl0 ð X e s are bounded, for some positive number K1. Additionally, since l0 , l00 and r are continuous and Xs, Xu and X l0 ð Xe s Þ; l0 ðX s Þ; rðX s Þ; l00 ð Xe s Þ, and l00 (Xh) are all bounded. Thus, without loss of generality, we can assume,
e s X s Þ2 þ ðu sÞ2 þ ðX u X s Þ4 g: A2u 6 K 1 fð X It is well known Es[(Xu Xs)2m] = O((u s)m) for any natural number m, and thus for some positive number K2,
e s X s Þ2 þ ðu sÞ2 þ Es ½ðX u X s Þ4 g 6 K 2 fe2 þ 2ðu sÞ2 g: Es ½A2u 6 K 1 fð X s Hence,
Es
Z
t s
n o A2u du 6 K 2 Dte2s þ ð2=3ÞðDtÞ3 :
ð9Þ
Taking some positive number K, we can assume,
2Es
Z s
t
eu ðl~ u lu Þdu 6 KEs
Z s
t
e2u du þ K Dte2s þ KðDtÞ3 :
ð10Þ
2670
I. Shoji / Commun Nonlinear Sci Numer Simulat 16 (2011) 2667–2671
On the other hand, using the Lipschitz condition of r,
~ u ru Þ 2 ¼ ðr
n
rð Xe s Þ rðX s Þ ðrðX u Þ rðX s ÞÞ
o2
2 e s Þ rðX s Þ þ ðrðX u Þ rðX s ÞÞ2 6 2 rð X
n o 6 L e2s þ ðX u X s Þ2 ;
ð11Þ
for some positive number L. Therefore, using the same argument on (Xu Xs)2 as above and rearranging L, we get,
Es
Z
t
s
~ u ru Þ2 du 6 LDt e2s þ LðDtÞ2 : ðr
ð12Þ
Using (10) and (12) and applying the unconditional expectation,
E½e2t e2s 6 bDtE½e2s þ bðDtÞ2 þ bE
Z
t s
e2u du ;
for some positive number b. Here define /(t) by /ðtÞ ¼ E½e2t and then,
/ðtÞ 6 aðtÞ þ b
Z
t
ð13Þ
/ðuÞdu: s
where a(t) = (1 + b(t s))/(s) + b(t s)2. From the Gronwall inequality,
/ðtÞ 6 aðtÞ þ b
Z
t
expðbðt sÞÞaðuÞdu
s
Z ¼ 1 þ bDt þ b
t
Z expðbðt uÞÞð1 þ bðu sÞÞdu /ðsÞ þ bðDtÞ2 þ b
s
¼
Z 1 þ bDt þ b expðbDtÞ
t
expðbðt uÞÞbðu sÞ2 du
s
Dt
Z expðbuÞð1 þ buÞÞdu /ðsÞ þ bðDtÞ2 þ b expðbDtÞ
0
Dt
expðbuÞbu2 du:
0
Here note,
Z
Dt
expðbuÞð1 þ buÞÞdu ¼ OðDtÞ;
0
Z
Dt
expðbuÞbu2 du ¼ OððDtÞ2 Þ;
0
since application of the l’Hôpital’s rule produces,
R Dt lim
Dt!0
lim
Dt!0
0
R Dt 0
expð2uÞð1 þ buÞÞdu ¼ 1; Dt expð2uÞbu2 du ¼ 0: ðDtÞ2
Without loss of generality, we can assume,
/ðtÞ 6 ð1 þ bDtÞ/ðsÞ þ bðDtÞ2 :
ð14Þ
According to lemma 1.3 in Milstein [9], we get,
/ðt n Þ 6 ðexpðbsÞ 1ÞDt:
ð15Þ
e t ¼ X t ¼ x by definition and thus /(t0) = 0. Therefore, Here, note that X 0 0
lim /ðt n Þ ¼ lim /ðsÞ ¼ 0:
Dt!0
Dt!0
This completes the proof. h Theorem 1. The rate of convergence in Proposition 1 is Dt. Proof. Immediately from (15). Zhao et al. [17] discusses the convergence rate of the other method when the diffusion coefficient is a time-homogeneous function of a diffusion process and it shows the rate is the same as the LL/Euler method. It can be seen in the above proof that the determinant of the convergence rate is the rate of the last term of the right-hand side of (14). This rate comes from the lower rate of K(D t)3 in (10) or L(Dt)2 in (12). In other words, the convergence rate becomes worse because of the
I. Shoji / Commun Nonlinear Sci Numer Simulat 16 (2011) 2667–2671
2671
approximation of the diffusion coefficient by the Euler method, that is L(D t)2 in (12). So, it may be possible to improve the convergence rate if a more efficient approximation method of the diffusion coefficient is used. For example, we could apply the similar linear approximation to r and then, as shown in (8), the difference in (11) can be formulated as,
~ u ru Þ 2 ¼ ðr
r0 ð Xe s Þeu þ C u
2
62
r0 ð Xe s Þ2 e2u þ C 2u ;
where,
Es
Z
t s
n o C 2u du 6 L Dte2s þ ð2=3ÞðDtÞ3 ;
for some positive number L. See the inequality (9). Therefore, the convergence rate would be improved to (Dt)2. This is the same rate as in the original setting of the LL method that the diffusion coefficient is constant. h Furthermore, we could apply higher order Taylor expansions to l and r, but the convergence rate would not be improved because the term jXu Xsj4 would be left, thereby the rate remains O((Dt)2). In addition, from a practical point of view, applying a higher order method to r as well as l may not be promising because such a method would produce a more complicated approximation of the diffusion coefficient so that the discretization is not so easy as the LL/Euler method and thus be quite intractable for manipulation. 3. Concluding remarks This note shows the mixed local linearization method guarantees the approximated process converges to the true one as time interval for discretization Dt going to zero with the rate of convergence Dt in the mean squared error. It might be possible to improve the convergence rate by using a more efficient approximation of the diffusion coefficient which probably would produce a more complicated approximation so that the discretization might not be so easy as the LL/Euler method although the higher rate, or (Dt)2, can be achieved in the original setting of the LL method. This note handles one-dimensional diffusions, but the same argument can be applied to multi-dimensional diffusions with somewhat cumbersome notations. Acknowledgements We are grateful to the editor and to the anonymous referees for their helpful comments and suggestions. References [1] Biscay RJ, Jimenez JC, Riera J, Valdes P. Local linearization method for the numerical solution of stochastic differential equations. Ann Inst Statist Math 1996;48:631–44. [2] Carbonell F, Jimenez JC, Biscay RJ. Weak local linear discretizations for stochastic differential equations: convergence and numerical schemes. J Comput Appl Math 2006;197:578–96. [3] Chan KC, Karolyi GA, Longstaff FA, Sanders AS. An empirical comparison of alternative models of the short-term interest rate. J Financ 1992;47:1209–27. [4] Heston SL. A closed-form solution for options with stochastic volatility with applications to bonds and currency options. Rev Financ Stud 1993;6:327–43. [5] Hurn AS, Jeisman JI, Lindsay KA. Seeing the wood for the trees: a critical evaluation of methods to estimate the parameters of stochastic differential equations. J Financ Econom 2007;5:390–455. [6] Iacus SM. Simulation and inference for stochastic differential equations: with r examples. Springer Series in Statistics. Springer; 2008. [7] Jones CS. The dynamics of stochastic volatility: evidence from underlying and options markets. J Econom 2003;116:181–224. [8] Jimenez JC, Biscay RJ. Approximation of continuous time stochastic processes by the local linearization method revisited. Stoch Anal Appl 2002;20:105–21. [9] Milstein GN. Numerical Integration of Stochastic Differential Equations. Kluwer Academic Publishers; 1988. [10] Pelcuhetti S, Roberts GO. An empirical study of the efficiency of EA for diffusion simulation, Working Paper No.08-14, Center for Research in Statistical Methodology (CRiSM), The University of Warwick, 2008. [11] Roberts GO, Stramer O. Langevin diffusions and metropolis-hastings algorithms. Meth Comput Appl Probab 2002;4:337–57. [12] Shoji I. Approximation of continuous time stochastic processes by a local linearization method. Math Comput 1998;67:287–98. [13] Shoji I, Ozaki T. A statistical method of estimation and simulation for systems of stochastic differential equations. Biometrika 1998;85:240–3. [14] Stramer O. The local linearization scheme for nonlinear diffusion models with discontinuous coefficients. Statist Probab Lett 1999;42:249–56. [15] Stramer O, Tweedie RL. Langevin-type models I: diffusions with given stationary distributions, and their discretizations. Meth Comput Appl Probab 1999;1:283–306. [16] Stramer O, Tweedie RL. Langevin-type models II: self-targeting candidates for MCMC algorithms. Meth Comput Appl Probab 1999;1:307–28. [17] Zhao W, Tian L, Ju L. Convergence analysis of a splitting method for stochastic differential equations. Int J Numer Anal Model 2008;5:673–92.