*
* MONTESVAR_MH.RPF
* Monte Carlo integration of a structural VAR using Random Walk
* Metropolis.
*
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)