GARCHGIBBS.RPF

GARCHGIBBS.RPF is an example of Metropolis-Hastings for a univariate GARCH model, estimating the same GARCH model used in the GARCHIMPORT.RPF example. This blocks the parameters into the mean AR(1) parameters and the GARCH parameters. To illustrate both methods, the AR(1) parameters are handled with a fixed proposal and the GARCH are done by Random Walk Metropolis.

Estimate a linear regression to get its residual variance for use as the pre-sample value for the GARCH models.

linreg ftbs3

# constant ftbs3{1}

compute h0 = %sigmasq

Now, estimate the 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. (GX will be used throughout to save the density at the current value of the draw for the regression coefficients).

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)

This sets the control parameters for the Gibbs sampler:

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)

These are some declarations for VECTORs that are used in the analysis.

dec vect betau(%nreg) ygarch

Initialize counters for the number of jumps in the two part of the sampler. (BJUMPS for the regression coefficients and GJUMPS for the GARCH coefficients).

compute bjumps=0, gjumps=0

The draw loop includes an INFOBOX to keep track of the progress.

infobox(action=define,progress,lower=-nburn,upper=ndraws) $

"Random Walk Metropolis"

The pseudo-code for the sampling loop is the following:

do draw=-nburn,ndraws

...draw regression parameters given the current GARCH parameters using a fixed proposal distribution.

...draw GARCH parameters given the current regression parameters using a random walk proposal

...if the value of draws>0, update the first and second moment accumulators

end do draw

The three parts to this are shown in detail below.

Drawing Regression Parameters Given The Current GARCH Parameters

First, evaluate the log likelihood (into FX) at the current values for XBETA and XGARCH.

garch(p=1,q=1,reg,method=eval,presample=h0,$

initial=xbeta~~xgarch) / ftbs

# constant ftbs3{1}

set u = %resids

compute fx = %logl

Next, draw a test vector (YBETA) for the regression block from the (fixed) proposal distribution. Fetch the log density at the draw (into GY). Recalculate the log likelihood (with the existing XGARCH, but YBETA replacing XBETA) (into FY)

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

Finally, compute the jump alpha. If we need to move, reset XBETA and GX, and copy the residuals into the series U.

compute alpha=exp(fy-fx+gy-gx)

if %ranflip(alpha) {

compute bjumps = bjumps + 1

compute xbeta=ybeta, gx=gy

set u = %resids

}

Draw GARCH Parameters Given the Current Regression Parameters

Because we are taking the regression coefficients (that is, the mean model) as given, we can evaluate the GARCH model on the residuals only, using the NOMEAN option. First, 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

Next, 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). If we need to move, save the YGARCH draw to XGARCH and update the jump counter.

compute alpha = %if(%valid(fy),exp(fy-fx),0.0)

if %ranflip(alpha) {

compute gjumps=gjumps+1

compute xgarch=ygarch

}

Update the accumulators for the moments

Glue the current values of the two blocks together into a single VECTOR and add the value and its %OUTERXX onto the running sums.

compute betau=xbeta~~xgarch

compute b=b+betau,bxx=bxx+%outerxx(betau)

After the Draw Loop

This first displays the percentage of jumps made by each of the block samplers. We would hope that the

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)

Then transform the accumulated first and second moments into mean and variance, and display the means and standard deviations in a REPORT.

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)

Full Program

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) 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 into u.

compute alpha=exp(fy-fx+gy-gx)

if %ranflip(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 %ranflip(alpha) {

compute gjumps=gjumps+1

compute xgarch=ygarch

}

infobox(current=draw)

if draw<=0

* 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)

Output

The final part of this depends upon random numbers, and so will not be exactly reproducible.

Linear Regression - Estimation by Least Squares

Dependent Variable FTBS3

Monthly Data From 1957:02 To 2006:12

Usable Observations 599

Degrees of Freedom 597

Centered R^2 0.9719842

R-Bar^2 0.9719373

Uncentered R^2 0.9941920

Mean of Dependent Variable 5.3828881469

Std Error of Dependent Variable 2.7551056598

Standard Error of Estimate 0.4615335100

Sum of Squared Residuals 127.16886894

Regression F(1,597) 20712.4001

Significance Level of F 0.0000000

Log Likelihood -385.7953

Durbin-Watson Statistic 1.3186

Variable Coeff Std Error T-Stat Signif

************************************************************************************

1. Constant 0.0816496008 0.0413816505 1.97309 0.04894649

2. FTBS3{1} 0.9853633822 0.0068466985 143.91803 0.00000000

GARCH Model - Estimation by BFGS

Convergence in 41 Iterations. Final criterion was 0.0000000 <= 0.0000100

Dependent Variable FTBS3

Monthly Data From 1957:02 To 2006:12

Usable Observations 599

Log Likelihood -61.3428

Variable Coeff Std Error T-Stat Signif

************************************************************************************

1. Constant 0.0631210293 0.0237229350 2.66076 0.00779646

2. FTBS3{1} 0.9900678792 0.0053806857 184.00404 0.00000000

3. C 0.0014556841 0.0005784613 2.51648 0.01185349

4. A 0.3065947388 0.0461751344 6.63982 0.00000000

5. B 0.7300827508 0.0325249623 22.44684 0.00000000

Percentage of jumps for mean parms 47.65930

Percentage of jumps for GARCH parms 47.83202

a 0.063282 0.013174

b 0.990070 0.002992

gamma 0.001803 0.000694

alpha 0.326477 0.051527

beta 0.715309 0.035293