# ::Free Statistics and Forecasting Software::

v1.1.23-r7
Secure website (SSL) | Donate to wessa.net

### :: Example 6.5 - Time Series Analysis and Its Applications (Edition 2: With R Examples) ::

All rights reserved. The non-commercial (academic) use of this software is free of charge. The only thing that is asked in return is to cite this software when results are used in publications.

This free online software (calculator) reproduces Example 6.5 (Prediction, Filtering, and Smoothing for the Local Level Model) in Shumway R.H., Stoffer D.S., Time Series Analysis and Its Applications (2nd ed. with R examples). The module includes all necessary functions to reproduce their results (Kfilter.R, Ksmooth.R, and ex65.txt). In addition, some hard-coded numbers have been replaced by user-defined parameters.

 Send output to: Browser Blue - Charts White Browser Black/White CSV MS Excel MS Word Number of values (?) Seed (?) Chart options Width: Height: Y-axis minimum Y-axis maximum

 Source code of R module par1 <- as.numeric(par1) par2 <- as.numeric(par2) Kfilter = function(num,y,A,mu0,Sigma0,Phi,cQ,cR){ Q=t(cQ)%*%cQ R=t(cR)%*%cR N=num+1 nseries=ncol(as.matrix(y)) yobs=matrix(0,N,nseries) yobs[2:N,]=y # yobs is y with first row=zeros xp=vector("list",N) # xp=x_t^{t-1} Pp=vector("list",N) # Pp=P_t^{t-1} xf=vector("list",N) # xf=x_t^t Pf=vector("list",N) # Pf=x_t^t innov=vector("list",N) # innovations sig=vector("list",N) # innov var-cov matrix like=0 # -log(likelihood) xf[[1]]=mu0 Pf[[1]]=Sigma0 for (i in 2:N){ xp[[i]]=Phi%*%xf[[i-1]] Pp[[i]]=Phi%*%Pf[[i-1]]%*%t(Phi)+Q siginv=A[[i]]%*%Pp[[i]]%*%t(A[[i]])+R sig[[i]]=(t(siginv)+siginv)/2 # make sure sig is symmetric siginv=solve(sig[[i]]) # now siginv is sig[[i]]^{-1} K=Pp[[i]]%*%t(A[[i]])%*%siginv innov[[i]]=as.matrix(yobs[i,])-A[[i]]%*%xp[[i]] xf[[i]]=xp[[i]]+K%*%innov[[i]] Pf[[i]]=Pp[[i]]-K%*%A[[i]]%*%Pp[[i]] like= like + log(det(sig[[i]])) + t(innov[[i]])%*%siginv%*%innov[[i]] } like=0.5*like list(xp=xp,Pp=Pp,xf=xf,Pf=Pf,like=like,innov=innov,sig=sig,Kn=K) } Ksmooth = function(num,y,A,mu0,Sigma0,Phi,cQ,cR){ kf=Kfilter(num,y,A,mu0,Sigma0,Phi,cQ,cR) N=num+1 xs=vector("list",N) # xs=x_t^n Ps=vector("list",N) # Ps=P_t^n J=vector("list",N) # J=J_t xs[[N]]=kf\$xf[[N]] Ps[[N]]=kf\$Pf[[N]] for(k in N:2) { J[[k-1]]=(kf\$Pf[[k-1]]%*%t(Phi))%*%solve(kf\$Pp[[k]]) xs[[k-1]]=kf\$xf[[k-1]]+J[[k-1]]%*%(xs[[k]]-kf\$xp[[k]]) Ps[[k-1]]=kf\$Pf[[k-1]]+J[[k-1]]%*%(Ps[[k]]-kf\$Pp[[k]])%*%t(J[[k-1]]) } list(xs=xs,Ps=Ps,J=J,xp=kf\$xp,Pp=kf\$Pp,xf=kf\$xf,Pf=kf\$Pf,like=kf\$like,Kn=kf\$K) } set.seed(par2) num=par1 #50 w = rnorm(num+1,0,1) v = rnorm(num,0,1) mu = ts(cumsum(w)) # Note: states mu_0=mu[1], mu_1=mu[2],...,mu_n=mu[num+1] np1 <- par1+1 y=mu[2:np1]+v # because I cannot figure out how to use 0 for an index A=matrix(1,np1,1) A=as.list(A,np1) # A has to be a list - need this for mu0=0 sigma0=1 phi=1 cQ=1 # These are the Cholesky decomps of Q and R cR=1 # Why? For optimization later. kf=Kfilter(par1,y,A,mu0,sigma0,phi,cQ,cR) mup=unlist(kf\$xp) # .. p for predictor Pp=unlist(kf\$Pp) muf=unlist(kf\$xf) # .. f for filter Pf=unlist(kf\$Pf) ks=Ksmooth(par1,y,A,mu0,sigma0,phi,cQ,cR) # you can also get mup, Pp, muf, Pf from here... e.g., mup=unlist(ks\$xp) mus=unlist(ks\$xs) Ps=unlist(ks\$Ps) # .. s for smoother bitmap(file="test1.png") if ((ylimmin == "") | (ylimmax == "")) plot.ts(mu[-1],type="p", main="Prediction", ylim=NULL) else plot.ts(mu[-1],type="p", main="Prediction", ylim=c(ylimmin,ylimmax)) lines(mup) lines(mup+2*sqrt(Pp), lty="dashed", col="blue") lines(mup-2*sqrt(Pp), lty="dashed", col="blue") dev.off() bitmap(file="test2.png") if ((ylimmin == "") | (ylimmax == "")) plot.ts(mu,type="p", main="Filter", ylim=NULL) else plot.ts(mu,type="p", main="Filter", ylim=c(ylimmin,ylimmax)) lines(muf) lines(muf+2*sqrt(Pf), lty="dashed", col="blue") lines(muf-2*sqrt(Pf), lty="dashed", col="blue") dev.off() bitmap(file="test3.png") if ((ylimmin == "") | (ylimmax == "")) plot.ts(mu,type="p", main="Smoother", ylim=NULL) else plot.ts(mu,type="p", main="Smoother", ylim=c(ylimmin,ylimmax)) lines(mus) lines(mus+2*sqrt(Ps), lty="dashed", col="blue") lines(mus-2*sqrt(Ps), lty="dashed", col="blue") dev.off()
 Top | Output | Charts | References | History | Feedback

 Cite this software as: Wessa P., (2008), Reproduction of Example 6.5 in Time Series Analysis and Its Applications (v1.0.5) in Free Statistics Software (v1.1.23-r7), Office for Research Development and Education, URL http://www.wessa.net/rwasp_shumwaystofferex65.wasp/ The R code is based on : Shumway R.H., Stoffer D.S., Time Series Analysis and Its Applications (2nd ed. with R examples), Springer, 2006 Shumway R.H., Stoffer D.S., Time Series Analysis and Its Applications (2nd ed. with R examples), URL http://www.stat.pitt.edu/stoffer/tsa2/#Rchapter6
 Top | Output | Charts | References | History | Feedback
 Top | Output | Charts | References | History | Feedback