*
* GARCHGIBBS.PRG
* Manual Example 13.5
*
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,nkeep=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
do draws=-nburn,nkeep
   *
   * 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
   }
   if draws<=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 draws
*
disp "Percentage of jumps for mean parms" 100.0*float(bjumps)/(nburn+nkeep+1)
disp "Percentage of jumps for GARCH parms" 100.0*float(gjumps)/(nburn+nkeep+1)
*
* Transform the accumulated first and second moments into mean and variance
*
compute b=b/nkeep,bxx=bxx/nkeep-%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)