#---- this is example 6.10 using ss1 ---# ######################################### # detailed comments are in that file # # if you need them # #!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!# #-- assumes ss1 has been sourced --# #!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!# # -- Read Data & Make A Matrix y=ts(scan("/mydata/jj.dat"),start=1960,freq=4) num=length(y) A=array(cbind(1,1,0,0), dim=c(1,4,num)) # this and a call to Kfilter1 and Ksmooth1 below # are the only differences between the 2 examples # -- Function to Calculate Likelihood -- Linn=function(para){ phi=para[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=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=Kfilter1(num,y,A,mu0,Sigma0,Phi,0,0,cQ,cR,0) return(kf\$like) } # -- Initial Parameters -- mu0=c(.7,0,0,0) Sigma0= diag(.04,4) initpar=c(1.03,.1,.1,.5) # initial parameters # -- Estimation -- 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 cR=est\$par[4] ks=Ksmooth1(num,y,A,mu0,Sigma0,Phi,0,0,cQ,cR,0) #-- 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)