*
* MONTESUR.RPF
* Uses Gibbs sampling to analyze the impulse response functions for a
* near-VAR.
*
compute lags=4			   ;*Number of lags
compute steps=16		   ;*Number of response steps
compute nburn =5000		;*Burn-in draws
compute ndraws=25000		;*Keeper draws
*
open data haversample.rat
calendar(q) 1959
data(format=rats) 1959:1 2006:4 fm1 gdph
*
log gdph
log fm1
*
equation gdpeq gdph
# gdph{1 to lags} fm1{1 to lags} constant
equation m1eq fm1
# fm1{1 to lags} constant
group surmodel m1eq gdpeq
compute nvar=%modelsize(surmodel)
*
* Do set up for Gibbs sampling on SUR model
*
@SURGibbsSetup surmodel
**********************************************
*
* Do a SUR to get a set of estimates to initialize the Gibbs sampler.
*
sur(model=surmodel)
compute ntotal=%nreg
compute bdraw=%beta
compute wishdof=%nobs
*
* For saving the IRF's
*
declare vect[rect] %%responses(ndraws)
*
infobox(action=define,progress,lower=-nburn,upper=ndraws) $
    "Gibbs Sampling"
do draw=-nburn,ndraws
   *
   * Compute covariance matrix of residuals at the current beta
   *
   compute covmat=SURGibbsSigma(bdraw)
   *
   * Do a draw for the precision matrix conditional on beta
   *
   compute hdraw=%ranwishartf(%decomp(inv(covmat)),wishdof)
   *
   * Compute the required information with the interaction between
   * the precision matrix and the data
   *
   @SURGibbsDataInfo hdraw hdata hbdata
   *
   * Draw coefficients given the precision matrix hdraw
   *
   compute hpost=hdata
   compute vpost=inv(hpost)
   compute bpost=vpost*hbdata
   compute bdraw=bpost+%ranmvnormal(%decomp(vpost))
   infobox(current=draw)
   if draw<=0
      next
   *
   * Do the bookkeeping here.
   *
   compute %modelsetcoeffs(surmodel,bdraw)
   impulse(noprint,model=surmodel,factor=%decomp(hdraw),results=impulses,steps=steps)
   *
   * Store the impulse responses
   *
   dim %%responses(draw)(nvar*nvar,steps)
   ewise %%responses(draw)(i,j)=ix=%vec(%xt(impulses,j)),ix(i)
end do draw
infobox(action=remove)
@MCGraphIRF(model=surmodel,center=median,page=all)