*
* Example 12.2 from page 564
* Note: at 1000 draws, this takes a LONG time. Depending upon machine speed
* it might take over an hour.
*
cal(w) 1962:1:5
all 1999:9:10
open data wgs3yr.dat
data(org=columns,format=free) / date rate3yr
*
set c3 = rate3yr-rate3yr{1}
*
* The estimates shown in the text were done by full sample maximum likelihood,
* although the remaining analysis is done conditionally on the pre-sample values.
* We'll use the second set of estimates (which can be done with LINREG) to
* initialize the Gibbs sampler.
*
boxjenk(noconst,ar=3,maxl) c3 1988:3:18 *
linreg c3 1988:3:18 *
# c3{1 to 3}
*
* Prior for deltas
*
compute g1=5.0,g2=95.0
*
* Prior for outliers
*
compute xisq=0.1
*
* Prior precision for phi's
*
dec symm hphi0 vphi
dec vect phim  phidraw
compute hphi0=4.0*%identity(3)
*
* Prior for sigmasq
*
compute pscale=.00256,pdf=5
*
compute phi1=%beta(1),phi2=%beta(2),phi3=%beta(3)
compute eps=.05,sigmasq=%sigmasq
*
set delta 1988:3:18 * = 0.0
set beta  1988:3:18 * = 0.0
*
*  The analysis is all done conditional on the three pre-sample values for c3
*
compute t0=1988:3:18
compute [vector] x0=||c3(t0-1),c3(t0-2),c3(t0-3)||
set ytilde = c3
*
compute ndraws=250
compute burnin=50
*
* Prepare bookkeeping locations. We'll keep five full reports (for the three
* phi's, sigmasq and eps) and just the sums for the deltas and betas.
*
dec vect[series] stats(5)
do i=1,5
   set stats(i) 1 ndraws = 0.0
end do i
set deltahat t0 * = 0.0
set betahat  t0 * = 0.0
*
infobox(action=define,progress,lower=1,upper=ndraws)
do draws=1,ndraws
   *
   * Draw betas. This is equivalent to the procedure described in the book. The delta x beta
   * term is treated as a measurement error in a state space model. When we do a conditional
   * simulation of the model, we'll get simulated values for the betas for the time periods
   * where delta is 1. For the time periods where delta is zero, we just need to do a Normal
   * draw.
   *
   dlm(start=%(a=||phi1,phi2,phi3|1.0,0.0,0.0|0.0,1.0,0.0||,sw=%diag(||sigmasq,0.0,0.0||)),$
      a=a,c=||1.0,0.0,0.0||,y=c3,sw=sw,x0=x0,sv=delta*xisq*100,type=csimulate,vhat=betav) t0 *
   set beta 1988:3:18 * = %if(delta,betav(t)(1),%ran(sqrt(xisq)))
   *
   * Draw deltas. This is computed by computed the log likelihood of the state space
   * model with delta(time) tested at both 0 and 1. delta(time) is then randomly
   * set to 0 or 1 based upon the values of eps x f given 1 and 1-eps x f given 0
   * %ranbranch returns 1 or 2, so we have to subtract 1 from it to get the simulated
   * delta.
   *
   do time=t0,1999:9:10
      compute delta(time)=1.0
      dlm(start=%(a=||phi1,phi2,phi3|1.0,0.0,0.0|0.0,1.0,0.0||,sw=%diag(||sigmasq,0.0,0.0||)),$
         a=a,c=||1.0,0.0,0.0||,y=c3-delta*beta,sw=sw,x0=x0,type=filter) t0 *
      compute logl1=%logl
      compute delta(time)=0.0
      dlm(start=%(a=||phi1,phi2,phi3|1.0,0.0,0.0|0.0,1.0,0.0||,sw=%diag(||sigmasq,0.0,0.0||)),$
         a=a,c=||1.0,0.0,0.0||,y=c3-delta*beta,sw=sw,x0=x0,type=filter) t0 *
      compute logl2=%logl
      compute delta(time)=%ranbranch(||(1-eps)*exp(logl2),eps*exp(logl1)||) - 1
   end do time
   *
   * Compute y adjusted for the outlier components. Drawing a phi is standard from that.
   *
   set ytilde t0 * = c3-delta*beta
   *
   cmom(noprint) t0 *
   # ytilde{1 2 3 0}
   compute vphi=inv(hphi0+%xsubmat(%cmom,1,3,1,3)/sigmasq)
   compute phim=vphi*%xsubmat(%cmom,1,3,4,4)/sigmasq
   compute phidraw=phim+%ranmvnormal(%decomp(vphi))
   compute phi1=phidraw(1),phi2=phidraw(2),phi3=phidraw(3)
   set u t0 * = ytilde-phi1*ytilde{1}-phi2*ytilde{2}-phi3*ytilde{3}
   *
   * Draw new values for sigmasq and eps
   *
   sstats t0 * u**2>>%rss delta>>count
   compute sigmasq=(pscale*pdf+%rss)/%ranchisqr(%nobs+pdf)
   *
   * The draw for eps isn't described in the text. Given the prior and the data,
   * the posterior for eps is Beta(g1+count,g2+# of observations-count).
   *
   compute eps=%ranbeta(g1+count,g2+%nobs-count)
   *
   * Save the statistics of interest
   *
   compute stats(1)(draws)=phi1
   compute stats(2)(draws)=phi2
   compute stats(3)(draws)=phi3
   compute stats(4)(draws)=sigmasq
   compute stats(5)(draws)=eps
   *
   if draws>burnin {
      set deltahat t0 * = deltahat+delta
      set betahat  t0 * = betahat+beta
   }
   infobox(current=draws)
end do draws
infobox(action=remove)
*
set deltahat t0 * = deltahat/(ndraws-burnin)
set betahat  t0 * = betahat/(ndraws-burnin)
*
spgraph(footer="Figure 12.2 Time plots of weekly change series of U.S. 3-year Treasury constant maturity interest rate",vfields=3)
graph(header="(a) Changes in weekly interest rate")
# c3 t0 *
graph(header="(b) Posterior probability of being an outlier",min=0.0,max=1.0,style=spike)
# deltahat
graph(header="(c) Posterior mean of outlier size",style=spike)
# betahat
spgraph(done)
*
report(action=define)
report(atrow=1,atcol=2) "Mean" "Std Dev"
stats(noprint) stats(1) burnin+1 ndraws
report(atrow=2,atcol=1) "phi1" %mean sqrt(%variance)
stats(noprint) stats(2) burnin+1 ndraws
report(atrow=3,atcol=1) "phi2" %mean sqrt(%variance)
stats(noprint) stats(3) burnin+1 ndraws
report(atrow=4,atcol=1) "phi3" %mean sqrt(%variance)
stats(noprint) stats(4) burnin+1 ndraws
report(atrow=5,atcol=1) "sigmasq" %mean sqrt(%variance)
report(action=show)