# ::Free Statistics and Forecasting Software::

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

### :: Example 6.6 - 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.6 (Newton-Raphson for Ex. 6.3) 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 ex66.txt). In addition, some hard-coded numbers have been replaced by user-defined parameters.
Note: many combinations of parameter values will not yield satisfactory results.

 Send output to: Browser Blue - Charts White Browser Black/White CSV MS Excel MS Word Number of values (?) Seed (?) phi (?) SD of state equation (?) SD of obs. equation (?)

 Source code of R module par1 <- as.numeric(par1) par2 <- as.numeric(par2) par3 <- as.numeric(par3) par4 <- as.numeric(par4) par5 <- as.numeric(par5) 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) } 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 # this makes sure Sigma0 is never negative kf = Kfilter(num,y,A,mu0,Sigma0,phi,sigw,sigv) return(kf\$like) } set.seed(par2) num=par1 #100 N=num+1 x =arima.sim(n=N, list(ar = par3, sd=par4)) # x[1]=x_0, x[2]=x_1,... v = rnorm(num,0,par5) y=x[2:N]+v data=ts(y) u=ts.intersect(data,lag(data,-1),lag(data,-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); warning = "" if (q < 0) { warning = "Warning: initial estimate for SD(w) could not be computed. Initial value is set to 1.
" q=1 } if (r < 0) { warning = paste(warning,"Warning: initial estimate for SD(v) could not be computed. Initial value is set to 1.
 Top | Output | Charts | References | History | Feedback

 Cite this software as: Wessa P., (2012), Reproduction of Example 6.6 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_shumwaystofferex66.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 Belisle, C. J. P. (1992) Convergence theorems for a class of simulated annealing algorithms on Rd., J Applied Probability, vol. 29, 885-895. Byrd, R. H., Lu, P., Nocedal, J. and Zhu, C. (1995) A limited memory algorithm for bound constrained optimization., SIAM J. Scientific Computing, vol. 16, 1190-1208. Fletcher, R. and Reeves, C. M. (1964) Function minimization by conjugate gradients., Computer Journal, vol. 7, 148-154. Nash, J. C. (1990), Compact Numerical Methods for Computers., Linear Algebra and Function Minimisation., Adam Hilger. Nelder, J. A. and Mead, R. (1965), A simplex algorithm for function minimization., Computer Journal, vol. 7, 308-313. Nocedal, J. and Wright, S. J. (1999), Numerical Optimization., Springer.
 Top | Output | Charts | References | History | Feedback
 Top | Output | Charts | References | History | Feedback