################################################### # --- MUST source ss0 and EM0 and install nlme ---# ################################################### ######################################### library(nlme) # loads package nlme ######################################### #-- generate data (same as ex66, see that file for more details) set.seed(999) num=100 N=num+1 # need 101 x's # generate x(t)=.8x(t-1)+w(t), t=1,...,100 # and x(0) from the stationary distribution: x = arima.sim(n=N, list(ar = .8, sd=1)) # here you have x0 to x100 v = rnorm(num,0,1) # obs noise y=ts(x[-1]+v) # observations y[1] ,..., y[100] #-- initial estimates u=ts.intersect(y,lag(y,-1),lag(y,-2)); varu=var(u); coru=cor(u); phi=coru[1,3]/coru[1,2]; q = (1-phi^2)*varu[1,2]/phi; r = varu[1,1] - q/(1-phi^2); cr=sqrt(r) cq=sqrt(q) mu0=0 Sigma0=2.8 #-- estimation (stop at 75 iterations or when -loglike changes by < .001) #--- EM0 will display the iteration number, but you have to scroll manually to see it (em=EM0(num,y,1,mu0,Sigma0,phi,cq,cr,75,.001)) # -- Standard Errors -- uses package nlme phi=em$Phi cq=chol(em$Q) cr=chol(em$R) mu0=em$mu0 Sigma0=em$Sigma0 para=c(phi,cq,cr) #-- evaluates likelihood at estimates Linn=function(para){ kf = Kfilter0(num,y,1,mu0,Sigma0,para[1],para[2],para[3]) return(kf$like) } emhess=fdHess(para, function(para) Linn(para)) stderr=sqrt(diag(solve(emhess$Hessian))) #-- display summary of estimation estimate=c(para, em$mu0, em$Sigma0) se=c(stderr,NA,NA) u=cbind(estimate,se) rownames(u)=c("phi","sigw","sigv","mu0","Sigma0") u