GIBBSVAR.RPF

GIBBSVAR.RPF is an example of Gibbs sampling applied to a Bayesian VAR with a "Minnesota" prior. This has been superseded by GIBBSVARBUILD.RPF which is effectively identical (except for a minor difference in choosing the scale factors for the variables) but uses a special purpose procedure to generate the matrices which describe the Minnesota prior.

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
*
* First, estimate the system with a prior via mixed estimation, to see if
* we come up with something similar via Gibbs sampling.
*
system(model=varmodel)
variables loggdp loginv logc ftb3
lags 1 to lags
specify(type=symmetric,tight=tight) other
det constant
end(system)
************************************************************************************
*
* Estimate with a prior
*
estimate
*
* Same system without a prior. This is needed to get the likelihood
* information.
*
estimate(ols,sigma)
compute nvar=%nvar
*
* Get the estimated standard errors from the regressions for use in
* scaling the prior.
*
dec vect olssee(nvar)
ewise olssee(i)=sqrt(%sigma(i,i))
*
* Create the prior mean and (square root of) precision. The olssee
* values are used to scale the precision. The prior mean is 1 for the
* first own lag in each equation, and 0 otherwise. The constant is given
* a mean of zero with a zero precision (that is, a flat prior).
*
compute ncoef=lags*nvar+1
dec rect minnmean(ncoef,nvar) minnprec(ncoef,nvar)
*
do i=1,nvar
   do j=1,nvar
      do l=1,lags
         compute minnmean((j-1)*lags+l,i)=(i==j.and.l==1)
         compute minnprec((j-1)*lags+l,i)=olssee(j)/olssee(i)*%if(i==j,1.0/tight,1.0/(other*tight))
      end do l
   end do j
   compute minnmean(ncoef,i)=0.0,minnprec(ncoef,i)=0.0
end do i
*
ewise minnprec(i,j)=minnprec(i,j)^2
compute hprior=%diag(%vec(minnprec))
compute bprior=minnmean
*
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