*
* MONTENEARSVAR.RPF
* Example of MCMC analysis of a combination of a near VAR for the lag
* coefficients and a structural VAR for the covariance matrix. This is
* based upon the SIMSIRF.RPF example from the Sims and Zha Econometrica
* 1999 replication files.
*
open data simszha.xls
calendar(q) 1948
data(format=xls,org=columns) 1948:01 1989:03 gin82 gnp82 gnp gd lhur fygn3 m1
*
set gin82 = log(gin82)
set gnp82 = log(gnp82)
set gd = log(gd)
set m1 = log(m1)
set fygn3 = fygn3*.01
set lhur = lhur*.01
*
compute nlags =4
compute nvar =6
compute nsteps=25
compute nburn =2000
compute ndraws=10000
*
* Create a near-VAR that treats gin82 as exogenous.
*
system(model=fivevar)
variables fygn3 m1 gnp82 gd lhur
lags 1 to nlags
det constant gin82{1 to nlags}
end(system)
*
equation gineqn gin82
# constant gin82{1 to nlags}
*
compute nearvar=fivevar+gineqn
*
* Set up the Gibbs sampler for the SUR system
*
@SURGibbsSetup nearvar
*
* Do a SUR to get a set of estimates to initialize the Gibbs sampler.
*
sur(model=nearvar,cvout=sigmaAtBeta)
compute tnobs =%nobs
*
* 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||
*
* 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. This will be used to initialize the M-H sampler and
* provide (with scaling) the covariance matrix for the increment for
* Random Walk Metropolis.
*
compute delta=3.5
cvmodel(parmset=simszha,pdf=delta,dmatrix=marginalized,method=bfgs) sigmaAtBeta afrml
*
compute shocklabels=||"MS","MD","Y","P","U","I"||
compute varlabels=||"R","M","Y","P","U","I"||
*
* This is the covariance matrix for the RW increment. We might have to
* change the scaling on this to achieve a reasonable acceptance
* probability.
*
compute f=%decomp(%xx)*.35
compute nudraw=6
*
* We need to save the draws in %%RESPONSES as described in the "RATS
* User's Guide" (for use with MCGRAPHIRF).
*
declare rect[series] impulses(nvar,nvar)
declare vect[rect] %%responses(ndraws)
*
compute thetaDraw =%parmspeek(simszha)
compute accept =0
declare vect lambdadraw(nvar) vdiag(nvar)
*
* This is for saving the draws for the structural VAR parameters.
*
dec series[vect] tgibbs
gset tgibbs 1 ndraws = thetaDraw
*
infobox(action=define,progress,lower=-nburn,upper=ndraws) $
"Gibbs Sampling"
do draw=-nburn,ndraws
*
* Unlike the full VAR in SIMSIRF, the likelihood for a near VAR
* doesn't isolate the structural model for the covariance matrix from
* the lag coefficients. For the full VAR, the structural coefficients
* can be drawn using only the covariance matrix from the OLS
* estimates; here, we need to use the recomputed covariance matrix at
* each Gibbs sweep. Also unlike SIMSIRF, we need to draw coefficients
* with each sweep, not just after the burn-in. (With the full VAR,
* the coefficients are needed only for doing the impulse response
* functions).
*
* So the procedure is:
*
* 1. Compute the log likelihood for the structural model given the
* covariance matrix of residuals from the current draw for the
* coefficients (<>).
*
* 2. Draw a candidate set of structural parameters and compute the
* log likelihood at those.
*
* 3. Do a Metropolis acceptance test using those two to determine
* whether to reject or accept the candidate draw.
*
* 4. Draw the diagonal elements for the structural covariance matrix
* using the chosen set of structural parameters and <>.
*
* 5. Draw the coefficients from the SUR, and recompute <>
* for use in the next sweep.
*
* Note that there is no DFC option in the CVMODEL---that's needed only
* when the beta's can be marginalized out.
*
cvmodel(parmset=simszha,pdf=delta,dmatrix=marginalized,method=evaluate,a=afrml) sigmaAtBeta
compute logplast=%funcval
*
* Draw a new theta based at previous value
*
compute theta=thetadraw+%ranmvt(f,nudraw)
*
* Evaluate the model there
*
compute %parmspoke(simszha,theta)
cvmodel(parmset=simszha,pdf=delta,dmatrix=marginalized,method=evaluate,a=afrml) sigmaAtBeta
compute logptest=%funcval
*
* Compute the acceptance probability
*
compute alpha =exp(logptest-logplast)
if alpha>1.0.or.%uniform(0.0,1.0)