GIBBSVARBUILD.RPF

GibbsVARBuild.rpf does Gibbs sampling for a Bayesian VAR. It is effectively identical to the (older) GibbsVAR.rpf except that it uses the @BVARBuildPriorMN procedure to create the prior mean and precision matrices, rather than creating them using a (long) series of calculations.

Gibbs sampling analyzes the behavior of the VAR as a whole rather than one equation at a time, as is done with ESTIMATE. This can be applied successfully to VAR's of modest size. This example has 4 variables x 4 lags (plus constant), so it has a total of 68 regression parameters (4 equations with 17 coefficients in each). This will probably run in a reasonable time with up to perhaps 500 total parameters. The sticking point is that a full-system draw for the VAR coefficients requires taking a factor of a # of coefficients x # of coefficients matrix with no helpful structure to reduce the calculation time: factoring a matrix is an \(O(size^3)\) calculation, so a size 1000 matrix will take roughly eight times as long as 500.

The control parameters are shown at the top of the program:

compute lags=4 ;*Number of lags

compute nstep=16 ;*Number of response steps

compute nburn=500 ;*Number of burn-in draws

compute ndraws=2500 ;*Number of keeper draws

These are the number of lags in the VAR, the number of forecast steps, the number of "burn-in" Gibbs sampler draws and the number of retained Gibbs sampler draws.

The model itself is a four variable VAR on logs of GDP, investment and consumption with a short-term interest rate.

open data haversample.rat

cal(q) 1959

data(format=rats) 1959:1 2006:4 ftb3 gdph ih cbhm

set loggdp = log(gdph)

set loginv = log(ih)

set logc = log(cbhm)

These are the controlling parameters for the symmetric "Minnesota" prior—the own lag tightness and the relative tightness on the other lags.

compute tight=.1

compute other=.5

We set up the (non-Bayesian) form of the VAR

system(model=varmodel)

variables loggdp loginv logc ftb3

lags 1 to lags

det constant

end(system)

This uses @BVARBuildPriorMN and @BVARFinishPrior to create the mean and precision matrices needed for Gibbs sampling on a multivariate linear model.

@BVARBuildPriorMN(model=varmodel,tight=tight,other=other) hbpriors hpriors

@BVARFinishPrior hbpriors hpriors bprior hprior

The Gibbs sampler (which will alternate between drawing the residual covariance matrix given the VAR coefficients, and drawing the VAR coefficients given the residual covarianc matrix) will start at the least squares estimates. The CMOM instruction is used to compare the cross product matrix of the regressors, which substantially reduces the calculation time. (The fact that the regressors are the same in each equation considerably simiplifies this). Simulations of distributions in multivariate linear models are covered in detail in the Bayesian Methods course.

estimate(ols,sigma)

compute nvar=%nvar

compute sigmad=%sigma

cmom(model=varmodel)

The program generates out-of-sample forecasts for the four series of the model, generating both the means and standard errors. Unlike the FORECAST instruction (which treats the input model as known), this takes into account that the coefficients of the model have to be estimated. The following initializes the series used for collecting statistics on the generated forecasts (over the range FSTART to FEND which are the first NSTEP steps outside of sample.

dec vect[series] forecast(nvar)

dec vect[series] forestderr(nvar)

compute fstart=%regend()+1

compute fend =fstart+nstep-1

do i=1,nvar

set forecast(i) fstart fend = 0.0

set forestderr(i) fstart fend = 0.0

end do i

This is the pseudo-code for the Gibbs draw loop:

infobox(action=define,progress,lower=-nburn,upper=ndraws) "Gibbs Sampler"

do draw=-nburn,ndraws

infobox(current=draw)

Draw b given sigma

Draw sigma given b

if draw<=0

simulate out-of-sample data given current coefficients

end do draw

The steps within this are described below.

Draw b given sigma:

This takes the current draw for the covariance matrix (SIGMAD) and the cross product matrix of the regressors (%CMOM) combined with the precision matrix (HPRIOR) and mean (BPRIOR) of the prior and returns a draw from the multivariate Normal regression (a VAR is simply a special case of a multivariate regression with identical explanatory variable). Note that this is the step in the process which is most time consuming because it has to take a factor of a large (effectively) unstructured covariance matrix.

RSSMAT is the sums of squared residuals for the model given the new draw for the coefficients (BDRAW).

compute bdraw =%ranmvkroncmom(%cmom,inv(sigmad),hprior,bprior)

compute rssmat=%sigmacmom(%cmom,bdraw)

Draw sigma given b:

This is a straightforward inverse Wishart draw given the sums of squares just computed.

compute sigmad=%ranwisharti(%decomp(inv(rssmat)),%nobs)

Simulate out-of-sample data given current coefficients

We need to "poke" the new coefficients into the model, and use the draw for SIGMAD to control the random simulations of the data. We then accumulate the simulated values (into forecast(i)) and their squares (into forestderrs(i)). For now those are the sums, which we will convert outside the loop into forecasts and standard errors.

compute %modelsetcoeffs(varmodel,bdraw)

simulate(model=varmodel,cv=sigmad,results=simresults,steps=nstep)

do i=1,nvar

set forecast(i) fstart fend = forecast(i)+simresults(i)

set forestderr(i) fstart fend = forestderr(i)+simresults(i)^2

end do i

That finishes the draw loop. The processing outside the loop requires converted the raw sums and sums of squares into means and standard deviations, and graphing the results. The graphs include the forecast (mean of the simulated out-of-sample data) , upper and lower 2 standard deviation limits and the last 12 period of actual data.

do i=1,nvar

set forecast(i) fstart fend = forecast(i)/ndraws

set forestderr(i) fstart fend = sqrt(forestderr(i)/ndraws-forecast(i)^2)

end do i

do i=1,nvar

set lower fstart fend = forecast(i)-2.0*forestderr(i)

set upper fstart fend = forecast(i)+2.0*forestderr(i)

graph(header="Forecasts of "+%modellabel(varmodel,i)) 4

# forecast(i) fstart fend

# lower fstart fend 2

# upper fstart fend 2

# %modeldepvars(varmodel)(i) fstart-12 fstart-1

end do i

Full Program

compute lags=4          ;*Number of lags
compute nstep=16        ;*Number of response steps
compute nburn=500       ;*Number of burn-in draws
compute ndraws=2500     ;*Number of keeper draws
*
open data haversample.rat
cal(q) 1959
data(format=rats) 1959:1 2006:4 ftb3 gdph ih cbhm
*
set loggdp = log(gdph)
set loginv = log(ih)
set logc   = log(cbhm)
*
* These are the controlling parameters for the symmetric "Minnesota"
* prior - the own lag tightness and the relative tightness on the other
* lags.
*
compute tight=.1
compute other=.5
*
* Set up the VAR
*
system(model=varmodel)
variables loggdp loginv logc ftb3
lags 1 to lags
det constant
end(system)
*
* Set up the prior
*
@BVARBuildPriorMN(model=varmodel,tight=tight,other=other) hbpriors hpriors
*
* Convert the working arrays to the proper form
*
@BVARFinishPrior hbpriors hpriors bprior hprior
*
estimate(ols,sigma)
compute nvar=%nvar
compute sigmad=%sigma
cmom(model=varmodel)
*
dec vect[series] forecast(nvar)
dec vect[series] forestderr(nvar)
compute fstart=%regend()+1
compute fend =fstart+nstep-1
*
do i=1,nvar
   set forecast(i)   fstart fend = 0.0
   set forestderr(i) fstart fend = 0.0
end do i
*
infobox(action=define,progress,lower=-nburn,upper=ndraws) "Gibbs Sampler"
do draw=-nburn,ndraws
   infobox(current=draw)
   *
   * Draw b given sigma
   *
   compute bdraw =%ranmvkroncmom(%cmom,inv(sigmad),hprior,bprior)
   compute rssmat=%sigmacmom(%cmom,bdraw)
   *
   * Draw sigma given b
   *
   compute sigmad=%ranwisharti(%decomp(inv(rssmat)),%nobs)
   if draw<=0
      next

   compute %modelsetcoeffs(varmodel,bdraw)
   simulate(model=varmodel,cv=sigmad,results=simresults,steps=nstep)
   do i=1,nvar
      set forecast(i)   fstart fend = forecast(i)+simresults(i)
      set forestderr(i) fstart fend = forestderr(i)+simresults(i)^2
   end do i
end do draw
infobox(action=remove)
*
do i=1,nvar
   set forecast(i) fstart fend   = forecast(i)/ndraws
   set forestderr(i) fstart fend = sqrt(forestderr(i)/ndraws-forecast(i)^2)
end do i
*
do i=1,nvar
   set lower fstart fend = forecast(i)-2.0*forestderr(i)
   set upper fstart fend = forecast(i)+2.0*forestderr(i)
   graph(header="Forecasts of "+%modellabel(varmodel,i)) 4
   # forecast(i) fstart fend
   # lower fstart fend 2
   # upper fstart fend 2
   # %modeldepvars(varmodel)(i) fstart-12 fstart-1
end do i

Graph (One of Four)