*
* MONTESVAR.RPF
* RATS Version 8, User's Guide, Described in Section 7.5.
*
open data haversample.rat
calendar(q) 1959
data(format=rats) 1959:1 2006:4 ffed fm1 fnh gdph lr dgdp
*
compute lags=4
compute nvar=6
compute nstep=16
compute ndraws=10000
*
set ffed = ffed*.01
set lr   = lr*.01
log gdph
log dgdp
log fnh
log fm1
*
* In this section, we set up a system for a structural model. First,
* estimate the VAR by OLS.
*
system(model=varmodel)
variables ffed fm1 gdph dgdp lr fnh
lags 1 to lags
det constant
end(system)
*
estimate(cvout=vmat)
compute xxx=%xx
*
* Save a factor of the inv(X'X) matrix from the VAR, and the OLS
* coefficients.
*
compute fxx    =%decomp(xxx)
compute betaols=%modelgetcoeffs(varmodel)
compute ncoef  =%rows(fxx)
*
* Create the set of non-linear parameters, and the "A" formula for
* CVMODEL.
*
nonlin(parmset=simszha) a12 a21 a23 a24 a31 a36 a41 a43 a46 a51 a53 a54 a56
compute nfree=13
*
dec frml[rect] afrml
frml afrml = ||1.0   ,a12  ,0.0   ,0.0   ,0.0 ,0.0|$
               a21   ,1.0  ,a23   ,a24   ,0.0 ,0.0|$
               a31   ,0.0  ,1.0   ,0.0   ,0.0 ,a36|$
               a41   ,0.0  ,a43   ,1.0   ,0.0 ,a46|$
               a51   ,0.0  ,a53   ,a54   ,1.0 ,a56|$
               0.0   ,0.0  ,0.0   ,0.0   ,0.0 ,1.0||
*
declare rect fsigmad betaols betadraw
declare symm sigmad
*
* Compute the maximum of the log of the marginal posterior density for
* the A coefficients with a prior of the form |D|**(-delta). Delta
* should be at least (nvar+1)/2.0 to ensure an integrable posterior
* distribution.
*
compute delta=3.5
cvmodel(parmset=simszha,dfc=ncoef,pdf=delta,dmatrix=marginalized,method=bfgs) vmat afrml
dec rect faxx
*
* axbase is the maximizing vector of coefficients.
*
compute [vector] axbase=%parmspeek(simszha)
*
* <<faxx>> is a factor of the (estimated) inverse Hessian at the final
* estimates. This gets scaled slightly to fatten up the tails a bit.
*
compute faxx=1.2*%decomp(%xx)
*
* <<scladjust>> is used to prevent overflows when computing the weight
* function.
*
compute scladjust=%funcval
*
* <<nu>> is the degrees of freedom for the multivariate Student used in
* drawing A's.
*
compute nu=30.0
*
* We need to save the draws in %%RESPONSES as described in the manual
* (for use with MCGRAPHIRF). We also need a series to hold the weights.
*
set weights 1 ndraws = 0.0
declare vect au(%nreg)
declare vect dbase d(nvar)
declare rect[series] impulses(nvar,nvar)
*
* This is the standard name used by the MCGraphIRF procedure
*
declare vect[rect] %%responses(ndraws)
*
* sumwt and sumwt2 will hold the sum of the weights and the sum of
* squared weights.
*
compute sumwt=0.0
compute sumwt2=0.0
infobox(action=define,progress,lower=1,upper=ndraws) "Monte Carlo Integration"
*
* Antithetic acceleration is used for the lag coefficients only, so on
* odd number draws, we generate a new sigma matrix, and a set of changes
* to the lag coefficients, then on the even ones, we generate a new draw
* with the sign switched on the deltas to the lag coefficients.
*
do draw=1,ndraws
   if %clock(draw,2)==1 {
      *
      * Do a draw for the coefficients from a multivariate t density and
      * "poke" it back into the parmset so the AFRML can get the new
      * values. Save the log kernel of the draw into <<idensity>>.
      *
      compute au      =%rant(nu)
      compute idensity=%ranlogkernel()
      compute %parmspoke(simszha,axbase+faxx*au)
      *
      * Evaluate the model at the draw for the coefficients. Save the
      * log likelihood.
      *
      compute a       =afrml(1)
      compute dhat    =a*vmat*tr(a)
      compute ddiag   =%xdiag(dhat)
      cvmodel(parmset=simszha,dfc=ncoef,pdf=delta,dmatrix=marginalized,method=evaluate) vmat afrml
      compute pdensity=%funcval
      *
      * Compute the weight value by exp'ing the difference between the
      * two densities, with a scale adjustment term to prevent overflow.
      *
      compute weight  =exp(pdensity-scladjust-idensity)
      *
      * Conditioned on A, make a draw for the D matrix
      *
      ewise d(i)      =(%nobs/2.0)*ddiag(i)/%rangamma(.5*(%nobs-ncoef)+delta+1)
      *
      * Combine D and A to generate the draw for a factor of sigma.
      *
      compute fsigmad =inv(a)*%diag(%sqrt(d))
		*
		* Compute the + draw for coefficients
		*
      compute betau   =%ranmvkron(fsigmad,fxx)
      compute betadraw=betaols+betau
   }
   else
		*
      * Compute the - draw for the coefficients
		*
      compute betadraw=betaols-betau

   compute %modelsetcoeffs(varmodel,betadraw)
   compute nshock=1
   impulse(noprint,model=varmodel,factor=fsigmad,results=impulses,steps=nstep) nvar
	*
	* Store the impulse responses and the draw weight
	*
   compute sumwt   =sumwt+weight
   compute sumwt2  =sumwt2+weight^2
   dim %%responses(draw)(nvar*3,nstep)
   compute weights(draw)=weight
   ewise %%responses(draw)(i,j)=xt=%xt(impulses,j),ix=%vec(%xcol(xt,3)~%xcol(xt,5)~%xcol(xt,6)),ix(i)
   infobox(current=draw)
end do draw
infobox(action=remove)
*
* The efficacy of importance sampling depends upon function being
* estimated, but the following is a simple estimate of the number of
* effective draws.
*
disp "Effective sample size" sumwt^2/sumwt2
*
* Normalize the weights to sum to 1.0
*
set weights 1 ndraws = weights/sumwt
*
@MCGraphIRF(model=varmodel,center=mean,page=all,weights=weights,$
  shocklabels=||"P","U","I"||)