*
* Durbin & Koopman
* Time Series Analysis by State Space Methods, 2nd edition.
* Illustration from section 12.4.4
* Bootstrap particle filter
*
open data nile.dat
calendar 1871
data(format=free,org=columns,skips=1) 1871:1 1970:1 nile
*
* Standard Kalman filter
*
dlm(a=1.0,c=1.0,x0=0.0,sx0=1.e+7,sv=15099.0,sw=1469.1,presample=x1,$
   y=nile,vhat=vhat,svhat=fhat) / xstates vstates
*
* Get the KF estimates of the state and variance
*
set ahatkf = xstates(t)(1)
set phatkf = vstates(t)(1,1)
*
* Bootstrap filter without resampling
*
compute sigmaeta=sqrt(1469.1)
compute sigmax0 =sqrt(1.e+7)
compute sigmaepssq=15099.0
*
* Set up the particles
*
compute nstate=1
compute npart=10000
*
* This can be done without saving the entire histories of alpha and logw
* - just keeping the vectors from the last time period.
*
dec vect[vect] alphanew(npart)
dec series[vect[vect]] alpha
ewise alphanew(i)=%zeros(nstate,1)
gset alpha   = alphanew
*
dec series[vect] logw
dec vect wt(npart)
gset logw    = %zeros(npart,1)
*
dec series[vect] alphahat
dec series[symm] alphahatxx
*
gset alphahat   = %zeros(nstate,1)
gset alphahatxx = %zeros(nstate,nstate)
*
set  ess        = 0.0
*
do time=1871:1,1970:1
   do part=1,npart
      if time==1871:1 {
         compute alpha(time)(part)=sigmax0*%ranmat(1,1)
         compute laggedlogw=0.0
      }
      else {
         compute alpha(time)(part)=alpha(time-1)(part)+%ran(sigmaeta)
         compute laggedlogw=logw(time-1)(part)
      }
      compute logw(time)(part)=laggedlogw+$
          %logdensity(sigmaepssq,nile(time)-alpha(time)(part)(1))
   end do part
   *
   * Do overflow-safe normalization of weights
   *
   compute maxlogw=%maxvalue(logw(time))
   compute logw(time)=logw(time)-maxlogw
   compute wt=%exp(logw(time))
   compute wt=wt/%sum(wt)
   *
   do part=1,npart
      compute alphahat(time)=alphahat(time)+wt(part)*alpha(time)(part)
      compute alphahatxx(time)=alphahatxx(time)+wt(part)*%outerxx(alpha(time)(part))
   end do part
   compute ess(time)=1.0/%normsqr(wt)
end do time
*
set phat    = alphahatxx(t)(1,1)-alphahat(t)(1)^2
set ahatpf  = alphahat(t)(1)
set lowerpf = ahatpf+sqrt(phat)*%invnormal(.05)
set upperpf = ahatpf+sqrt(phat)*%invnormal(.95)
*
spgraph(vfields=2,hfields=2,$
   footer="Figure 12.1 Bootstrap filter with N=10000 without resampling")
 graph(footer="Filtered Estimates, 90% Band and Data",overlay=dots,ovsame) 4
 # ahatpf
 # lowerpf / 2
 # upperpf / 2
 # nile
 graph(row=1,col=2,footer="Kalman filter and Particle filter estimates") 2
 # ahatkf
 # ahatpf
 graph(row=2,col=1,footer="Kalman filter and Particle filter variances") 2
 # phatkf
 # phat
 graph(row=2,col=2,footer="Effective Sample Size") 1
 # ess
spgraph(done)
*
* With resampling
*
* Clear the accumulators
*
gset alphahat   = %zeros(nstate,1)
gset alphahatxx = %zeros(nstate,nstate)
set  ess        = 0.0
*
do time=1871:1,1970:1
   do part=1,npart
      if time==1871:1 {
         compute alpha(time)(part)=sigmax0*%ranmat(1,1)
         compute laggedlogw=0.0
      }
      else {
         compute alpha(time)(part)=alpha(time-1)(part)+%ran(sigmaeta)
         compute laggedlogw=logw(time-1)(part)
      }
      compute logw(time)(part)=laggedlogw+$
          %logdensity(sigmaepssq,nile(time)-alpha(time)(part)(1))
   end do part
   compute maxlogw=%maxvalue(logw(time))
   compute logw(time)=logw(time)-maxlogw
   compute wt=%exp(logw(time))
   compute wt=wt/%sum(wt)
   *
   do part=1,npart
      compute alphahat(time)=alphahat(time)+wt(part)*alpha(time)(part)
      compute alphahatxx(time)=alphahatxx(time)+wt(part)*%outerxx(alpha(time)(part))
   end do part
   compute ess(time)=1.0/%normsqr(wt)
   *
   * Resample, using %RANBRANCH to stratify based upon the current
   * weights.
   *
   do i=1,npart
      compute redraw=%ranbranch(wt)
      compute alphanew(i)=alpha(time)(redraw)
   end do i
   compute alpha(time)=alphanew
   *
   * Reset the probabilities to be even across particles
   *
   compute logw(time)=%zeros(npart,1)
end do time
*
set phat    = alphahatxx(t)(1,1)-alphahat(t)(1)^2
set ahatpf  = alphahat(t)(1)
set lowerpf = ahatpf+sqrt(phat)*%invnormal(.05)
set upperpf = ahatpf+sqrt(phat)*%invnormal(.95)
*
spgraph(vfields=2,hfields=2,$
   footer="Figure 12.2 Bootstrap filter with N=10000 with resampling")
 graph(footer="Filtered Estimates, 90% Band and Data",overlay=dots,ovsame) 4
 # ahatpf
 # lowerpf / 2
 # upperpf / 2
 # nile
 graph(row=1,col=2,footer="Kalman filter and Particle filter estimates") 2
 # ahatkf
 # ahatpf
 graph(row=2,col=1,footer="Kalman filter and Particle filter variances") 2
 # phatkf
 # phat
 graph(row=2,col=2,footer="Effective Sample Size") 1
 # ess
spgraph(done)