sign restrictions: Rubio-Waggoner-Zha approach

[tclark]

In light of the broad interest in sign restrictions on VARs, before writing my own code I thought it might be worth asking if anybody has already written code for the sign restriction algorithm developed (see pp.45-46) in "Structural Vector Autoregressions: Theory of Identification and Algorithms for Inference," by J. Rubio-Ramirez, D. Waggoner, and T. Zha (http://home.earthlink.net/~tzha01/worki ... TIFICATION). The same algorithm is used in the paper "Markov-Switching Structural Vector Autoregressions: Theory and Application" by the same authors.

If anyone has code already and would be willing to share it, I would very much appreciate it.

Todd

[TomDoan]

There's a %QRDECOMP function in RATS, so you can do that with

Code: Select all

compute qr=%qrdecomp(%ranmat(n,n))

The "Q" matrix is qr(1). %QRDECOMP doesn't control the signs on the diagonal of the R matrix, but that doesn't seem to be important. (Flipping signs on columns of Q doesn't affect its distribution).

Note, by the way, that that procedure is actually no different (in practice) from the sequential generation that's done in our Uhlig examples.

[tclark]

Thanks, Tom. The coding seems easy enough now.

[tclark]

The program in the box modifies Tom Doan's Uhlig-with historical decompositions program to use the sign restrictions algorithm of Rubio, Waggoner, and Zha. As Tom noted, the RWZ approach is the same in practice, except that the RWZ approach is faster with large models. This code handles the impulse responses and historical decompositions (which I need), but not the forecast error variance decomposition (which is simply commented out).

Code: Select all

*
* This program applies the sign restriction algorithm of "Structural Vector Autoregressions: Theory of Identification 
* and Algorithms for Inference" and "Markov-Switching Structural Vector Autoregressions: Theory and Application"
* by J. Rubio-Ramirez, D. Waggoner, and T. Zha to estimate impulse responses as in Uhlig (2005), "What are the
* effects of monetary policy on output? Results from an agnostic identification procedure." Journal of Monetary
* Economics, 52, pp 381-419. 
*
* Todd Clark wrote the program by modifying Tom Doan's Uhlig code.
*
* This version omits the computation of the forecast error variance but does inclue an historical decomposition 
* of the data showing the cumulative effects
* of the monetary policy shocks. Note that this is not part of the paper.
*
grparm(portrait) header 15 subheader 12 keylabeling 16 footer 10 vlabel 12 hlabel 12
env nowshowgraphs gsave="Myresults*.rgf" gformat=rgf

open data uhligdata.xls
calendar(m) 1965
data(format=xls,org=columns) 1965:01 2003:12 gdpc1 gdpdef cprindex totresns bognonbr fedfunds
*
set gdpc1 = log(gdpc1)*100.0
set gdpdef = log(gdpdef)*100.0
set cprindex = log(cprindex)*100.0
set totresns = log(totresns)*100.0
set bognonbr = log(bognonbr)*100.0
*
system(model=varmodel)
variables gdpc1 gdpdef cprindex fedfunds bognonbr totresns
lags 1 to 12
end(system)
estimate(noprint)
*
dec vect[strings] vl(6)
compute vl=||"real GDP","GDP price defl","Comm. Price Ind.","Fed Funds Rate","Nonborr. Reserv.","Total Reserves"||
*
* n1 is the number of draws from the posterior of the VAR
* n2 is the number of draws from the unit sphere for each draw for the VAR
* nvar is the number of variables
* nstep is the number of IRF steps to compute
* KMAX is the "K" value for the number of steps constrained
*
compute n1=10000
compute n2=10000
compute nkeep=20000
compute nvar=6
compute nstep=60
compute KMAX=5

* set up the number of shocks considered, along with signs
* Note that the computations portion of the program is set up to allow
*  for multiple shocks, but the portion of the program that constructs the charts
*  assumes only one shock of interest
comp nshock = 1
dec rec signs(nvar,nshock)   ;* a matrix with the signs to be imposed, for # shocks = nshock
                              * use a 1 for >= 0, -1 for <= 0, and %NA for unrestricted
comp signs = ||%NA|-1|-1|1|-1|%NA||

*************************************************************************************
compute hstart = 1998:1
compute hend = 2002:12
compute nhor = hend-hstart+1
dec vect[series] upaths(nvar) u(nvar) onesteps(nvar) base(nvar) baseplus(nvar)
clear upaths
*************************************************************************************
*
* This is the standard setup for MC integration of an OLS VAR
*
compute sxx =%decomp(%xx)
compute svt =%decomp(inv(%nobs*%sigma))
compute betaols=%modelgetcoeffs(varmodel)
compute ncoef =%rows(sxx)
compute wishdof=%nobs-ncoef
*
* Most draws are going to get rejected. We allow for up to nkeep good ones. The
* variable accept will count the number of accepted draws. GOODRESP will be a
* RECT(nvar,nshock) at each accepted draw.
*
declare rect[rect] goodresp(nkeep,nstep) ;* goodfevd(nkeep)
declare vector ik a ones(nvar)
dec rec ikmat
declare series[rect] irfsquared
compute ones=%const(1.0)

*************************************************************************************
declare vector hk(nvar)
declare rect[rect] goodhist(nkeep,nshock)
*************************************************************************************
* set up a vector of selection matrices that correspond to e(i) e(i)' in Tom Doan's notes on hist. decomp.
dec vec[rec] selmat(nshock)
do i = 1,nshock
 comp selmat(i) = %zeros(nvar,nvar)
 comp selmat(i)(i,i) = 1.
end do i

*
compute accept=0
infobox(action=define,progress,lower=1,upper=n1) "Monte Carlo Integration"
do draws = 1,n1
  *
  * Make a draw from the posterior for the VAR and compute its impulse
  * responses.
  *
  compute sigmad =%ranwisharti(svt,wishdof)
  comp [symm] invsigmad = inv(sigmad)
  compute swish =%decomp(sigmad)
  compute betau =%ranmvkron(swish,sxx)
  compute betadraw=betaols+betau
  compute %modelsetcoeffs(varmodel,betadraw)
  **********************************************************************************
  *
  * Use static forecasts over the period for historical decomposition to get the
  * residuals. Use dynamic forecasts over the same period to get the base forecast
  * series for the historical decomposition.
  *
  forecast(model=varmodel,from=hstart,to=hend,results=onestep,errors=u,static)
  forecast(model=varmodel,from=hstart,to=hend,results=base)
  **********************************************************************************
  *
  impulse(noprint,model=varmodel,decomp=swish,results=impulses,steps=nstep)
  *gset irfsquared 1 1 = %xt(impulses,t).^2
  *gset irfsquared 2 nstep = irfsquared{1}+%xt(impulses,t).^2
  *
  * Do the subdraws over the unit sphere. These give the weights on the
  * orthogonal components.
  *
  do subdraws = 1,n2
    * first compute the QR decomposition and normalize signs as in RWZ Matlab code
    comp qr=%qrdecomp(%ranmat(nvar,nvar))
    do i = 1,nvar
     if qr(2)(i,i)<0.
      comp %psubmat(qr(1),1,i,-1.*%xcol(qr(1),i))
    end do i
    comp pmat = qr(1)
    *
    * Check that the responses have the correct signs for steps 1 to KMAX+1
    * (+1 because in the paper, the steps used 0-based subscripts, rather than
    * the 1 based subscripts used by RATS).
    *
    do k=1,KMAX+1
     comp ikmat = %xt(impulses,k)*pmat  ;* matrix of responses of i (row) to shock to j (col)
     do i = 1,nshock  ;* loop across shocks to check signs
      comp ik = %xcol(ikmat,i).*%xcol(signs,i)  ;* if the sign is right, the product should be >= 0
      if %minvalue(ik)<0
       branch 105
     end do i
    end do k
    *
    * This is an accepted draw. Copy the information out. If we have enough
    * good ones, drop out of the overall loop.
    *
    compute accept=accept+1
    dim goodresp(accept,nstep)(nvar,nshock) ;*goodfevd(accept)(nstep,nvar)
    do k = 1,nstep
      comp ikmat = %xt(impulses,k)*pmat  ;* matrix of responses of i (row) to shock to j (col), at horizon k
      comp goodresp(accept,k) = %xsubmat(ikmat,1,nvar,1,nshock)  ;* pull out just responses to shocks of interest (subset of cols)
    end do i
    *ewise goodresp(accept)(i,j)=ik=%xt(impulses,i)*a,ik(j)
    *ewise goodfevd(accept)(i,j)=ik=(irfsquared(i)*(a.^2))./(irfsquared(i)*ones),ik(j)
    *******************************************************************************
    *
    * Compute the historical paths of shocks to the system corresponding to the
    * <<a>> weights on the Choleski factor.
    *
    comp PAmat = swish*pmat
    do k = 1,nshock
     dim goodhist(accept,k)(nhor,nvar)
     compute remap = %mqform(selmat(k),tr(PAmat))*invsigmad ;* this uses the middle term of equation 4 of Tom Doan's notes
                                                             * on historical decompositions, with F = PAmat
     do t=hstart,hend
      compute %pt(upaths,t,remap*%xt(u,t))
     end do t
     *
     * Forecast over the decomposition period with those added shocks
     *
     forecast(model=varmodel,from=hstart,to=hend,results=baseplus,paths)
     # upaths
     *
     * Save the difference between the forecasts with and without and added shocks
     * (and shift them into the slots 1,...,nhor) of <<goodhist>>.
     *
     ewise goodhist(accept,k)(i,j)=baseplus(j)(i+hstart-1)-base(j)(i+hstart-1)
    end do k
    *******************************************************************************
    if accept>=nkeep
     break
  :105
  end do subdraws
  if accept>=nkeep
    break
  infobox(current=draws)
end do draws

infobox(action=remove)
*
* Post-processing. Graph the mean of the responses along with the 16% and 84%-iles
*
clear upper lower resp
*
spgraph(vfields=3,hfields=2,hlabel="Figure 6. Impulse Responses with Pure-Sign Approach")

do i=1,nvar
 compute minlower=maxupper=0.0
 smpl 1 accept
 do k=1,nstep
  set work = goodresp(t,k)(i,1)   ;* only 1 shock of interest, so we only pull out the first column
  compute frac=%fractiles(work,||.16,.84||)
  compute lower(k)=frac(1)
  compute upper(k)=frac(2)
  compute resp(k)=%avg(work)
 end do k

 smpl 1 nstep
 graph(ticks,number=0,picture="##.##",header="Impulse Responses for "+vl(i)) 3
 # resp
 # upper / 2
 # lower / 2
end do i

spgraph(done)

/*
spgraph(vfields=3,hfields=2,hlabel="Figure 8. Fraction of Variance Explained with Pure-Sign Approach")

do i=1,nvar
 compute minlower=maxupper=0.0
 smpl 1 accept
 do k=1,nstep
  set work = goodfevd(t)(k,i)
  compute frac=%fractiles(work,||.16,.50,.84||)
  compute lower(k)=frac(1)
  compute upper(k)=frac(3)
  compute resp(k)=frac(2)
end do k

 smpl 1 nstep
 graph(ticks,number=0,min=0.0,picture="##.##",header="Fraction Explained for "+vl(i)) 3
 # resp
 # upper / 2
 # lower / 2
end do i

spgraph(done)
*/

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

declare series upperH lowerH respH
*
spgraph(vfields=3,hfields=2,hlabel="Figure X. Historical Decomposition with Pure-Sign Approach")

do i=1,nvar
 compute minlower=maxupper=0.0
 smpl hstart hend
 clear upperH lowerH respH
 smpl 1 accept
 do k=hstart,hend
  set work 1 accept = goodhist(t,1)(k-hstart+1,i) ;* only 1 shock of interest, so we only pull out the first column
  compute frac=%fractiles(work,||.16,.50,.84||)
  compute lowerH(k)=frac(1)
  compute upperH(k)=frac(3)
  compute respH(k)=frac(2)
 end do k

 smpl hstart hend
 graph(ticks,dates,picture="##.##",header="Contribution of Monetary Policy shock to "+vl(i)) 3
 # respH hstart hend
 # upperH / 2
 # lowerH / 2
end do i

spgraph(done)

[marko3499]

i would like to do forward error variance decomposition with this code, but I’m not quite sure how to do this. As far as I can see, the goodfevd matrix would be included in a loop similar to goodresp in Todd Clark’s code (where ikmatv would be the holder for FEVD below):

dim goodresp(accept,nstep)(nvar,nshock) goodfevd(accept,nstep)(nvar,nshock)
do k = 1,nstep
comp ikmat = %xt(impulses,k)*pmat ;* matrix of responses of i (row) to shock to j (col)
comp goodresp(accept,k) = %xsubmat(ikmat,1,nvar,1,nshock) ;* pull out just responses to shocks of interest (subset of cols)
comp ikmatv = ???
comp goodfevd(accept,k) = %xsubmat(ikmatv,1,nvar,1,nshock)
end do i

How could this be done so that goodfevd could then be used for graphs analogously to goodresp?

Thanks in advance!

[istiak]

Dear Todd
From your code posted on October 27, 2009; I can understand that the VAR model uses monthly data from 1965.01-2003.12. But the code also shows the following part:

compute hstart = 1998:1
compute hend = 2002:12

My question is, what does 'hstart' and 'hend' stand for.

Thank you so much.

KI

[TomDoan]

Dear Todd
From your code posted on October 27, 2009; I can understand that the VAR model uses monthly data from 1965.01-2003.12. But the code also shows the following part:

compute hstart = 1998:1
compute hend = 2002:12

My question is, what does 'hstart' and 'hend' stand for.

Thank you so much.

KI

That's the range for the historical decomposition.

[istiak]

Dear Todd/Tom
I will really appreciate the following two helps from the code.

1. In the code we are getting impulse responses for 1 SD increase in fedrate. Can you please show how I can get impulse responses for 1% increase in fedrate?
2. Besides the IRFs from 1% increase in monetary shock, how can I get IRFs from 1% increase in total reserve (the fourth variable in the VAR model)?

Thank you so much.
Regards,

KI

[KOBE24]

Dear Tom,

I am using the above program by Todd Clark to achieve historical decomposition after identifying 4 shocks.
The strange thing, at least to me, is that the baseline looks to be different depending on the number of draws from the posterior and / or from the unit sphere (I mean, very different, if I take
whatever n1 or n2 belonging to the following intervals

n1 [10000, 50000]
n2 [5000, 10000]

Is it supposed to be like that? What am I missing?
My intuition is that only baseplus(j)(i+hstart-1) should be changing across draws, not base(j)(i+hstart-1).

Thank you in advance for taking the time of answer this (basic) question.

Code: Select all

*******************************************************************************
         ** Forecast over the decomposition period with those added shocks
*******************************************************************************

         forecast(model=varmodel,from=hstart,to=hend,results=baseplus,paths)
         # upaths

         * Save the difference between the forecasts with and without and added shocks
         * (and shift them into the slots 1,...,nhor) of <<goodhist>>.

         ewise goodhist(accept,k)(i,j) = baseplus(j)(i+hstart-1)-base(j)(i+hstart-1)

         end do k

         *******************************************************************************
         if accept>=nkeep
         break
         :105
      end do subdraws

      if accept>=nkeep
        break
      infobox(current=draws)
}
end do draws

infobox(action=remove)

[TomDoan]

base changes with each value of draw. Given that, it's fixed across each subdraw.

[KOBE24]

thanks a lot, Tom.

your help is invaluable!

[lali62]

Hi,
Could you tell me if I understand this logic of branch 105 in the loop correctly:
lets say I have 3 shocks and KMAX=5
do subdraws = 1,n2
.
.
.
do k=1,KMAX+1
comp ikmat = %xt(impulses,k)*pmat ;* matrix of responses of i (row) to shock to j (col)
do i = 1,nshock ;* loop across shocks to check signs
*disp %minvalue(%xcol(ikmat,2).*%xcol(signs3,2))
comp ik = %xcol(ikmat,i).*%xcol(signs,i) ;* if the sign is right, the product should be >= 0
if %minvalue(ik)<0
branch 105
end do i
end do k

if %minvalue(ik)<0, I iterate through next value of nshock.
lets say the value of i=3 and %minvalue(ik)<0, then I would go to the next iteration of k.
Also lets say if I am at i=3 and k=6, and %minvalue(ik)<0, then I would go to next iteration of subdraw.
Is this understanding correct.
if %minvalue(ik)>=0, I break out of inner loops i and k,and treat this as an accepted draw

[TomDoan]

The point of this is to reject the subdraw completely whenever the calculation fails to produce the proper sign (>=0).

if %minvalue(ik)<0
branch 105

It's only if you fully complete the I and K loops that you get an accepted subdraw.