#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!# #-- assumes ss0 has been sourced --# #!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!# #-- generate data 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 (method discussed in the text) 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); initpar=c(phi,sqrt(q),sqrt(r)) initpar # view the initial estimates #-- function to evaluate innovations likelihood Linn=function(para){ phi=para[1] sigw=para[2] # this is the standard dev sigv=para[3] # this is the standard dev mu0=0 Sigma0=(sigw^2)/(1-phi^2) Sigma0[Sigma0<0] =0 # make sure Sigma0 is never negative kf = Kfilter0(num,y,1,mu0,Sigma0,phi,sigw,sigv) # run filter under present parameters return(kf$like) # return -log likelihood } #-- do the estimation #-- view help(optim) for details here... # the iteration number is written to the screen but you have to manually scroll down to see it est=optim(initpar, Linn, NULL, method = "BFGS", hessian = TRUE, control=list(trace=1,REPORT=1)) stderr=sqrt(diag(solve(est$hessian))) # standard errors #-- display summary of estimation estimate=est$par u=cbind(estimate,stderr) rownames(u)=c("phi","sigw","sigv") u