*
* This demonstrates a technique for drawing from a VAR when the data
* have changed so new data are used in each Gibbs sweep. In this case,
* the data stay the same, so it could be done (much) more efficiently by
* using the sample statistics from ESTIMATE, as is done in the standard
* VAR Monte Carlos.
*
compute lags=4			;*Number of lags
compute nstep=16		;*Number of response steps
compute ndraw=10000	;*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)
*
system(model=varmodel)
variables loggdp loginv logc ftb3
lags 1 to lags
det constant
end(system)
******************************************************************
estimate
compute nvar   =%nvar
compute wishdof=%nobs-%nreg
*
dec vect[rect] %%responses(ndraw)
*
infobox(action=define,progress,lower=1,upper=ndraw) "Monte Carlo Integration"
do draws=1,ndraw
   *
   * This computes the cross product matrix in the order X~Y.
   *
   cmom(model=varmodel)
   *
   * This is for a sampler with the standard Jeffrey's prior, so sigma
   * can be drawn unconditionally and beta can be drawn conditional on
   * sigma. sweepxx gives the inv(X'X) matrix in the top block, the
   * coefficients in the upper right and T x sigma in the bottom right
   * corner.
   *
   compute sweepxx =%sweeptop(%cmom,%nreg)
   compute fwish   =%decomp(inv(%xsubmat(sweepxx,%nreg+1,%ncmom,%nreg+1,%ncmom)))
   compute sigmad  =%ranwisharti(fwish,wishdof)
   compute fsigma  =%decomp(sigmad)
   compute fxx     =%decomp(%xsubmat(sweepxx,1,%nreg,1,%nreg))
   compute betaols =%xsubmat(sweepxx,1,%nreg,%nreg+1,%ncmom)
   compute betadraw=betaols+%ranmvkron(fsigma,fxx)
   compute %modelsetcoeffs(varmodel,betadraw)
   *
   impulse(noprint,model=varmodel,factor=fsigma,results=impulses,steps=nstep)
   *
   * Store the impulse responses
	*
   dim %%responses(draws)(nvar*nvar,nstep)
   ewise %%responses(draws)(i,j)=ix=%vec(%xt(impulses,j)),ix(i)
   infobox(current=draws)
end do draws
infobox(action=remove)
*
@mcgraphirf(model=varmodel)