*
* 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. Since these are all off-diagonal, we can try just zeros
* for the guess values.
*
nonlin(parmset=simszha,zeros) 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 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) $
"Metropolis within Gibbs"
do draw=-nburn,ndraws
*
* Unlike the full VAR in SIMSIRF.RPF, 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.RPF, 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 %ranflip(alpha)
compute thetadraw=theta,logplast=logptest,accept=accept+1
infobox(current=draw) %strval(100.0*accept/(draw+nburn+1),"##.#")
*
* Conditioned on theta, make a draw for lambda.
*
compute %parmspoke(simszha,thetadraw)
compute a=afrml(1)
compute vdiag =%mqformdiag(sigmaAtBeta,tr(a))
ewise lambdadraw(i)=(tnobs/2.0)*vdiag(i)/%rangamma(.5*(tnobs)+delta+1)
*
* Combine Lambda and A to generate the factor for the SVAR (for
* use in computing the impulse responses), and the draw for the
* structural precision (inverse covariance) matrix.
*
compute factor=inv(a)*%diag(%sqrt(lambdadraw))
compute hdraw=tr(a)*inv(%diag(lambdadraw))*a
*
* 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 bdraw=%ranmvposthb(hdata,hbdata)
*
* Compute covariance matrix of residuals at the current beta.
* (This is used only for the next draw for the SVAR parameters).
*
compute sigmaAtBeta=SURGibbsSigma(bdraw)/tnobs
if draw<=0
next
*
* Do the bookkeeping here.
*
compute tgibbs(draw)=thetadraw
compute %modelsetcoeffs(nearvar,bdraw)
impulse(noprint,model=nearvar,factor=factor,steps=nsteps,$
flatten=%%responses(draw))
*
* Store the impulse responses
*
end do draw
infobox(action=remove)
*
@mcgraphirf(model=nearvar,$
shocklabels=shocklabels,varlabels=varlabels,$
percent=||.025,.16,.84,.975||,center=median,$
footer="Pointwise 68% and 95% Posterior Bands, Six Variable Near-VAR")