* * 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)