#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!# #-- 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