#-- this is example 6.5 using ss1 --# #!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!# #-- assumes ss1 has been sourced --# #!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!# #~~~~~ detailed comments are in ex65.txt ~~~~~# # the only difference here is that # A is an array of num 1x1 matrices (see line 19) #-- generate data -- set.seed(1) num=50 w = rnorm(num+1,0,1) v = rnorm(num,0,1) mu = array(cumsum(w)) # states: mu[0], mu[1],...,mu[50] is a random walk, no drift y=mu[-1]+v # obs: y[1],..., y[50] mu[-1] is mu vector w/o mu[0] A=array(1,dim=c(1,1,num)) # array of num 1x1 matrices of 1s # set parameters mu0=0 sigma0=1 phi=1 cQ=1 cR=1 # kf = Kfilter1(50,y,A,mu0,sigma0,phi,0,0,cQ,cR,0) # commented out ks = Ksmooth1(50,y,A,mu0,sigma0,phi,0,0,cQ,cR,0) Time=1:num par(mfrow=c(3,1)) plot(Time, mu[-1], main="Prediction", ylim=c(-5,10)) lines(ks$xp) lines(ks$xp+2*sqrt(ks$Pp), lty="dashed", col="blue") lines(ks$xp-2*sqrt(ks$Pp), lty="dashed", col="blue") plot(Time, mu[-1], main="Filter", ylim=c(-5,10)) lines(Time, ks$xf) lines(Time, ks$xf+2*sqrt(ks$Pf), lty="dashed", col="blue") lines(Time, ks$xf-2*sqrt(ks$Pf), lty="dashed", col="blue") plot(Time, mu[-1], main="Smoother", ylim=c(-5,10)) lines(Time, ks$xs) lines(Time, ks$xs+2*sqrt(ks$Ps), lty="dashed", col="blue") lines(Time, ks$xs-2*sqrt(ks$Ps), lty="dashed", col="blue") # if you want to see the initial stuff use these lines: mu[1] # actual initial value ks$x0n # smoothed value sqrt(ks$P0n) # and its stan error