MONTESVAR_MH.RPF is an example of Monte Carlo integration of a structural VAR (SVAR) using Random Walk Metropolis. MONTESVAR.RPF is similar, but uses importance sampling. Importance sampling is generally superior if it can be made to work, however, that typically requires a structural model where the parameters are well-estimated by CVMODEL, which usually means the model is a restriction of a Cholesky factor, or something close to it. Random Walk Metropolis likely requires more experimentation to come up with a good method for generating candidates, but can much more easily handle likelihood surfaces with odd shapes.

This uses the six-variable SVAR used by Sims in a series of papers. The six variables are "R" (interest rate), "M" (log money), "Y" (log real GDP"), "P" (log GDP deflator), "U" (unemployment) and "I" (log real investment) (in order in the system below).

system(model=varmodel)

variables ffed fm1 gdph dgdp lr fnh

lags 1 to lags

det constant

end(system)

*

estimate(cvout=vmat)

This takes information from the estimation that will be needed for drawing lag coefficients later (the inverse X'X matrix and the coefficients).

compute xxx=%xx

compute fxx    =%decomp(xxx)

compute betaols=%modelgetcoeffs(varmodel)

compute ncoef  =%rows(fxx)

compute nvar   =%nvar

This sets up the structural model (A form). The six shocks are labeled as Money Supply (MS), Money Demand (MD), Output, Price, Unemployment and Investment. The biggest statistical problem with the model turns out to be separating MS and MD. While the point estimates show the expected signs for A12 and A21, they're both statistically insignificant, so it's possible that the signs could flip.

nonlin(parmset=simszha) a12 a21 a23 a24 a31 a36 $a41 a43 a46 a51 a53 a54 a56 * 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(a=afrml,parmset=simszha,dfc=ncoef,pdf=delta,$dmatrix=marginalized,method=bfgs) vmat compute nfree=%nreg This starts the chain in the MCMC procedure at the maximizer for the log posterior. LOGPLAST is the value of the log posterior at that point. compute [vector] thetadraw=%parmspeek(simszha) compute logplast=%funcval This is the (factor of the) covariance matrix for the RW increment. We might have to change the scaling on this to achieve a reasonable acceptance probability. What's chosen here is a diagonal covariance matrix (thus the components are chosen independently) using the estimated variances from the CVMODEL estimation. The .5 scale factor means that we're taking .5 of an estimated standard deviation on each parameter. compute f=0.5*%diag(%sqrt(%xdiag(%xx))) This sets the number of burn-in and saved draws and initializes the number of acceptances. compute nburn =5000 compute ndraws=10000 compute accept=0 We're saving the draws for the free parameters, and the impulse responses. The latter uses %%RESPONSES so @MCGRAPHIRF (or @MCPROCESIRF) can be used to process the draws for the IRF's. declare vect lambdadraw(nvar) vdiag(nvar) dec series[vect] tgibbs gset tgibbs 1 ndraws = %zeros(nfree,1) * declare vect[rect] %%responses(ndraws) This uses the standard structure for doing Metropolis-Hastings MCMC, with an INFOBOX to let the user know about the progress, non-positive numbers of DRAW used for the burn-in draws and positive values for the saved draws. infobox(action=define,progress,lower=-nburn,upper=ndraws)$

"Random Walk MH"

*

do draw=-nburn,ndraws

...do draws

infobox(current=draw) %strval(100.0*accept/(draw+nburn+1),"##.#")

if draw<=0

next

...save information

end do draw

infobox(action=remove)

Inside the loop, this takes a new draw, evaluates the posterior there, and does the Metropolis test for whether to retain the previous draw or accept the candidate.

*

* Draw a new theta based at previous value

*

compute theta=thetadraw+%ranmvnormal(f)

*

* Evaluate the model there

*

compute %parmspoke(simszha,theta)

cvmodel(parmset=simszha,dfc=ncoef,pdf=delta,$dmatrix=marginalized,method=evaluate,a=afrml) vmat compute logptest=%funcval compute alpha =exp(logptest-logplast) if %ranflip(alpha) compute thetadraw=theta,logplast=logptest,accept=accept+1 We only have to draw the VAR coefficients and compute the impulse responses when we are saving draws. We need to draw the "lambda" (diagonal matrix of the variance of the shocks) conditional on the "theta" (structural coefficients). The two combined give us FSIGMAD, which is a draw for the factor of the covariance matrix. This is the proper calculation for an "A" type SVAR. * * Save the current draw for theta * compute tgibbs(draw)=thetadraw * * Conditioned on theta, make a draw for lambda. (We don't need to do * this until we're saving results, since theta is being drawn * unconditionally). * compute %parmspoke(simszha,thetadraw) compute a=afrml(1) compute vdiag =%mqformdiag(vmat,a) ewise lambdadraw(i)=$

(%nobs/2.0)*vdiag(i)/%rangamma(.5*(%nobs-ncoef)+delta+1)

compute fsigmad =inv(a)*%diag(%sqrt(lambdadraw))

Given FSIGMAD, this draws the lag coefficients for the VAR by standard method. The new coefficients are "poked" into the model. This program is only keeping tracking of three of the shocks (P, U and I, which are the 3rd, 5th and 6th as the model is constructed), so the "factor" is reduced to just a 6 × 3 array which those columns. IMPULSE computes the responses to that set of shocks and "flattens" them into the form needed by %%RESPONSES.

*

* Compute the draw for the lag coefficients

*

compute betadraw=betaols+%ranmvkron(fsigmad,fxx)

compute %modelsetcoeffs(varmodel,betadraw)

*

* Extract the P, U and I shocks

*

compute shockpui=%xcol(fsigmad,3)~%xcol(fsigmad,5)~%xcol(fsigmad,6)

impulse(noprint,model=varmodel,factor=shockpui,$flatten=%%responses(draw),steps=nstep) nvar Outside the loop, @MCGraphIRF is used to generate graphs of the responses with error bands: @MCGraphIRF(model=varmodel,center=mean,page=all,shocklabels=||"P","U","I"||) and @MCMCPostProc is used to analyze the draws for the structural VAR coefficients: @mcmcpostproc(ndraws=ndraws,mean=bmean,$

stderrs=bstderrs,cd=bcd,nse=bnse) tgibbs

report(action=define)

report(atrow=1,atcol=1,align=center) "Variable" "Coeff" $"Std Error" "NSE" "CD" do i=1,%nreg report(row=new,atcol=1) %parmslabels(simszha)(i) bmean(i)$

bstderrs(i) bnse(i) bcd(i)

end do i

report(action=format,atcol=2,tocol=3,picture="*.###")

report(action=format,atcol=4,picture="*.###")

report(action=show)

Full Program

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
*
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)
compute nvar   =%nvar
*
* Create the set of non-linear parameters, and the "A" formula for
* CVMODEL. Since these are all off-diagonal coefficients, all zeros
* has a good chance of working as guess values.
*
nonlin(parmset=simszha,zeros) a12 a21 a23 a24 a31 a36 $a41 a43 a46 a51 a53 a54 a56 * 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(a=afrml,parmset=simszha,dfc=ncoef,pdf=delta,$dmatrix=marginalized,method=bfgs) vmat compute nfree=%nreg * * Start the chain at the maximizer * compute [vector] thetadraw=%parmspeek(simszha) compute logplast=%funcval * * 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=0.5*%diag(%sqrt(%xdiag(%xx))) * compute nburn =5000 compute ndraws=10000 compute accept=0 declare vect lambdadraw(nvar) vdiag(nvar) dec series[vect] tgibbs gset tgibbs 1 ndraws = %zeros(nfree,1) * * We need to save the draws in %%RESPONSES as described in the manual * (for use with @MCGRAPHIRF). * declare vect[rect] %%responses(ndraws) * infobox(action=define,progress,lower=-nburn,upper=ndraws)$
"Random Walk MH"
*
do draw=-nburn,ndraws
*
* Draw a new theta based at previous value
*
compute theta=thetadraw+%ranmvnormal(f)
*
* Evaluate the model there
*
compute %parmspoke(simszha,theta)
cvmodel(parmset=simszha,dfc=ncoef,pdf=delta,$dmatrix=marginalized,method=evaluate,a=afrml) vmat compute logptest=%funcval 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),"##.#") if draw<=0 next * * Save the current draw for theta * compute tgibbs(draw)=thetadraw * * Conditioned on theta, make a draw for lambda. (We don't need to do * this until we're saving results, since theta is being drawn * unconditionally). * compute %parmspoke(simszha,thetadraw) compute a=afrml(1) compute vdiag =%mqformdiag(vmat,a) ewise lambdadraw(i)=$
(%nobs/2.0)*vdiag(i)/%rangamma(.5*(%nobs-ncoef)+delta+1)

compute fsigmad =inv(a)*%diag(%sqrt(lambdadraw))
*
* Compute the draw for the lag coefficients
*
compute betadraw=betaols+%ranmvkron(fsigmad,fxx)
compute %modelsetcoeffs(varmodel,betadraw)
*
* Extract the P, U and I shocks
*
compute shockpui=%xcol(fsigmad,3)~%xcol(fsigmad,5)~%xcol(fsigmad,6)
impulse(noprint,model=varmodel,factor=shockpui,$flatten=%%responses(draw),steps=nstep) nvar end do draw infobox(action=remove) * @MCGraphIRF(model=varmodel,center=mean,page=all,shocklabels=||"P","U","I"||) * * Monte Carlo statistics/diagnostics on the SVAR coefficients * @mcmcpostproc(ndraws=ndraws,mean=bmean,$
stderrs=bstderrs,cd=bcd,nse=bnse) tgibbs
report(action=define)
report(atrow=1,atcol=1,align=center) "Variable" "Coeff" $"Std Error" "NSE" "CD" do i=1,%nreg report(row=new,atcol=1) %parmslabels(simszha)(i) bmean(i)$
bstderrs(i) bnse(i) bcd(i)
end do i
report(action=format,atcol=2,tocol=3,picture="*.###")
report(action=format,atcol=4,picture="*.###")
report(action=show)


Output

(The VAR output isn't shown here.)

This is the output from the maximum posterior estimates of the model.

Covariance Model-Marginal Posterior - Estimation by BFGS

Convergence in    32 Iterations. Final criterion was  0.0000043 <=  0.0000100

Observations                              188

Log Posterior                       3867.5518

Prior DF                               3.5000

Variable                        Coeff      Std Error      T-Stat      Signif

************************************************************************************

1.  A12                          -0.210169018  0.328503265     -0.63978  0.52231722

2.  A21                           0.239494036  0.266017749      0.90029  0.36796411

3.  A23                          -0.130823823  0.104316657     -1.25410  0.20980460

4.  A24                          -0.376127996  0.222527257     -1.69026  0.09097902

5.  A31                          -0.007454838  0.068674529     -0.10855  0.91355690

6.  A36                          -0.232790224  0.028240889     -8.24302  0.00000000

7.  A41                          -0.046540455  0.030370724     -1.53241  0.12542085

8.  A43                           0.057057884  0.034831917      1.63809  0.10140251

9.  A46                          -0.014429582  0.015065034     -0.95782  0.33815382

10. A51                           0.098851055  0.017241422      5.73335  0.00000001

11. A53                           0.156988068  0.021989863      7.13911  0.00000000

12. A54                           0.035038164  0.049088412      0.71378  0.47536528

13. A56                           0.007954866  0.009414004      0.84500  0.39810888

These are the corresponding results from the Markov Chain Monte Carlo process. (Note that these are subject to simulation error). Note that the coefficients on the first two structural equations are considerably different from the point estimates, and have poor "CD" diagnostics. This particular model is known to have difficulties separating "money supply" and "money demand."

Variable Coeff  Std Error  NSE    CD

A12      -0.821     0.631 0.054  4.298

A21       0.692     0.440 0.037 -4.175

A23      -0.304     0.221 0.019  3.754

A24      -0.645     0.545 0.046  4.808

A31       0.059     0.095 0.007 -1.202

A36      -0.242     0.028 0.002  1.329

A41      -0.020     0.048 0.004 -2.675

A43       0.054     0.035 0.002  1.550

A46      -0.018     0.015 0.001 -0.065

A51       0.099     0.018 0.001 -0.449

A53       0.154     0.022 0.001 -4.069

A54       0.029     0.057 0.004 -1.680

A56       0.009     0.009 0.001  3.114

Graphs