*
* Boat race example from section 14.5
* (Note-this examples takes several minutes to run).
*
open data boat.dat
cal 1829
data(format=free,skip=1,missing=-9999.99) 1829:1 2000:1 cambridge
*
dec series p theta
*
* The likelihood is quite flat as a function of sigsq. Even with a very large
* number of draws, the approximation error in computing the E(p/g) is large
* enough relative to the likelihood that the estimate of sigsq can vary by about
* +/-.02
*
compute ndraws=4000
set weight 1 ndraws = 0.0
*
nonlin sigsq
compute sigsq=1.0
*
frml ctest    = p*(1-p)
frml svtest   = p*(1-p)
frml yadjtest = (cambridge-p)+p*(1-p)*theta
*
find(method=bfgs) max %logl
   if sigsq<0.0 {
      compute %logl=%na
      next
   }
   *
   *  Compute the approximating Gaussian model
   *
   set theta = 0.0
   do iters=1,10
      set p     = exp(theta)/(1+exp(theta))
      set adj   = yadjtest
      set c     = ctest
      set sv    = svtest
      dlm(a=1.0,c=c,sv=sv,y=adj,sw=sigsq,exact,type=smooth) / xstates
      set theta = xstates(t)(1)
   end do iters
   compute loglg=%logl
   *
   *  Use the same seed for all simulations
   *
   seed 56511
   do draws=1,ndraws
      *
      * Do a conditional draw on the odd draws and an additive antithete on
      * the even ones
      *
      if %clock(draws,2)==1 {
         dlm(a=1.0,c=c,sv=sv,y=adj,sw=sigsq,exact,type=csim) / xstates
         set xdraw = xstates(t)(1)
      }
      else
         set xdraw = 2*theta-xdraw
      set phat   = exp(xdraw)/(1+exp(xdraw))
      set truel  = cambridge*log(phat)+(1-cambridge)*log(1-phat)
      set fakel  = %logdensity(sv,adj-c*xdraw)
      sstats(smpl=%valid(cambridge)) / fakel>>fakelogl truel>>truelogl
      *
      * This is kept in log scale for now
      *
      compute weight(draws)=truelogl-fakelogl
   end do draws
   *
   * Compute the mean of the difference in log likelihoods. Subtract that
   * prior to exp'ing to avoid problems with overflows. log(wbar)+center
   * will be the correction term for the log likelihood
   *
   sstats(mean) 1 ndraws weight>>center
   sstats(mean) 1 ndraws exp(weight-center)>>wbar
   compute %logl=loglg+log(wbar)+center
end find
*
*  Do Kalman smoothed estimates at the final set of parameters
*
dlm(a=1.0,c=c,sv=sv,y=adj,sw=sigsq,exact,type=smooth) / xstates vstates
set theta  = xstates(t)(1)
set lower  = theta+%invnormal(.25)*sqrt(vstates(t)(1,1))
set upper  = theta+%invnormal(.75)*sqrt(vstates(t)(1,1))
set phat   = exp(theta)/(1+exp(theta))
set plower = exp(lower)/(1+exp(lower))
set pupper = exp(upper)/(1+exp(upper))
graph(footer="Fig 14.7 Results (dots) with predicted probabilities and 50% confidence intervals",$
  overlay=dots,ovsame) 4
# phat
# plower / 2
# pupper / 2
# cambridge / 1