*
* Example 12.3 from pp 569-571
* (Note - even though we've reduced the number of kept iterations to 2000 from
* 5000, this still takes a long time to run)
*
open data m-sp6299.dat
calendar(m) 1962
data(format=free,org=columns) 1962:1 1999:12 sp500
*
graph(footer="Figure 12.3 Monthly log returns of the S&P 500 index")
# sp500
*
* GARCH model
*
garch(p=1,q=1,hseries=garchvol,resids=u) / sp500
set ustd = u/sqrt(garchvol)
set ustdsq = ustd**2
@regcorrs(title="Standardized Residuals") ustd
@regcorrs(dfc=2,title="Standardized Residuals Squared") ustdsq
*
* SV model
*
compute nkeep=2000
compute nburn=100
compute ngrid=500
*
* Set up the grid for drawing h's. (This allows for a quite a bit wider band than
* shown in the book).
*
stats sp500
dec vect hgrid(ngrid) fgrid(ngrid)
ewise hgrid(i)=i*(20.0*%variance)/ngrid
compute mu=%mean
*
* Start with the variances from the GARCH model
*
set h = garchvol
*
* Prior precision (inverse of covariance matrix) and mean for the alpha's
*
compute [symm] halpha0=%diag(||1.0/.09,1.0/.04||)
compute [vect] alpha0=||0.4,0.8||
dec symm valpha
dec vect malpha dalpha
*
* Prior precision and mean for mu
*
compute [symm] hmu0=||1.0/9.0||
compute [vect] mu0 =||0.0||
dec symm vmu
dec vect mmu
*
* Prior for the sigma squared. pscale is the "prior mean" for sigma squared
* (actually the reciprocal of mean of the density for 1/sigma-squared) and pdf is
* the prior degrees of freedom for the chi-square.
*
compute pscale=.2,pdf=5.0
compute sigmasq=.2
*
compute gstart=%regstart(),gend=%regend()
*
* Initialize bookkeeping. We're keeping four complete sets (for mu, a0, a1 and sigmasq)
* and the sum of the h's series
*
dec vect[series] stats(4)
do i=1,4
   set stats(i) 1 nkeep = 0.0
end do i
set hstats gstart gend = 0.0
*
infobox(action=define,progress,lower=-nburn+1,upper=nkeep) "Gibbs Sampler"
do draws=-nburn+1,nkeep
   infobox(current=draws)
   *
   *  Draw alphas given other information
   *
   set z = log(h)
   cmom(noprint)
   # constant z{1} z
   *
   *  Step one - compute the inverse of the sums of the precision matrices. The
   *  precision for the data is sigma**-2 X'X (inverse of the standard OLS
   *  covariance matrix).
   *
   compute valpha=inv(halpha0+%xsubmat(%cmom,1,2,1,2)/sigmasq)
   *
   *  Step two - compute the mean of the posterior.
   *
   compute malpha=valpha*(halpha0*alpha0+%xsubmat(%cmom,1,2,3,3)/sigmasq)
   *
   *  Make a random draw from the Normal distribution for the alpha's.
   *
   compute dalpha=malpha+%ranmvnormal(%decomp(valpha))
   compute a0=dalpha(1),a1=dalpha(2)
   *
   * Draw sigma's given alpha's. Compute the sum of squares of the residuals,
   * and draw from the appropriate inverse chi-squared density.
   *
   set u = z-a0-a1*z{1}
   sstats gstart+1 gend u**2>>%rss
   compute sigmasq=(pscale*pdf+%rss)/%ranchisqr(%nobs+pdf)
   *
   * Draw h's given everything else
   *
   do time=gstart,gend
      if time==gstart
         compute meanlogh=a0+a1*log(h(time+1)),sigmalogh=sigmasq
      else
      if time==gend
         compute meanlogh=a0+a1*log(h(time-1)),sigmalogh=sigmasq
      else {
         compute meanlogh  = (a0*(1-a1)+a1*(log(h(time+1))+log(h(time-1))))/(1+a1**2)
         compute sigmalogh = sigmasq/(1+a1**2)
      }
      ewise fgrid(i)=hgrid(i)**(-1.5)*exp(-(sp500(time)-mu)**2/(2*hgrid(i))-(log(hgrid(i))-meanlogh)**2/(2*sigmalogh))
      compute h(time)=%rangrid(hgrid,fgrid)
   end do time
   *
   * Draw mu's given h's
   *
   cmom(spread=h)
   # constant sp500
   compute vmu=inv(hmu0+%xsubmat(%cmom,1,1,1,1))
   compute mmu=vmu*(hmu0*mu0+%xsubmat(%cmom,1,1,2,2))
   compute mu=mmu(1)+%ran(sqrt(vmu(1,1)))

   if draws<=0
      next

   compute stats(1)(draws)=mu
   compute stats(2)(draws)=a0
   compute stats(3)(draws)=a1
   compute stats(4)(draws)=sigmasq

   set hstats gstart gend = hstats+h

end do draws
infobox(action=remove)
*
stats stats(1) 1 nkeep
stats stats(2) 1 nkeep
stats stats(3) 1 nkeep
stats stats(4) 1 nkeep
set hstats gstart gend = hstats/nkeep
*
spgraph(vfields=2,footer="Figure 12.5 Time plots of fitted volatilies")
graph(header="(a) Stochastic Volatility model",min=0.0)
# hstats
graph(header="(b) GARCH model",min=0.0)
# garchvol
spgraph(done)