*
* GARCHGIBBS.RPF
* RATS Version 8, User's Guide, Example 16.5
* Random walk MH for GARCH model
*
open data haversample.rat
calendar(m) 1957
data(format=rats) 1957:1 2006:12 ftbs3
*
* Estimate a linear regression to get its residual variance for use as
* the pre-sample value.
*
linreg ftbs3
# constant ftbs3{1}
compute h0=%sigmasq
*
* Estimate a GARCH
*
garch(p=1,q=1,reg,presample=h0,hseries=h) / ftbs3
# constant ftbs3{1}
*
* Pull out the regression coefficients (1 and 2). This forms block one
* of the Metropolis/Gibbs. Get a decomp of its submatrix of the
* covariance matrix for use in the proposal density. Do a draw from this
* to initialize the Gibbs sampler and save the log kernel of the density
* at the draw.
*
compute xbeta0=%xsubvec(%beta,1,2)
compute fbeta =%decomp(%xsubmat(%xx,1,2,1,2))
compute xbeta =xbeta0+%ranmvnormal(fbeta)
compute gx    =%ranlogkernel()
*
* Pull out the GARCH coefficients (3,4 and 5). This forms block two.
* Again, get the decomp of its covariance matrix. Because these are
* being done by Random Walk Metropolis, we don't need the proposal
* densities.
*
compute fgarch=%decomp(%xsubmat(%xx,3,5,3,5))
compute xgarch=%xsubvec(%beta,3,5)
*
compute nburn=1000,ndraws=10000
*
* Initialize vectors for the first and second moments of the coefficients
*
compute sumwt=0.0,sumwt2=0.0
compute [vect] b=%zeros(%nreg,1)
compute [symm] bxx=%zeros(%nreg,%nreg)
dec vect betau(%nreg) betadraw ygarch
*
* Counters for the number of jumps
*
compute bjumps=0,gjumps=0
infobox(action=define,progress,lower=-nburn,upper=ndraws) $
  "Random Walk Metropolis"
do draw=-nburn,ndraws
   *
   * Drawing regression parameters. Evaluate the log likelihood (into
   * FX) at the values for XBETA and XGARCH.
   *
   garch(p=1,q=1,reg,method=eval,presample=h0,$
      initial=xbeta~~xgarch) / ftbs3
   # constant ftbs3{1}
   set u = %resids
   compute fx=%logl
   *
   * Draw a test vector from the (fixed) proposal distribution.
   * Recalculate the log likelihood (into FY) and fetch the log density
   * at the draw (into GY).
   *
   compute ybeta=xbeta0+%ranmvnormal(fbeta)
   compute gy=%ranlogkernel()
   garch(p=1,q=1,reg,method=eval,presample=h0,initial=ybeta~~xgarch) / ftbs3
   # constant ftbs3{1}
   compute fy=%logl
   *
   * Compute the jump alpha. If we need to move, reset xbeta and gx, and
   * copy the residuals from Y into u.
   *
   compute alpha=exp(fy-fx+gy-gx)
   if %uniform(0.0,1.0)<alpha {
      compute bjumps=bjumps+1
      compute xbeta=ybeta,gx=gy
      set u = %resids
   }
   *
   * Evaluate the log likelihood of a GARCH model (into FX) on the
   * residuals at the current settings for XGARCH.
   *
   garch(p=1,q=1,nomean,method=eval,presample=h0,initial=xgarch) / u
   compute fx=%logl
   *
   * Draw from the proposal distribution centered around XGARCH and
   * evaluate the log likelihood into FY.
   *
   compute ygarch=xgarch+%ranmvnormal(fgarch)
   garch(p=1,q=1,nomean,method=eval,presample=h0,initial=ygarch) / u
   compute fy=%logl
   *
   * Compute the jump alpha. (Because draws for the GARCH parameters can
   * turn the model explosive, check for an invalid FY).
   *
   compute alpha=%if(%valid(fy),exp(fy-fx),0.0)
   if %uniform(0.0,1.0)<alpha {
      compute gjumps=gjumps+1
      compute xgarch=ygarch
   }
   infobox(current=draw)
   if draw<=0
      next
   *
   * Update the accumulators if we're past the burn-in
   *
   compute betau=xbeta~~xgarch
   compute b=b+betau,bxx=bxx+%outerxx(betau)
end do draw
infobox(action=remove)
*
disp "Percentage of jumps for mean parms" $
  100.0*float(bjumps)/(nburn+ndraws+1)
disp "Percentage of jumps for GARCH parms" $
  100.0*float(gjumps)/(nburn+ndraws+1)
*
* Transform the accumulated first and second moments into mean and variance
*
compute b=b/ndraws,bxx=bxx/ndraws-%outerxx(b)
report(action=define)
report(atrow=1,atcol=1,fillby=cols) "a" "b" "gamma" "alpha" "beta"
report(atrow=1,atcol=2,fillby=cols) b
report(atrow=1,atcol=3,fillby=cols) %sqrt(%xdiag(bxx))
report(action=show)