#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!# #-- assumes ss0 has been sourced --# #!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!# #-- 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] # set initial parameters mu0=0 sigma0=1 phi=1 cQ=1 # These are the Cholesky decomps of Q and R, cR=1 # which are the standard devs in this case. # you don't have to run the filter because Ksmooth0 returns # the filters and predictions as well as the smoothers... # ... that's why the kf= line below is commented out ... # it's just an example of what to do if you don't want to smooth kf = Kfilter0(50,y,1,mu0,sigma0,phi,cQ,cR) ks = Ksmooth0(50,y,1,mu0,sigma0,phi,cQ,cR) #-- pictures ## notes: ## mu_t^t-1=ks\$xp, P_t^t-1=ks\$Pp, t= 1,...,n and "p" for prediction ## mu_t^t=ks\$xf, P_t^t=ks\$Pf, t= 1,...,n and "f" for filter ## mu_t^n=ks\$xs, P_t^n=ks\$Ps, t= 1,...,n and "s" for smoother 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