* * GARCHGIBBS.RPF * Random walk Metropolis-Hastings for GARCH model * * RATS User's Guide, Example from Section 16.8. * 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)