WZSAMPLER.RPF

WZSAMPLER.RPF uses the Waggoner-Zha(2003) Gibbs sampler to analyze the A form model from the CVMODEL.RPF example. (The WZ sampler only applies to A form models).

Their sampler draws each row of A separately. The advantage that this has over importance sampling or Metropolis-Hastings is that it doesn't require experimentation with "tuning" parameters---it's possible to draw the free coefficients in each row directly (conditional on the others).

This does not assume a normalization. This has both good and bad points. The good point is that it can't create a problem by choosing a normalization on a coefficient whose sign isn't determined. The bad point is that interpreting a structural VAR without a normalization can be difficult---does it make sense for a "price'' shock to be of undetermined sign or for a "money demand'' shock to not hit money? The sampler by construction picks a random sign fairly deep in the calculation, so this produces rows with roughly 50-50 signs. How this gets converted to a usable normalization is not clear. Based upon the examples that I've done, it looks like picking a normalization on a row that has the smallest (most negative) kurtosis tends to work well---that's the parameter that is most "bimodal".

This is a rather straightforward example of a structural VAR model without any real identification issues, so the WZ sampler gives basically the same results (once normalized) as Metropolis-Hastings draws.

A large block of the program (between the full lines of ***) is a generic implementation of the sampler and doesn't need to change. As it says in the comments after that,

* WZBDRAWS now has the sign-corrected draws for all the elements. What you do with

* this will obviously depend upon the application.

This does a scatter graph of the MR and MM coefficients (from the money equation) which are the 10th and 11th free parameters in the model. (The SMPL option limits this to only 1 in 10 of the points, as drawing 25000 dots takes a while).

set mrx 1 ndraws = wzbdraws(t)(10)

set mmx 1 ndraws = wzbdraws(t)(11)

scatter(footer="Coefficients from Money Demand Equation",$

hlabel="M coefficient",vlabel="R coefficient",smpl=%clock(t,10)==10)

# mmx mrx

The equation gets sign-normalized on the M coefficient, so those are all positive. Note that the money equation (as is the case with all the equations in a typical SVAR "A" model) has all coefficients on the left side, rather than M on the left and R on the right, so if this represents money demand, we would expect the R coefficients to be positive as well. That's is largely (but not exclusively) true.

This does a density graph of the normalized (now to a unit coefficient on M, not just positive sign) of the R coefficient in the M equation.

set mrstd 1 ndraws = mrx/mmx

density(smoothing=1.5,grid=automatic,maxgrid=100) mrstd / mrv mrf

scatter(style=line,vmin=0.0,footer="Posterior Density of Normalized MR coefficient")

# mrv mrf

The final part of the program takes the draws for the SVAR model as given and generates error bands by drawing lag coefficients and doing impulse responses on the combination of the SVAR and VAR and does a standard MC graph matrix for impulse responses.

Full Program

compute nburn =10000

compute ndraws=25000

compute nsteps=24

open data oecdsample.rat

calendar(q) 1981

data(format=rats) 1981:1 2006:4 can3mthpcp canexpgdpchs $

canexpgdpds canm1s canusxsr usaexpgdpch

set logcangdp = 100.0*log(canexpgdpchs)

set logcandefl = 100.0*log(canexpgdpds)

set logcanm1 = 100.0*log(canm1s)

set logusagdp = 100.0*log(usaexpgdpch)

set logexrate = 100.0*log(canusxsr)

system(model=canmodel)

variables logusagdp logcangdp can3mthpcp logexrate logcanm1 logcandefl

lags 1 to 4

det constant

end(system)

estimate(noprint)

* Set up A matrix with no sign normalizations

dec frml[rect] afrml

nonlin(parmset=aparms) uu ur cu cc cr rr rx xx mc mr mm mp pc pr pp

frml afrml = ||uu ,0.0,ur ,0.0,0.0,0.0|$

cu ,cc ,cr ,0.0,0.0,0.0|$

0.0,0.0,rr ,rx ,0.0,0.0|$

0.0,0.0,0.0,xx ,0.0,0.0|$

0.0,mc ,mr ,0.0,mm ,mp |$

0.0,pc ,pr ,0.0,0.0,pp||

* Estimate the model by ML. This uses DMATRIX=IDENTITY. The guess values are the

* reciprocals of the standard deviations down the diagonal, and zeros off of it.

compute uu=1.0/sqrt(%sigma(1,1)),cc=1.0/sqrt(%sigma(2,2))

compute rr=1.0/sqrt(%sigma(3,3)),xx=1.0/sqrt(%sigma(4,4))

compute mm=1.0/sqrt(%sigma(5,5)),pp=1.0/sqrt(%sigma(6,6))

compute ur=cu=cr=rx=mc=mr=mp=pc=pr=0.0

cvmodel(a=afrml,parmset=aparms,dmatrix=identity,$

method=bfgs,pmethod=simplex,piters=10) %sigma

* This is the ML estimate of the A matrix. You have to use ABASE for this, since

* it gets used in setting up the sampler.

compute abase=afrml(0)

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

* From here to the next full line of ***'s doesn't require change

* Set up mappings that show the free parameters in each row. You have to

* use the name WZSLOTS for this. WZB is the vector of free parameters in each row.

compute nvar=%nvar,wztnobs=%nobs

dec vect[vect[int]] wzslots(nvar)

dec vect[vect] wzb(nvar)

* Build up the wzslots list by examining the non-zero elements in the estimated A

* matrix. Set the WZB vectors by pulling information out of A.

dec vect[int] empty

do i=1,nvar

local list[integer] active

compute active=empty

do j=1,nvar

if abase(i,j)<>0.0

compute active=active+j

end do j

compute wzslots(i)=active

dim wzb(i)(%size(active))

ewise wzb(i)(j)=abase(i,active(j))

end do i

* Create a parameter set with the subvectors out of each row.

nonlin(parmset=allwz) wzb

****

* WZAMatrix reconstructs the A matrix given the current values in the

* WZB vectors

function WZAMatrix

type rect WZAMatrix

local int i j

compute WZAMatrix=%zeros(nvar,nvar)

do i=1,nvar

do j=1,%size(WZSlots(i))

compute WZAMatrix(i,WZSlots(i)(j))=WZB(i)(j)

end do j

end do i

end

****

* WZPerp returns a vector perpendicular to all rows in the A matrix

* except the focus row.

function WZPerp row

type vector WZPerp

type integer row

local rect a

local vect arow

local int i j

local vect[rect] qr

compute a=WZAMatrix()

* Bubble the focus row to the bottom

compute arow=%xrow(a,row)

ewise a(i,j)=%if(i<row,a(i,j),%if(i==nvar,arow(j),a(i+1,j)))

* Do the QR decomp on the transpose and pick the final row

compute qr=%qrdecomp(tr(a))

compute WZPerp=%xcol(qr(1),nvar)

end

****

function WZQForm row

type symmetric WZQForm

type integer row

local integer i j

dim WZQForm(%size(wzslots(row)),%size(wzslots(row)))

ewise WZQForm(i,j)=%sigma(wzslots(row)(i),wzslots(row)(j))

compute WZQForm=inv(WZQForm)

end

****

dec vect WZUW

dec vect wzdraw

dec vect wzbeta(nvar)

infobox(action=define,progress,lower=-nburn,upper=ndraws) "Waggoner-Zha Sampler"

dec series[vect] wzbdraws

gset wzbdraws 1 ndraws = %parmspeek(allwz)

set wzlogl 1 ndraws = 0.0

do draw=-nburn,ndraws

* Do a sweep through the A matrix parameters

do row=1,nvar

compute qrow=%size(wzslots(row))

compute WZTMatrix=%decomp(WZQForm(row))

compute w=WZPerp(row)

dim WZUW(qrow)

ewise WZUW(i)=w(wzslots(row)(i))

compute WZTUW=tr(WZTMatrix)*WZUW

compute w1=wztuw/sqrt(%normsqr(wztuw))

compute wortho=%perp(w1)

dim wzdraw(qrow)

compute wzbeta(1)=sqrt(2.0/wztnobs*%rangamma(.5*(wztnobs+1)))*%if(%ranflip(.5),-1.0,1.0)

compute wzdraw=w1*wzbeta(1)

do i=1,qrow-1

compute wzbeta(i+1)=%ran(1.0/sqrt(wztnobs))

compute wzdraw=wzdraw+wzbeta(i+1)*%xcol(wortho,i)

end do i

compute wzb(row)=WZTMatrix*wzdraw

end do testrow

infobox(current=draw)

if draw<=0

compute wzbdraws(draw)=%parmspeek(allwz)

cvmodel(dmatrix=identity,parmset=allwz,method=evaluate,a=WZAMatrix()) %sigma

compute wzlogl(draw)=%logl

end do idraw

infobox(action=remove)

* Figure out which element within a row has the smallest kurtosis (which would

* mean the most bimodal in sign). Normalize the sign based upon on it.

compute ki=0

do i=1,nvar

compute minkurtosis=%na

do j=1,%size(wzslots(i))

set test 1 ndraws = wzbdraws(t)(ki+j)

stats(noprint) test

if .not.%valid(minkurtosis).or.%kurtosis<minkurtosis

compute minkurtosis=%kurtosis,bestj=j

end do j

?"Equation" i "sign-normalized on" wzslots(i)(bestj)

do t=1,ndraws

compute alla=wzbdraws(t)

compute flip=%sign(alla(ki+bestj))

do j=1,%size(wzslots(i))

compute alla(ki+j)=alla(ki+j)*flip

end do j

compute wzbdraws(t)=alla

end do t

compute ki=ki+%size(wzslots(i))

end do i

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

* WZBDRAWS now has the sign-corrected draws for all the elements. What you do with

* this will obviously depend upon the application.

* This does a scatter graph of the MR and MM coefficients (from the money equation).

set mrx 1 ndraws = wzbdraws(t)(10)

set mmx 1 ndraws = wzbdraws(t)(11)

scatter(footer="Coefficients from Money Demand Equation",$

hlabel="M coefficient",vlabel="R coefficient",smpl=%clock(t,10)==10)

# mmx mrx

* This does a density graph of the normalized (to a unit coefficient on M) of the

* R coefficient in the M equation.

set mrstd 1 ndraws = mrx/mmx

density(smoothing=1.5,grid=automatic,maxgrid=100) mrstd / mrv mrf

scatter(style=line,vmin=0.0,footer="Posterior Density of Normalized MR coefficient")

# mrv mrf

* This does MC integration of the impulse responses taking the keeper draws for

* the A matrix coefficients as given. Since the lag coefficients can be drawn

* directly given A, this doesn't need any burn-in.

declare vect[rect] %%responses(ndraws)

* Re-estimate the model so we can pull out the information required for drawing

* the coefficients.

estimate(noprint)

compute fxx =%decomp(%xx)

compute betaols=%modelgetcoeffs(canmodel)

infobox(action=define,progress,upper=ndraws) "Monte Carlo Integration"

do draw=1,ndraws

* Poke the draw back into the parmset and compute the approximate factor (which

* is the inverse of the A matrix).

compute %parmspoke(allwz,wzbdraws(draw))

compute factor=inv(WZAMatrix())

* Draw the VAR coefficients and put the draw back into the model.

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

compute %modelsetcoeffs(canmodel,betadraw)

impulse(model=canmodel,factor=factor,flatten=%%responses(draw),steps=nsteps,noprint)

infobox(current=draw)

end do draws

infobox(action=remove)

@mcgraphirf(model=canmodel,footer="IRF Error Bands using WZ Sampler",$

shocklabels=||"Real 1","Real 2","Fin 1","Fin 2","Nom 1","Nom 2"||)

Output

Covariance Model-Likelihood - Estimation by BFGS

Convergence in 24 Iterations. Final criterion was 0.0000032 <= 0.0000100

Observations 100

Log Likelihood -581.0251

Log Likelihood Unrestricted -578.6729

Chi-Squared(6) 4.7044

Significance Level 0.5822397

Variable Coeff Std Error T-Stat Signif

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

1. UU 2.498729939 0.176726758 14.13895 0.00000000

2. UR -0.248354254 0.161798078 -1.53496 0.12479263

3. CU -1.272425088 0.265922843 -4.78494 0.00000171

4. CC 2.519950035 0.183451004 13.73637 0.00000000

5. CR 0.369561196 0.168431453 2.19413 0.02822577

6. RR 1.662514287 0.117613062 14.13546 0.00000000

7. RX 0.160339770 0.061437998 2.60978 0.00906000

8. XX 0.602263595 0.042976564 14.01377 0.00000000

9. MC -0.350913490 0.229761667 -1.52729 0.12668806

10. MR 0.271351530 0.164502448 1.64953 0.09903934

11. MM 1.048247458 0.074180565 14.13103 0.00000000

12. MP -0.637429494 0.230929212 -2.76028 0.00577517

13. PC 0.400736908 0.229054274 1.74953 0.08019976

14. PR 0.175737590 0.164250832 1.06993 0.28464891

15. PP 2.259102755 0.163082068 13.85255 0.00000000

Equation 1 sign-normalized on 1

Equation 2 sign-normalized on 2

Equation 3 sign-normalized on 3

Equation 4 sign-normalized on 4

Equation 5 sign-normalized on 5

Equation 6 sign-normalized on 6

Graphs