Bayesian Static Parameter Estimation for Partially Observed Diffusions via Multilevel Monte Carlo 7 1 0 Ajay Jasra∗1, Kengo Kamatani†2, Kody J. H. Law‡3, and Yan Zhou§4 2 n a 1,4Department of Statistics Applied Probability, National University of Singapore, J 0 SG 2 2Department of Engineering Science, Osaka University, JP ] O 3Computer Science and Mathematics Division, Oak Ridge National Laboratory, C . t USA a t s [ 1 v 2 Abstract 9 8 InthisarticleweconsiderstaticBayesianparameterestimationforpartiallyobserveddif- 5 0 fusionsthatarediscretelyobserved. Weworkundertheassumptionthatonemustresortto . 1 discretizingtheunderlyingdiffusionprocess,forinstanceusingtheEulerMaruyamamethod. 0 7 Given this assumption, we show how one can use Markov chain Monte Carlo (MCMC) and 1 : particularly particle MCMC [Andrieu, C., Doucet, A. & Holenstein, R. (2010). Particle v i Markov chain Monte Carlo methods (with discussion). J. R. Statist. Soc. Ser. B, 72, 269– X r 342]toimplementanewapproximationofthemultilevel(ML)MonteCarlo(MC)collapsing a sum identity. Our approach comprises constructing an approximate coupling of the poste- rior density of the joint distribution over parameter and hidden variables at two different discretization levels and then correcting by an importance sampling method. The variance of the weights are independent of the length of the observed data set. The utility of such a method is that, for a prescribed level of mean square error, the cost of this MLMC method ∗[email protected] †[email protected] ‡[email protected] §[email protected] 1 is provably less than i.i.d. sampling from the posterior associated to the most precise dis- cretization. However the method here comprises using only known and efficient simulation methodologies. The theoretical results are illustrated by inference of the parameters of two prototypicalprocessesgivennoisypartialobservationsoftheprocess: thefirstisanOrnstein Uhlenbeck process and the second is a more general Langevin equation. Key words: Multilevel Monte Carlo, Markov chain Monte Carlo, Diffusion Processes 1 Introduction The Hidden Markov Model (HMM) is widely used in many disciplines, including applied mathe- matic, statistics, economics and finance; see [2] for an overview. In this article, we are interested in HMMs given by diffusions which are partially observed, discretely in time. In particular, we assumethatinordertofitthemodeltothedata,onemust resorttoadiscretizationofthediffu- sion, for instance, using Euler-Maruyama. In addition, we assume that associated to the model isastatic(non-time-varying)finitedimensionalparameter, whichoneisinterestedtoinfergiven afixedlengthdatarecord. Insimpleterms, thediscretization, oflevelhsay, whereash→0one obtainstheexactdiffusion,inducesaposteriorsayπ onthestaticparameterθandhiddenstates h at the observation times, say X . We seek to approximate E [ϕ(θ,X )] for appropriately 0:n πh 0:n defined real-valued functions. Ultimately, one might seek to remove the dependence upon h and gettheexactexpectationwithnodiscretizationbias. Weremarkthatthemodelwillbeformally introduced in the next section. This framework is relevant to a broad range of applications in science and engineering; see [2, 17] The task of computing the expectation for any fixed h>0 is a non-trivial task, which often requires quite advanced Monte Carlo methods. As has been remarked in many articles in the literature, ofen the joint correlation between θ and X means even standard MCMC methods 0:n may produce very inaccurate of inefficient approximations of the expectation of interest, despite theirtheoreticalvalidity. Animportantalgorithmthathas,toanextent,helpedtoalleviatethese difficulties is the particle MCMC (PMCMC) methods of [1] and their subsequent developments (e.g. [4]). Intrinsically, this method uses a sequential Monte Carlo (SMC) (e.g. [7]) method to help move the samples around the state-space, for instance, inside a Metropolis-Hastings acceptance/rejection scheme, although Gibbs versions also exist. PMCMC delivers a Markov 2 chain which provides consistent estimates of expectations of the form E [ϕ(θ,X )], for any πh 0:n fixed h SMC methods are well-known as being efficient techniques for filtering, when the state- variable at time k, X , is of moderate to low dimension and all the static parameters are fixed. k In the context of this article, there is an additional degree of freedom, which can be utilized to further enhance the PMCMC method. This is associated to the discretization level h. We consider using the multilevel Monte Carlo (MLMC) framework [8, 9, 11]. This allows one to leverage in an optimal way the nested problems arising in this context, hence minimizing the necessarycosttoobtainagivenlevelofmeansquareerror. Setπ astheposterioronθ,X with 0:n nodiscretizationbiasandπ asthetime-discretizedposterioronθ,X withtimediscretization hl 0:n h , one has for an intergrable, real-valued function ϕ and +∞ > h > h > ··· > h > 0 (the l 0 1 L levels) L (cid:88) E [ϕ(θ,X )]= {E [ϕ(θ,X )]−E [ϕ(θ,X )]} (1) πhL 0:n πhl 0:n πhl−1 0:n l=0 where E is the expectation operator and E [ϕ(θ,X )] := 0. The idea of MLMC is then πh−1 0:n to approximate each summand by independently simulating N samples from a dependent cou- l pling of (π ,π ). In such scenarios, one can show that the overall mean square error (MSE) hl hl−1 associated to the approximation of E [ϕ(θ,X )] is: π 0:n L MSE=Bias(L,ϕ)2+(cid:88) Vl , (2) N l l=0 where Bias(L,ϕ)=|E [ϕ(θ,X )]−E [ϕ(θ,X )]| , (3) πhL 0:n π 0:n and 0 < V < +∞ are a collection of constants. It is remarked that it is the coupled samples l which induce V to be a function of h which is often critical as we explain below. Assuming the l l cost of C per level, per sample, the cost of the algorithm is then (cid:80)L C N . Fixing (cid:15) > 0 and l l=0 l l given an appropriate parameterization of h (e.g. h = 2−l), one then chooses L to ensure that l l Bias(L,ϕ)2 =O((cid:15)2) and then given C ,V characterised as a function of h optimizes N ,...,N l l l 0 L tominimizethecostsothattheterm(cid:80)L Vl =O((cid:15)2); [8]givesthesolutiontothisconstrained l=0 Nl optimization problem. In many scenarios of practical interest the associated MLMC algorithm can achieve a MSE of O((cid:15)2) at a cost which is less than i.i.d. sampling from π ; note that this hL has not yet been established in the problem under study here. The main issue is that sampling independently from the couples (π ,π ) is not possible in our context. hl hl−1 Inthispaperweshowhowtoimplementanewapproximationofthemultilevelcollapsingsum 3 identity. Our approach comprises constructing an approximate coupling of the posterior density of the joint on the parameter and hidden space at two different discretization levels and then correcting by an importance sampling method, whose variance of the weights are independent of the length of the observed data set. The utility of such a method is that it comprises using known and efficient simulation methodologies, instead of coupling algorithms as explored in [13, 14, 15, 19]. In particular, our approach facilitates a mathematical analysis which allows us toestablishthatourapproachcanbebetterthansampling(e.g.byPMCMC)fromtheposterior associated to the most precise discretization. The algorithm presented here is distinct from either of the previously introduced multilevel MCMC (MLMCMC) algorithms [12, 16], and may be generalized. This article is structured as follows. In Section 2 the model is described. In Section 3 we describe our approach and give a mathematical result associated to the MSE of the method. In Section 4 we give practical simulations to establish the theory. The appendix contains some of the proofs for the result of Section 3. 2 Model We consider the following partially-observed diffusion process: dX = a (X )dt+b (X )dW (4) t θ t θ t t with X ∈ Rd = X, t ≥ 0, X has initial probability density f and {W } a Brownian t 0 θ t t∈[0,T] motion of appropriate dimension. θ ∈ Θ ⊆ Rdθ is a static parameter of interest. The following assumptions will be made on the diffusion process. Assumption 2.1. a :Rd →Rd, b :Rd →Rd×d satisfy θ θ (i) global Lipschitz property: there is a C >0 such that |a (x)−a (y)|+|b (x)−b (y)|≤ θ θ θ θ C|x−y| for all x,y ∈X and all θ ∈Θ; (ii) bounded moments: sup E |X |p <∞ for all p≥1. θ∈Θ θ 0 Notice that (i) and (ii) together imply that E |X |p <∞ for all n. θ n It will be assumed that the data are regularly spaced (i.e. in discrete time) observations y ,...,y , y ∈ Rm = Y. It is assumed that conditional on X , for discretization δ > 0, Y is 1 n k kδ k 4 independent of all other random variables with density g (x ,y ). For simplicity of notation let θ kδ k δ = 1 (which can always be done by rescaling time), so X = X . It is noted that we assume k kδ that one does not have access to a non-negative and unbiased estimate of the transition density of the diffusion and we are forced to work with a discretized process. The above formulation can then summarized as follows, on discretizing the diffusion process with discretization level h. We have a pair of discrete-time stochastic processes, {X } and n n≥0 {Y } , where X ∈ X (with associated σ−algebra X) is an unobserved process and y ∈ Y n n≥1 n n (withassociatedσ−algebraY)isobserved. Letθ ∈Θ⊆Rdθ beaparameter. Thehiddenprocess {X } is a Markov chain with initial density f at time 0 and transition density f (x ,x ), n θ θ,h p−1 p i.e. for each θ ∈Θ (cid:90) (cid:90) P (X ∈A)= f (x)dx and P (X ∈A|X =x )= f (x ,x )dx p≥1 θ,h 0 θ θ,h p p−1 p−1 θ,h p−1 p p A A (5) whereP denotesprobability,A∈X anddx isadominatingσ-finitemeasure. Inaddition,the θ,h n observations{Y } conditionedupon{X } arestatisticallyindependentandhavemarginal n n≥1 n n≥0 density g (x ,y ), i.e. θ n n (cid:90) P (Y ∈B|{X } ={x } )= g (x ,y )dy n≥1 (6) θ,h n k k≥0 k k≥0 θ n n n B with B ∈ Y and dy the dominating σ-finite measure. The HMM is given by equations (5)-(6) n and is often referred to in the literature as a state-space model. In our context θ ∈ Θ is a parameter of interest with prior π . θ Given the joint density on U:=Θ×Xn+1 n (cid:89) π (θ,x )∝π (θ)f (x ) g (x ,y )f (x ,x ) , h 0:n θ θ 0 θ p p θ,h p−1 p p=1 for ϕ ∈ B (U)∩Lip(U), where B (U) are the bounded and real-valued measurable functions on b b U and Lip(U) arethe Lipschitz, measurable functions on U, and for+∞>h >···>h >0 we 0 L would like to compute L (cid:88)(cid:110) (cid:111) E [ϕ(θ,X )]= E [ϕ(θ,X )]−E [ϕ(θ,X )] (7) πhL 0:n πhl 0:n πhl−1 0:n l=0 where E [·]=0. We will use the MLMC approach. πh−1 Consider only a single pair E [ϕ(θ,X )]−E [ϕ(θ,X )], h<h(cid:48). It is well known that if πh 0:n πh(cid:48) 0:n onecansamplefromadependentcouplingof(π ,π ),suchasthemaximalcoupling,thenMonte h h(cid:48) 5 Carlo estimation of such a difference can be performed at a lower cost than i.i.d sampling from the independent coupling of (π ,π ) [8, 9]. The main issue is that such couplings are typically h h(cid:48) not available up-to a non-negative and unbiased estimator. We consider the scenario where one samples from a sensible, approximate, coupling and corrects via importance sampling. 3 Method and Analysis 3.1 Method We are to approximate the identity (7). Our procedure, when considering the summands from 1,...,L will be to run L independent pairs of the idea to be described below. The case l = 0 is simply using (e.g.) PMCMC to approximate E [ϕ(θ,X )]; we refer the reader to [1] for πh0 0:n details on PMCMC - a simple decsription is below. We only consider a pair E [ϕ(θ,X )]− πh 0:n E [ϕ(θ,X )], h<h(cid:48). The methodology and analysis in this context of one pair will suffice to πh(cid:48) 0:n justify our approach as we will explain below. Letz =(x,x(cid:48))∈X×X=ZandQ (z,z¯)beanycoupling(otherthantheindependentone) θ,h,h(cid:48) of (f (x,x¯),f (x(cid:48),x¯(cid:48))). For instance, in the context of an Euler discretization a description θ,h θ,h(cid:48) can be found in [15] (see also appendix B). Let G (z) = max{g (x,y ),g (x(cid:48),y )} (note that p,θ θ p θ p alternative choices of G are possible). We propose to sample from the probability density on p,θ V=Θ×X2n+2 (write the associated σ−algebra as V) n (cid:89) π (θ,z )∝π (θ)ν (z ) G (z )Q (z ,z ). h,h(cid:48) 0:n θ θ 0 p,θ p θ,h,h(cid:48) p−1 p p=1 Then for ϕ∈B (U)∩Lip(U): b E [ϕ(θ,X )]−E [ϕ(θ,X )]= πh 0:n πh(cid:48) 0:n E [ϕ(θ,X )H (θ,Z )] E [ϕ(θ,X(cid:48) )H (θ,Z )] πh,h(cid:48) 0:n 1,θ 0:n − πh,h(cid:48) 0:n 2,θ 0:n (8) E [H (θ,Z )] E [H (θ,Z )] πh,h(cid:48) 1,θ 0:n πh,h(cid:48) 2,θ 0:n where n H (θ,z ) = (cid:89) gθ(xp,yp) 1,θ 0:n G (z ) p,θ p p=1 H (θ,z ) = (cid:89)n gθ(x(cid:48)p,yp). 2,θ 0:n G (z ) p,θ p p=1 We note that our choice of G (z) ensures that H and H are uniformly upper-bounded by p,θ 1,θ 2,θ 1 and hence that the variance w.r.t. any probability is independent of n. 6 3.1.1 Particle MCMC Let (W,W) be a measurable space such that V ⊆ W. Let K : W × W → [0,1] be any er- godic Markov kernel ofinvariant measure η such that onecanconsistentlyestimate expectations w.r.t. π . For instance, if for every A∈V h,h(cid:48) (cid:90) (cid:90) η(dw)= π (θ,z )d(θ,z ). h,h(cid:48) 0:n 0:n A×(W\V) A Our construction allows a particle MCMC approach to be adopted, which is not quite as the displayed equation, but nonetheless allows one to infer π . We focus on one particle MCMC h,h(cid:48) method for completeness, but, we reiterate that one can use the analysis here for more advanced versions of the algorithm, or indeed, any MCMC of the form above. We will now describe the particle marginal Metropolis-Hastings (PMMH) algorithm. Let M ≥1 and θ be fixed, and introduce random variables a ∈{1,...,M}n, which will denote 0:n−1 the indices of the selected particles upon resampling at the given steps. One can run a particle filter [5] to approximate n (cid:89) π (z |θ)∝ν (z ) G (z )Q (z ,z ) h,h(cid:48) 0:n θ 0 p,θ p θ,h,h(cid:48) p−1 p p=1 by sampling from the following joint, on the space {1,...,M}n×ZM(n+1) p(a1:M ,z1:M|θ)=(cid:16)(cid:89)M ν (zi)(cid:17)(cid:89)n (cid:89)M (cid:16) Gp−1,θ(zpa−ip−11) Q (zaip−1,zi)(cid:17) , (9) 0:n−1 0:n θ 0 (cid:80)M G (zj ) θ,h,h(cid:48) p−1 p i=1 p=1i=1 j=1 p−1,θ p−1 where G := 1. Note that better algorithms can be constructed, but we just present the most 0,θ simple approach. We remark that n M (cid:89)(cid:16) 1 (cid:88) (cid:17) pM (y |θ)= G (zj) (10) h,h(cid:48) 0:n M p,θ p p=1 j=1 is an unbiased estimator of p (y |θ) = (cid:82) ν (z )(cid:81)n G (z )Q (z ,z )dz ; see h,h(cid:48) 0:n Zn+1 θ 0 p=1 p,θ p θ,h,h(cid:48) p−1 p 0:n [5]. The PMMH algorithm works as follows. The superscripts for (θ,k) are the iteration (time) counter of the MCMC. 1. Initialize: Sampleθ0 fromthepriorandthensample(a1:M ,z1:M)fromp(a1:M ,z1:M|θ0) 0:n−1 0:n 0:n−1 0:n as in (9), and store pM (y |θ0) as in (10). Select a path zj , constructed by drawing zj h,h(cid:48) 0:n 0:n n with probability proportional to G (zj), and setting (zj(cid:48) |zj(cid:48)) = zajp(cid:48)−1; set k0 as the n,θ0 n p−1 p p−1 index of the selected path. Set i=1. 7 2. Iterate: Sampleθ(cid:48)|θi−1accordingtoaproposalwithconditionaldensityq(θ(cid:48)|θi−1)thenfrom p(a1:M ,z1:M|θ(cid:48)) as in (9). Select a path zj with probability proportional to G (zj) 0:n−1 0:n 0:n n,θ(cid:48) n and constructed as described above; set k(cid:48) as the index of the selected path. Set θi = θ(cid:48), ki =k(cid:48) with probability: 1∧ pMh,h(cid:48)(y0:n|θ(cid:48)) πθ(θ(cid:48))q(θi−1|θ(cid:48)) pM (y |θi−1)π (θi−1)q(θ(cid:48)|θi−1) h,h(cid:48) 0:n θ otherwise θi =θi−1, ki =ki−1. Set i=i+1 and return to the start of 2. We denote by K the PMMH kernel and denote by (W,W) the measurable space for which it is defined upon. The invariant measure is denoted η. For the analysis, we assume the MCMC algorithm is started in stationarity. Then one estimates (8) by 1 (cid:80)N ϕ(θi,xki )H (θi,zki ) 1 (cid:80)N ϕ(θi,x(cid:48)ki)H (θi,zki ) N i=1 0:n 1,θi 0:n − N i=1 0:n 2,θi 0:n . 1 (cid:80)N H (θi,zki ) 1 (cid:80)N H (θi,zki ) N i=1 1,θi 0:n N i=1 2,θi 0:n This estimate is consistent in the limit as N grows; see [1]. To simplify the notation we replace ki in the superscripts by i from here on. 3.2 Multilevel Considerations As described for MLMC in the introduction, we will approximate the expectation using the telescopic sum identity given in (1). We will establish error estimates for L (cid:88)E¯Nl(ϕ), E¯Nl(ϕ)=ENl(ϕ)−E (ϕ) , (11) l l l l l=0 where 1 (cid:80)Nl ϕ(θi,xi )H (θi,zi ) 1 (cid:80)Nl ϕ(θi,x(cid:48)i )H (θi,zi ) ENl(ϕ)= Nl i=1 0:n 1,θi 0:n − Nl i=1 0:n 2,θi 0:n (12) l 1 (cid:80)Nl H (θi,zi ) 1 (cid:80)Nl H (θi,zi ) Nl i=1 1,θi 0:n Nl i=1 2,θi 0:n is a consistent estimator of E (ϕ) := E [ϕ(θ,X )]−E [ϕ(θ,X )]. Therefore (11) is a l πhl 0:n πhl−1 0:n consistent estimator of E [ϕ(θ,X )] and the the MSE (2) can be bounded, up to a constant, πhL 0:n by the sum of the squared error of (11) and Bias(L,ϕ)2, as given by (3), which is O(h ) for L example using Euler Maruyama. Using E to denote the expectation w.r.t. the law associated to our algorithm, assuming the Markov chain is started in stationarity, our objective is therefore to investigate L L E[((cid:88)E¯Nl(ϕ))2]=(cid:88)E[E¯Nl(ϕ)2] (13) l l l=0 l=0 8 soastooptimallyallocateN ,...,N asdescribedintheintroduction. Thuswemustinvestigate 0 L terms such as E[E¯Nl(ϕ)2] for a given l. l 3.3 Analysis Below P(W) are the collection of probability measures on (W,W). (A1) For every y ∈Y there exist 0<C <C <+∞ such that for every x∈X, θ ∈Θ, C ≤g (x,y)≤C. θ For every y ∈Y, g (x,y) is globally Lipschitz on X×Θ. θ (A2) For any 0 ≤ k ≤ n, q ∈ {1,2} there exists a β > 0 such that for any ϕ ∈ B (Θ×Xk+1)∩ b Lip(Θ×Xk+1) there exists a C <+∞ (cid:32)(cid:90) k (cid:33)3−q (cid:89) |ϕ(θ,x )−ϕ(θ,x(cid:48) )|q Q (z ,z )π (θ)ν (z )dθdz ≤C(h(cid:48))β. 0:k 0:k θ,h,h(cid:48) k−1 k θ θ 0 0:k Θ×X2k+2 p=1 (A3) Supposethatforanyn>0thereexistaξ ∈(0,1)andν ∈P(W)suchthatforeachw ∈W, ϕ∈B (W)∩Lip(W), h,h(cid:48): b (cid:90) (cid:90) ϕ(w(cid:48))K(w,dw(cid:48))≥ξ ϕ(w)ν(dw). W W (cid:82) (cid:82) K is η-reversible, that is, η(dw)K(w,A)= η(dw)K(w,B) for any A,B ∈W. w∈B w∈A We note that (A1) can be verified for some state-space models (especially if Y and Θ are compact) and (A3) can be verified for a PMCMC kernel, if Θ,X are compact - indeed, the constants would all be independent of n under appropriate settings of the algorithm. Theorem 3.1. Assume (A1-3). Then for any n > 0, there exists a β > 0 such that for any ϕ∈B (Θ×Xn+1)∩Lip(Θ×Xn+1) there exists a C <+∞ such that b (cid:34)(cid:32) 1 (cid:80)N ϕ(θi,xi )H (θi,zi ) 1 (cid:80)N ϕ(θi,x(cid:48)i )H (θi,zi ) E N i=1 0:n 1,θi 0:n − N i=1 0:n 2,θi 0:n 1 (cid:80)N H (θi,zi ) 1 (cid:80)N H (θi,zi ) N i=1 1,θi 0:n N i=1 2,θi 0:n −(cid:32)Eπh,h(cid:48)[ϕ(θ,X0:n)H1,θ(θ,Z0:n)] − Eπh,h(cid:48)[ϕ(θ,X0(cid:48):n)H2,θ(θ,Z0:n)](cid:33)(cid:33)2(cid:35)≤ C(h(cid:48))β. E [H (θ,Z )] E [H (θ,Z )] N πh,h(cid:48) 1,θ 0:n πh,h(cid:48) 2,θ 0:n Proof. The result follows by using Lemma C.3. of [14], the C −inequality, the boundedness of 2 certain quantities and Proposition A.1.The proof is omitted as it is similar to the calculations in [14]. 9 3.4 A Return to Multilevel Considerations Returning to Section 3.2, we assume that h =2−l and introduce the further assumption l Assumption 3.1. The cost to simulate ENl in (12) is controlled by C(ENl)≤CN h−γ, and the l l l l bias is controlled by |E (ϕ(θ,X ))−E (ϕ(θ,X ))|≤Chα , πhL 0:n π 0:n L for γ,α,C >0. Following assumption (A2), α=β/2 satisfies the above, but it may be larger, e.g. for Euler- Maruyama in which α=β. Given (cid:15) > 0, in order to ensure the MSE is O((cid:15)2), the term (3) must be O((cid:15)2). Following from Assumption (A2), it suffices to let L∝2|log((cid:15))|/β so that h =(cid:15). L Following from Theorem 3.1, (cid:88)L E[E¯Nl(ϕ)2]≤C(cid:88)L hβl , l N l l=0 l=0 andnotethattheconstantC maydependuponthetimeparametern,whichhasbeensuppressed from the notation; we return to this point below. Suppose we minimize COST = (cid:80)Ll=0h−l γNl subject to (cid:80)Ll=0 hNβll = O((cid:15)2) as a function of N ,...,N . This is exactly considered in [8] for γ =1 and later in [3] for γ (cid:54)=1, and yields that 0 L N ∝ε−2K h(β+γ)/2, (14) l L l whereK =(cid:80)L h(β−γ)/2 (seealso[14,6]). ThisgivesacostofO(ε−2K2)pertimestep. Hence L l=1 l L the following corollary is immediate. Corollary 3.1 (ML Cost). Given (A1-3) and Assumption 3.1, for any n > 0 and any ϕ ∈ B (Θ×Xn+1)∩Lip(Θ×Xn+1), (L,{N }L ) can be chosen such that the estimator (cid:80)L ENl(ϕ), b l l=1 l=1 l with ENl given in (12), satisfies l (cid:34) L (cid:35) (cid:88) E | ENl(ϕ)−E (ϕ(θ,X ))|2 ≤C(cid:15)2 , l π 0:n l=1 for some C >0, for a total cost controlled by  (cid:15)−2, if β >γ, COST≤C (cid:15)−2|log((cid:15))|2, if β =γ, (15) (cid:15)−(2+γ−αβ), if β <γ. 10