In this paper we report ongoing work on parametric estimation for diffusions

using Monte Carlo EM algorithms. The work presented here has already been

extended to high dimensional problems and models with observation errors

Papaspiliopoulos, Omiros

Sermaidis, Giorgos

English

Warwick Research Archives Portal Repository

 University of Warwick institutional repository: http://go.warwick.ac.uk/wrap  This paper is made available online in accordance with publisher policies. Please scroll down to view the document itself. Please refer to the repository record for this item and our policy information available from the repository home page for further information.  To see the final version of this paper please visit the publisher’s website. Access to the published version may require a subscription. Author(s):  O Papaspiliopoulos and G Sermaidis Article Title: Parametric Estimation of Discretely Observed Diffusions using the EM Algorithm Year of publication: 2007 Link to published article: http://www2.warwick.ac.uk/fac/sci/statistics/crism/research/2007/paper07-21 Publisher statement:  None   1Parametric estimation of discretely observeddiffusions using the EM algorithmOmiros Papaspiliopoulos, Giorgos SermaidisDepartment of Statistics, University of Warwick, UKKeywords: EM algorithm, diffusion processes, random weightsAMS: 60J80AbstractIn this paper we report ongoing work on parametric estimation for diffusionsusing Monte Carlo EM algorithms. The work presented here has already beenextended to high dimensional problems and models with observation errors.1. IntroductionDiffusion processes have become an extensively used tool to model con-tinuous time phenomena in many scientific areas. A diffusion process V isdefined as the solution of a SDE of the formdVs = b(Vs; θ) ds+ σ(Vs; θ) dBs, (1)driven by the scalar Brownian motion B. The functionals b( · ; θ) and σ( · ; θ)are called the drift and the diffusion coefficient respectively and are allowed todepend on some parameters θ ∈ Θ. These functionals are presumed to satisfythe regularity conditions that guarantee a weakly unique, global solution ofequation (1); see chapter 4 of [1] for further details. In this paper we consideronly one-dimensional diffusions, although multivariate extensions are possible.Altough the process is defined in continuous time, the available dataare always sampled in discrete time. We assume that the process is observedwithout error at a given collection of time instancesv = {Vt0 , Vt1 , . . . , Vtn}, 0 = t0 < t1 < · · · < tnand denote the time increments between consecutive observations by ∆ti =ti − ti−1 for 1 ≤ i ≤ n.The log-likelihood of the dataset v is given byl(θ | v) =n∑i=1log{p∆ti (Vti , Vti−1 ; θ)},where pt(v, w; θ) is the transition density:pt(v, w; θ) = P (Vt ∈ dw | V0 = v; θ)/dw, t > 0, w, v ∈ R.15th European Young Statisticians Meeting 15thEYSMSeptember 10-14, 2007 Castro Urdiales (Spain)CRiSM Paper No. 07-21, www.warwick.ac.uk/go/crismEM for diffusions 2Estimation of the parameters of the diffusion process via maximum likelihood(ML) is hard since the transition density is typically not analytically available.On the other hand, ML estimation for continuous observed v is straightfor-ward since a likelihood can be easily derived (see (3)). This suggests treatingthe paths between consecutive observations as missing data and proceed toparameter estimation using the Expectation-Maximization algorithm.2. An Expectation-Maximization approachThe EM algorithm is a general purpose algorithm for maximum likeli-hood estimation in a wide variety of situations where the likelihood of the ob-served data is intractable but the joint likelihood of the observed and missingdata has a simple form (see [2]). The algorithm works for any augmentationscheme however appropriately defined missing data can lead to easier and moreefficient implementation. Without loss of generality we assume two data pointsand denote the observed data by vobs = {V0 = u, Vt = w}.The aim is to appropriately transform the missing paths so that thecomplete likelihood can be written in an explicit way. This is achieved byapplying two transformations. First, we transform Vs → η(Vs) =: Xs, whereη(u; θ) =u∫1σ(z; θ)dz .Applying Itô’s rule to the transformed process we get thatdXs = α(Xs; θ) ds+ dBs, X0 = x, s ∈ [0, t], (2)whereα(u; θ) =b{η−1(u; θ); θ}σ{η−1(u; θ); θ} − σ′{η−1(u; θ); θ} / 2, u ∈ R;η−1 denotes the inverse transformation and σ′denotes the derivative w.r.t. thespace variable. Note that the starting and ending point are functions of θ asx = η(u; θ) and y = η(w; θ). LetA(u; θ) :=∫ uα(z; θ)dzbe any anti-derivative of α and p˜t(x, y; θ) be the transition density of X. Se-cond, we apply the transformation Xs → X˙s where:X˙s := Xs − (1− st)x(θ)− sty(θ), s ∈ [0, t].X˙ is a diffusion bridge starting from X˙0 = 0 and finishing at X˙t = 0 and itsdynamics depend on θ. The inverse transformation from X˙ to X is:gθ(X˙s) := X˙s + (1− st)x(θ) +sty(θ), s ∈ [0, t].15th European Young Statisticians Meeting 15thEYSMSeptember 10-14, 2007 Castro Urdiales (Spain)CRiSM Paper No. 07-21, www.warwick.ac.uk/go/crismEM for diffusions 3The augmentation scheme for the EM algorithm is vmis = (X˙s, s ∈ [0, t]) forthe missing data and vcom = {vobs,vmis} for the complete data.LetW(t,x,y) be the standard Brownian Bridge (BB) measure and Q(t,x,y)θthe corresponding measure induced by the diffusion process (2). The conditio-nal density of X˙ on the observed data (w.r.t. the standard BB measure) isderived using Lemma 2 in [3] as:pi(X˙ | vobs, θ) = dQ(t,x,y)θdW(t,x,y){gθ(X˙s)}= Gθ{gθ(X˙s)}Gθ(ω) =Nt(y − x)p˜t(x, y; θ)exp{A(y; θ)−A(x; θ)−∫ t012(α2 + α′)(ωs; θ)ds}, (3)where ω is any path starting at x at time 0 and finishing at y at time t andNt(u) is the density of the normal distribution with mean 0 and variance tevaluated at u ∈ R. The complete log-likelihood can now be written asl(θ |vcom) = log |η′(w; θ)|+ log[Nt{y(θ)− x(θ)}] +A{y(θ); θ} −A{x(θ); θ}−∫ t012(α2 + α′){gθ(X˙s); θ}ds.It is helpful to introduce the random variable U ∼ Un[0, t], which is independentof vcom. The E-step requires analytic evaluation ofQ(θ, θ′) = Evmis|vobs;θ′ [ l(θ |vcom) ] = log |η′(w; θ)|+ logNt{y(θ)− x(θ)}+A{y(θ); θ} −A{x(θ); θ} − t2E(vmis,U)|vobs;θ′ [ (α2 + α′)(gθ(X˙U ); θ) ]. (4)The expression above cannot be evaluated analytically so we will estimate itby Monte Carlo. This suggests maximizing the estimator Q˜(θ, θ′) of Q(θ, θ′)which leads to a Monte Carlo Expectation Maximization algorithm (MCEM).2.1. Exact SimulationOne solution to the problem would be given by simulating M samplesfrom X˙U and estimating the intractable expectation with its sample average.Recent advances in simulation methodology have made exact simulation (ES)of diffusion processes feasible. The technique is called Retrospective Samplingand the resulting algorithm Exact Algorithm (EA). The simulation is based onrejection sampling; details are not presented in this paper but the reader cansee [3] for further details. So the θ(s+1) estimate of the MCEM algorithm isfound as follows:1. Simulate ui ∼ Un[0, t], for i = 1, 2, . . . ,M.2. Simulate xui ∼ Q(t,x,y)θ(s) using the EA.3. Apply the transformation in (2).15th European Young Statisticians Meeting 15thEYSMSeptember 10-14, 2007 Castro Urdiales (Spain)CRiSM Paper No. 07-21, www.warwick.ac.uk/go/crismEM for diffusions 44. Get θ(s+1) by maximizingQ˜(θ, θ(s)) = log |η′(w; θ)|+ logNt{y(θ)− x(θ)}+A{y(θ); θ}−A{x(θ); θ} − t2MM∑i=1(α2 + α′){gθ(x˙ui); θ}.2.2. Importance SamplingA drawback of estimating Q(θ, θ′) with the EA is that the latter is basedon rejection sampling which implies the simulation of a larger number of sam-ples than M that we actually use. This observation suggests to proceed to anImportance Sampling (IS) estimation for the unknown expectation in (4). Theidea is to simulate M samples from a standard Brownian Bridge and weightthem by the likelihood ratio (3). This expression cannot be calculated analyti-cally. However, following the approach by [4] we only need unbiased estimatesof the weights which can be derived using the Poisson Estimator (see [3]). Letf(X˙s; θ) =12(α2 + α′){gθ(X˙s); θ}.Let c ∈ R, λ > 0 be user-specified constans, and Ψ = {ψ1, ...ψκ} be a Poissonprocess of rate λ on [0, t]. Then it can be checked thatexp{−∫ t0f(X˙s)ds}= exp{(λ− c)t}E(wθ | X˙), wθ =κ∏j=1c− f(X˙ψj ; θ)λ.So, the θ(s+1) estimate of the MCEM algorithm using IS is derived by:1. Simulate κi ∼ Po(λt), for i = 1, 2, . . . ,M.2. Simulate ψ(i)j ∼ Un[0, t], for j = 1, 2, . . . , κi + 1.3. Simulate x˙ψ(i)j∼W(t,0,0).4. For every i use first κi samples to estimate the weight w(i)θ(s).5. Get θ(s+1) by maximizingQ˜(θ, θ(s)) = log |η′(w; θ)|+ logNt{y(θ)− x(θ)}+A{y(θ); θ}−A{x(θ); θ} −tM∑i=1(α2 + α′){gθ(x˙ψ(i)κi+1); θ}w(i)θ(s)2M∑i=1w(i)θ(s).15th European Young Statisticians Meeting 15thEYSMSeptember 10-14, 2007 Castro Urdiales (Spain)CRiSM Paper No. 07-21, www.warwick.ac.uk/go/crismEM for diffusions 53. An illustrative exampleThe algorithms are applied to an Ornstein-Uhlenbeck (OU) process forwhich b(v; θ) = −θv, σ(v; θ) = 1.0, θ ∈ (0,+∞). This process has been chosenfor illustration since the likelihood is available in this case. A sample size of1000 observations was simulated from the process using θ = 2.0 and V0 = 0.0.The algorithms were performed using M = 100 for the first five iterations andM = 1000 for the last five. The initial value was θ(0) = 1.0, the MLE was2.00274 and θ(10) was 2.00107 and 1.99918 for the ES and IS implementationsrespectively. The results are shown in Figure 1.2 4 6 8 10−8 00−7 80−7 60−7 40IterationL og −l i ke li ho od3 4 5 6 7 8 9 10−7 32 .5−7 32 .0−7 31 .5−7 31 .0−7 30 .5IterationL og −l i ke li ho odFigure 1: The log-likelihood evaluated at every EM iteration. The dotted linerepresents the ES method, the solid line the IS method and the dashed line thelog-likelihood evaluated at the true MLE. The right panel zooms at the lastseven iterations.4. Bibliography[1] Kloeden, P. and Platen, E. (1995) Numerical Solution of Stochastic Diffe-rential Equations. New York: Springer.[2] Dempster, A. P., Laird, N. M. and Rubin, D. B. (1997)Maximum likelihoodfrom incomplete data via the em algorithm (with discussion). J. R. Statist.Soc. B, 39, 1–38.[3] Beskos, A., Papaspiliopoulos, O., Roberts, G. O. and Fearnhead, P. (2006)Exact and computationally efficient likelihood-based estimation for discre-tely observed diffusion processes. J. R. Stat. Soc. Ser. B Stat. Methodol.,68, 333–382.[4] Fearnhead, P., Papaspiliopoulos, O. and Roberts, G. O. (2006) Particlefilters for partially observed diffusions. In revision.15th European Young Statisticians Meeting 15thEYSMSeptember 10-14, 2007 Castro Urdiales (Spain)CRiSM Paper No. 07-21, www.warwick.ac.uk/go/crism

Parametric estimation of discretely observed diffusions using the EM algorithm

Abstract

Similar works

Full text

Available Versions

Warwick Research Archives Portal Repository