*
* Volatility model from section 14.4
* Note - this can take quite a while (15-30 minutes) to run
*
open data sv.dat
data(format=free,skip=1) 1 945 xrate
*
* The data are already in log differenced form
*
set clipx = %if(abs(xrate)<.01,%tsign(.01,xrate),xrate)
*
nonlin psi1 psi2 psi3
compute psi3=log(.9/.1)
compute psi2=log(0.25)
compute psi1=log(0.25)
*
compute ndraws=200
set weight 1 ndraws = 0.0
*
find(method=simplex) max %logl
   *
   * Compute approximating Gaussian model
   *
   compute sigmasq=exp(2.0*psi1)
   compute phi=exp(psi3)/(1+exp(psi3))
   compute sigetasq=exp(2.0*psi2)
   set htilde = 2.0
   set ytilde = log(clipx**2/sigmasq)
   do iters=1,20
      dlm(a=phi,c=1.0,x0=0.0,sx0=sigetasq/(1-phi**2),$
         y=ytilde,sv=htilde,sw=sigetasq,type=smooth) / xstates
      set theta = xstates(t)(1)
      set htilde = 2*sigmasq*exp(theta)/clipx**2
      set ytilde = theta-.5*htilde+1
   end do iters
   *
   compute loglg=%logl
   *
   *  Use the same seed for all simulations
   *
   seed 23431
   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=phi,c=1.0,x0=0.0,sx0=sigetasq/(1-phi**2),$
            y=ytilde,sv=htilde,sw=sigetasq,type=csim) / xstates
         set xdraw = xstates(t)(1)
      }
      else
         set xdraw = 2*theta-xdraw
      *
      set truel  = %logdensity(sigmasq*exp(xdraw),xrate)
      set fakel  = %logdensity(htilde,ytilde-xdraw)
      sstats / 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