#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!# #-- assumes ss0 has been sourced --# #!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!# # -- Read Data -- y1=scan("/mydata/HL.dat") y2=scan("/mydata/Folland.dat") # -- Function to Calculate Likelihood -- Linn=function(para){ cQ=para[1] # sigma_w cR1=para[2] # 11 element of chol(R) cR2=para[3] # 22 element of chol(R) cR12=para[4] # 12 element of chol(R) cR=matrix(c(cR1,0,cR12,cR2),2) # put the matrix together kf=Kfilter0(num,y,A,mu0,Sigma0,Phi,cQ,cR) return(kf$like) } # -- Setup -- y=cbind(y1,y2) # y is data matrix (108 x 2) num=nrow(y) A=matrix(1,2,1) # A is a 2x1 matrix of 1s - same for all t mu0=-.35 Sigma0= .01 Phi=1 initpar=c(.1,.1,.1,0) # initial parameter values in order, para[1] to para[4] # -- 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))) # #-- display summary of estimation estimate=est$par u=cbind(estimate,stderr) rownames(u)=c("sigw","cR11", "cR22", "cR12") u #-- Smooth-- #-- first set parameters to their final estimates cQ=est$par[1] cR1=est$par[2] # 11 element of chol(R) cR2=est$par[3] # 22 element of chol(R) cR12=est$par[4] # 12 element of chol(R) cR=matrix(c(cR1,0,cR12,cR2),2) ks=Ksmooth0(num,y,A,mu0,Sigma0,Phi,cQ,cR) #-- Plot -- plot.ts(as.vector(ks$xs),lwd=2,ylim=c(-.7,.4), ylab="Temp Deviations") lines(y1,col="blue",lty="dashed") # color helps here lines(y2,col="red", lty="dashed")