#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!# #-- assumes ss0 has been sourced --# #!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!# # -- Read Data & Make Measurement Matrix y=ts(scan("/mydata/jj.dat"),start=1960,freq=4) num=length(y) A=cbind(1,1,0,0) # -- Function to Calculate Likelihood -- Linn=function(para){ phi=para[1] Phi=diag(0,4); # Phi is 4x4 but only one element is a parameter Phi[1,1]=phi; Phi[2,]=c(0,-1,-1,-1) Phi[3,] = c(0,1,0,0); Phi[4,] = c(0,0,1,0) cQ1=para[2] # sqrt q11 cQ2=para[3] # sqrt q22 cQ=diag(0,4); cQ[1,1]=cQ1; cQ[2,2]=cQ2 cR=para[4] # sqrt r11 kf=Kfilter0(num,y,A,mu0,Sigma0,Phi,cQ,cR) return(kf$like) } # -- Initial Parameters -- mu0=c(.7,0,0,0) Sigma0= diag(.04,4) initpar=c(1.03,.1,.1,.5) # initial parameters for Phi[1,1], the 2 Q's and R # -- Estimation -- # the iteration number is printed to the screen but you have to manually scroll 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("Phi11","sigw1","sigw2","sigv") u #-- Smooth-- phi=est$par[1] Phi=diag(0,4); Phi[1,1]=phi; Phi[2,]=c(0,-1,-1,-1) Phi[3,] = c(0,1,0,0); Phi[4,] = c(0,0,1,0) cQ1=est$par[2] cQ2=est$par[3] cQ=diag(1,4); cQ[1,1]=cQ1; cQ[2,2]=cQ2 # note lower diag is 2x2 ident for inversions (as a device, they're not used) cR=est$par[4] ks=Ksmooth0(num,y,A,mu0,Sigma0,Phi,cQ,cR) #-- Plot -- Tsm=ts(ks$xs[1,,],start=1960,freq=4) Ssm=ts(ks$xs[2,,],start=1960,freq=4) p1=2*sqrt(ks$Ps[1,1,]) p2=2*sqrt(ks$Ps[2,2,]) par(mfrow=c(3,1)) plot(Tsm, main="Trend Component", ylab="Trend") lines(Tsm+p1,lty="dashed", col="blue") lines(Tsm-p1,lty="dashed", col="blue") plot(Ssm, main="Seasonal Component", ylim=c(-5,4), ylab="Season" ) lines(Ssm+p2,lty="dashed", col="blue") lines(Ssm-p2,lty="dashed", col="blue") plot(y, type="p", main="Data (points) and Trend+Season (line)") lines(Tsm+Ssm)