|
Examples / MONTESVAR_MH.RPF |
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 @MCPROCESSIRF) 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 \times 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
Copyright © 2024 Thomas A. Doan