*
* Example 12.5 from page 586
*
open data m-sp5-6204.txt
calendar(m) 1962
data(format=free,org=columns) 1962:1 2004:11 sp500
*
set spret = sp500-sp500{1}
*
spgraph(footer="Figure 12.12 Time plots of monthly S&P 500 index",vfields=2)
graph(header="(a) log SP 500 index")
# sp500
graph(header="(b) returns")
# spret
spgraph(done)
*
compute meanx2=%digamma(0.5)-log(0.5)
compute varx2 =%trigamma(0.5)
*
nonlin(parmset=svparms) mu a1 lsigeps lsigeta
*
diff(center) spret / cret
compute mu=%mean
compute a1=.6,lsigeps=0.0,lsigeta=-2.0
*
dlm(parmset=svparms,method=bfgs,a=a1,c=1.0,y=log((spret-mu)**2)-meanx2-lsigeps,$
   sx0=exp(lsigeta)/(1-a1**2),sv=varx2,sw=exp(lsigeta)) 2 *
compute logl0=%logl
compute sxx=%decomp(%xx)
compute betax=%beta
dlm(type=smooth,vhat=vhat,a=a1,c=1.0,y=log((spret-mu)**2)-meanx2-lsigeps,$
   sx0=exp(lsigeta)/(1-a1**2),sv=varx2,sw=exp(lsigeta)) 2 *
set gap = .5*(vhat(t)(1)+meanx2)-exp((vhat(t)(1)+meanx2))/2.0-log(2.5)-%logdensity(varx2,vhat(t)(1))
sstats 2 * gap>>logpoverq0
*
compute ndraws=1000
*
dec vect[series] stats(4)
do i=1,4
   set stats(i) 1 ndraws = 0.0
end do i
set weight 1 ndraws = 0.0
*
dec vect udraw(4)
*
do draws=1,ndraws
   compute udraw=%ran(1.0)
   compute loggdensity=-.5*%normsqr(udraw)
   compute %parmspoke(svparms,betax+sxx*udraw)

   compute stats(1)(draws)=mu
   compute stats(2)(draws)=a1
   compute stats(3)(draws)=lsigeps
   compute stats(4)(draws)=lsigeta
   if abs(a1)>=1.0 {
      compute weight(draws)=0.0
      next
   }
   dlm(method=solve,$
      a=a1,c=1.0,y=log((spret-mu)**2)-meanx2-lsigeps,sx0=exp(lsigeta)/(1-a1**2),sv=varx2,sw=exp(lsigeta))
   compute logglike=%logl-logl0

   dlm(type=csim,$
      a=a1,c=1.0,y=log((spret-mu)**2)-meanx2-lsigeps,sx0=exp(lsigeta)/(1-a1**2),sv=varx2,sw=exp(lsigeta),$
      vhat=vhat)
   set gap = .5*(vhat(t)(1)+meanx2)-exp((vhat(t)(1)+meanx2))/2.0-log(2.5)-%logdensity(varx2,vhat(t)(1))
   sstats 2 * gap>>logpoverq
   compute weight(draws)=exp(logpoverq-logpoverq0+logglike-loggdensity)
   if weight(draws)>1.e+3
      break
end do draws