Identifying VARs with sign restrictions

[econdawg]

Note: most of this thread is now obsolete because the Mountford-Uhlig replication shows how to handle multiple shocks.

Hello,

RATS provides a code for replicating Uhlig(2005)'s paper which essentially identify VARs using sign restrictions.The code though allows for the identification of one shock. I was thinking about identifying more than one shock at the same time (say 2 shocks). i am kinda puzzled as to how to incorporate that in the code. Would identifying a second shock implies doing a second inner-loop subdraw for another vector over the unit sphere? That is i would still draw initially from the posterior of the VAR then compute impulse responses, but instead of doing one subdraw (and compute vector a) for each draw from the posterior, i would do two subdraws (and compute a and another vector--say b). Then i would be essentially be constructing 2 impulse vectors simultaneously and go from there to check sign restrictions to identify the two shocks of interests. Am i thinking right about this? it would be very helpful (if any of you have time) to maybe guide me as to how i would modify the code in order to identify that second shock. Thanks a lot.

[TomDoan]

Hello,

RATS provides a code for replicating Uhlig(2005)'s paper which essentially identify VARs using sign restrictions.The code though allows for the identification of one shock. I was thinking about identifying more than one shock at the same time (say 2 shocks). i am kinda puzzled as to how to incorporate that in the code. Would identifying a second shock implies doing a second inner-loop subdraw for another vector over the unit sphere? That is i would still draw initially from the posterior of the VAR then compute impulse responses, but instead of doing one subdraw (and compute vector a) for each draw from the posterior, i would do two subdraws (and compute a and another vector--say b). Then i would be essentially be constructing 2 impulse vectors simultaneously and go from there to check sign restrictions to identify the two shocks of interests. Am i thinking right about this? it would be very helpful (if any of you have time) to maybe guide me as to how i would modify the code in order to identify that second shock. Thanks a lot.

This shows an example of the subdraw procedure for doing multiple mutually orthogonal sign-restricted impulse vectors. This relatively simple logic will only work if no impulse vector can meet more than one set of restrictions. (Here ik(2) needs to be positive for the first set and negative for the second). If the restrictions don't separate the space like that (they probably should if the whole idea is to make sense), you would have to check to see if either of the restrictions are satisfied by your first draw, then see if the remaining restriction is met by the second.

Code: Select all

   do subdraws=1,n2
   ************************************************
   * First Set of Restrictions - Unique Impulse Vector
   ************************************************
      compute v1=%ran(1.0),v1=v1/sqrt(%normsqr(v1))
      compute v=i1=p*v1
      do k=1,KMAX+1
          compute ik=%xt(impulses,k)*v
          if ik(4)<0.or.ik(6)<0.or.ik(2)>0.or.ik(5)>0
             branch 105
      end do k
      *
      * Meets the first restriction
      * Draw from the orthogonal complement of i1 (last five columns of the
      * factor "f").
      *
      @forcedfactor(force=column) sigmad i1 f
      compute v2=%ran(1.0),v2=v2/sqrt(%normsqr(v2))
      compute v=i2=%xsubmat(f,1,6,2,6)*v2
      *****************************************************
      *  Second Set of Restrictions - Demand Shock
      *****************************************************
      do k=1,KMAX+1
         compute ik=%xt(impulses,k)*v
         if ik(1)<0.or.ik(2)<0
            branch 105
      end do k
      *
      * Meets both restrictions
      *
      compute accept=accept+1
      dim goodrespa(accept)(nstep,nvar)
      dim goodresp(accept)(nstep,nvar)
      ewise goodresp(accept)(i,j)=(ik=%xt(impulses,i)*i1),ik(j)
      ewise goodrespa(accept)(i,j)=(ik=%xt(impulses,i)*i2),ik(j)
      if accept>=1000
         break
   :105
   end do subdraws

[econdawg]

Thank you very much Tom for your help. And I totally see your point.

P.S what is exactly p here (as in "v=i1=p*v1") ?

[TomDoan]

p is (any) factor of the inverse Wishart draw for the covariance matrix.

[lhlee0506]

Hi~Tom,

I am confused about these two lines " @forcedfactor(force=column) sigmad i1 f"& "compute v=i2=%xsubmat(f,1,6,2,6)*v2".
If you could tell me their meaning and why we need them, I will be very appreciated.
Thank you very much!

[TomDoan]

@forcedfactor(force=column) sigmad i1 f

SIGMAD is the draw for the covariance matrix. I1 is the first identified shock. This computes a matrix F (there are many of them), for which FF'=SIGMAD and the first column of F is a scale multiple of I1. u=Fv is an orthogonalized representation of the residuals and I1 is proportional to F x ||1,0,0,...0||', that is, I1 is the first component in v. We want to identify a orthogonal component in the factorization. The space of those is the span of columns 2 to 6 of F.

compute v=i2=%xsubmat(f,1,6,2,6)*v2

v2 is a random direction in 5-space. 0.0~v2 is a random direction in 6-space (thus a potential impulse vector) that's orthogonal to I1 since it has a zero component in the 1st component. A simpler (and equivalent) way of looking at this would be

Code: Select all

compute v2=%ranmat(6,1)            ;* Draw a random 6-vector
compute v2(1)=0.0                  ;* Zero the 1st component
compute v2=v2/sqrt(%normsqr(v2))   ;* Normalize to unit length
compute v=i2=f*v2                  ;* Transform to impulse vector

[lhlee0506]

Dear Tom,

Tahnk you for your kind instruction.
But when I ran this program as following, it showed the error message like

## SX11. Identifier P is Not Recognizable. Incorrect Option Field or Parameter Order?
>>>> compute v=i1=p*v<<<<

## SX11. Identifier GOODRESPA is Not Recognizable. Incorrect Option Field or Parameter Order?
>>>> dim goodrespa(<<<<
why is that? did I miss any point?
If I could have your instruction, i will be very grateful.

Code: Select all

*
*  Replication File for 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. Pure sign restriction approach
*
open data uhligdata.xls
calendar 1965 1 12
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=200
compute n2=200
compute nvar=6
compute nstep=60
compute KMAX=5
*
* 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
dec rect ranc(ncoef,nvar)
*
* Most draws are going to get rejected. We allow for up to 1000
* good ones. The variable accept will count the number of accepted
* draws. GOODRESP will be a RECT(nsteps,nvar) at each accepted
* draw.
*
declare vect[rect] goodresp(1000)
declare vector ik a(nvar)
*
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)
   compute swish   =%decomp(sigmad)
   compute ranc    =%ran(1.0)
   compute betau   =sxx*ranc*tr(swish)
   compute betadraw=betaols+betau
   compute %modelsetcoeffs(varmodel,betadraw)
   impulse(noprint,model=varmodel,decomp=swish,results=impulses,steps=nstep)
   *
   * Do the subdraws over the unit sphere. These give the weights on the
   * orthogonal components.
   *
   do subdraws=1,n2
   ************************************************
   * First Set of Restrictions - Unique Impulse Vector
   ************************************************
      compute v1=%ran(1.0),v1=v1/sqrt(%normsqr(v1))
      compute p=inv(sigmad)
      do k=1,KMAX+1
          compute ik=%xt(impulses,k)*v
          if ik(4)<0.or.ik(6)<0.or.ik(2)>0.or.ik(5)>0
             branch 105
      end do k
      *
      * Meets the first restriction
      * Draw from the orthogonal complement of i1 (last five columns of the
      * factor "f").
      *
      @forcedfactor(force=column) sigmad i1 f
      compute v2=%ran(1.0),v2=v2/sqrt(%normsqr(v2))
      compute v=i2=%xsubmat(f,1,6,2,6)*v2

      *****************************************************
      *  Second Set of Restrictions - Demand Shock
      *****************************************************
      do k=1,KMAX+1
         compute ik=%xt(impulses,k)*v
         if ik(1)<0.or.ik(2)<0
            branch 105
      end do k
      *
      * Meets both restrictions
      *
      compute accept=accept+1
      dim goodrespa(accept)(nstep,nvar)
      dim goodresp(accept)(nstep,nvar)
      ewise goodresp(accept)(i,j)=(ik=%xt(impulses,i)*i1),ik(j)
      ewise goodrespa(accept)(i,j)=(ik=%xt(impulses,i)*i2),ik(j)
      if accept>=1000
         break
   :105
   end do subdraws

   if accept>=1000
      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)
      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)

[TomDoan]

This is the full running program from which the code above was taken. It uses different names for a few of the matrices than the Uhlig example, hence the confusion. The P in this is the same as the SWISH in your program.

Code: Select all

OPEN DATA fiscal.RAT
CALENDAR(q) 1980
ALL 2005:04
compute missc=1.0e+32
DATA(FORMAT=RATS,missing=missc) / GDP GOVE DEFLA CRUDE INT M1

set gdpc1    = log(gdp)*100.0
set gdpdef   = log(defla)*100.0
set cprindex = log(crude)*100.0
set bognonbr = log(m1)*100.0
set eerr = log(gove)*100.0
*
system(model=varmodel)
variables gdpc1 gdpdef cprindex int bognonbr gove
lags 1 to 4
end(system)
estimate(noprint,resid=resids)

dec vect[strings] vl(6)
compute vl=||'real GDP','GDP price defl','Comm. Price Ind.',$
     'Bank Rate','Money Stock M1','Real Effective Exchange Rate'||

compute n1=200
compute n2=200
compute nvar=6
compute nstep=40
compute KMAX=1

* This is the standard setup for MC integration of an OLS VAR
*
dec symm s(6,6)
dec vect v1(6)		;* For the unit vector on the 1st draw
dec vect v2(5)	   ;* For the unit vector on the 2nd draw
dec vect v(6)     ;* Working impulse vector

compute sxx    =%decomp(%xx)
compute svt    =%decomp(inv(%nobs*%sigma))
compute betaols=%modelgetcoeffs(varmodel)
compute ncoef  =%rows(sxx)
compute wishdof=%nobs-ncoef
dec rect ranc(ncoef,nvar)
*
* Most draws are going to get rejected. We allow for up to 1000
* good ones. The variable accept will count the number of accepted
* draws. GOODRESP will be a RECT(nsteps,nvar) at each accepted
* draw.
*
declare vect[rect] goodresp(1000)
declare vect[rect] goodrespa(1000)

declare vector ik a(nvar)
source forcedfactor.src
*
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)
   compute p       =%decomp(sigmad)
   compute ranc    =%ran(1.0)
   compute betau   =sxx*ranc*tr(p)
   compute betadraw=betaols+betau
   compute %modelsetcoeffs(varmodel,betadraw)
   *
   * This is changed to unit shocks rather than orthogonalized shocks.
   *
   impulse(noprint,model=varmodel,decomp=%identity(6),results=impulses,steps=nstep)
   *
   * Do the subdraws over the unit sphere. These give the weights on the
   * orthogonal components.
   *
   do subdraws=1,n2
   ************************************************
   * First Set of Restrictions - Unique Impulse Vector
   ************************************************
      compute v1=%ran(1.0),v1=v1/sqrt(%normsqr(v1))
      compute v=i1=p*v1
      do k=1,KMAX+1
          compute ik=%xt(impulses,k)*v
          if ik(4)<0.or.ik(6)<0.or.ik(2)>0.or.ik(5)>0
             branch 105
      end do k
      *
      * Meets the first restriction
      * Draw from the orthogonal complement of i1 (last five columns of the
      * factor "f").
      *
      @forcedfactor(force=column) sigmad i1 f
      compute v2=%ran(1.0),v2=v2/sqrt(%normsqr(v2))
      compute v=i2=%xsubmat(f,1,6,2,6)*v2
      *****************************************************
      *  Second Set of Restrictions - Demand Shock
      *****************************************************
      do k=1,KMAX+1
         compute ik=%xt(impulses,k)*v
         if ik(1)<0.or.ik(2)<0
            branch 105
      end do k
      *
      * Meets both restrictions
      *
      compute accept=accept+1
      dim goodrespa(accept)(nstep,nvar)
      dim goodresp(accept)(nstep,nvar)
      ewise goodresp(accept)(i,j)=(ik=%xt(impulses,i)*i1),ik(j)
      ewise goodrespa(accept)(i,j)=(ik=%xt(impulses,i)*i2),ik(j)
      if accept>=1000
         break
   :105
   end do subdraws
   if accept>=1000
      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='Impulse Responses with Pure-Sign Approach',subhea='Italy')
do i=1,nvar
   compute minlower=maxupper=0.0
   smpl 1 accept
   do k=1,nstep
      set work = goodresp(t)(k,i)
      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(pattern,ticks,number=0,picture='##.##',header='Impulse Responses for '+vl(i)) 3
   # resp
   # upper / 2
   # lower / 2
end do i
*
spgraph(done)

[lhlee0506]

Dear Tom,

How do we show the impluse response of shock seperately in your 2-shocks case?
and how to show the value of impluse response to allow us calculate the impacr multipliers?
Many thanks for your helps.

[TomDoan]

Code: Select all

      dim goodrespa(accept)(nstep,nvar)
      dim goodresp(accept)(nstep,nvar)
      ewise goodresp(accept)(i,j)=(ik=%xt(impulses,i)*i1),ik(j)
      ewise goodrespa(accept)(i,j)=(ik=%xt(impulses,i)*i2),ik(j)

This is a snippet out of the earlier program. It defines both goodresp, which are the responses to the accepted first schocks, and goodrespa which are the responses to the accepted second shocks. (The pair has to pass all tests in order to get there). If you just take the graphing code later and do the same thing with "goodrespa" replacing "goodresp", you'll get the responses to the second shock.

The impact multipliers are just the 1st entries of the goodresp's. goodresp(d)(1,j) is the impact response of variable "j" to the first shock for accepted draw "d", and goodrespa(d)(1,j) is the same for the second shock.

[lhlee0506]

Dear Tom,

Thank you for your kind instruction; I have successful modified my code.
But there is another confusing problem that the result seems not very stable.
When I run the code in second time, monete carlo integration will result in different impulse response graphs.
How to solve this problem?
If I wanna do further variation such as identify a year delayed shock, how do I reach it?
I will be very grateful to have your general response.

[TomDoan]

Regarding the instability, you may need to send us the code and data (to support@estima.com). However, I'm not surprised that it might not behave well. You might want to look at what percentage of draws end up working. I'm guessing that it could be quite small.

While it's possible to set this up to hit a true zero at some point beyond period 0 (you can construct a restriction matrix based upon the period "d" responses), I'm not sure that that's likely to be a fruitful exercise. After all, big positive at d-1, 0 at d, big positive at d+1 and big positive at d-1, 0 at d, big negative at d+1 would both pass, even though they're clearly very different. Since you don't have enough freedom to force 0's from d on, just hitting the one point doesn't really identify the behavior very well.

It's possible that you really should be looking at conditional forecasting, which uses restrictions at shocks at multiple steps to constrain the behavior of the data going forward.

[fabio fornari]

Hi Tom
thanks for the extension of the sign restriction to 2 impulse responses.
One question. I have three shocks to identify, then I can simply extend the process you did?
That is for the third shocks i could

comp v3 = %rannmat(whatever)

and then

comp v3(1) = 0
comp v3(2) = 0

...
and finally comp v=i3=f*v3 where f comes from forcefactor (the same f i used to identify the first ortghoonal complemet of i1, as called in your code)?

Thanks a lot for support.

F

[TomDoan]

That's correct. This is a (standalone) example which shows the process of drawing multiple orthogonal shocks. (It's abstracted from the process of determining whether the shocks actually meet your criteria).

Note that this uses the %RANSPHERE function that was introduced in version 6.35 precisely to simplify the process for drawing random directions in N space. v=%RANSPHERE(n) is the equivalent of v=%RANMAT(n,1),v=v/sqrt(%normsqr(v))

Code: Select all

dec symm s(4,4)
dec vect v(4)
compute s=%ranwishart(4,10.0)
compute p=%decomp(s)
*
* Draw an impulse vector
*
compute v=%ransphere(4)
compute i1=p*v
*
* Do a forcedfactor
*
@forcedfactor(force=column) s i1 f
*
* Draw an impulse vector from the span of the last three columns
*
compute v=%ransphere(3)
compute i2=%xsubmat(f,1,4,2,4)*v
*
* Do a forced factor with i1/i2
*
@forcedfactor(force=column) s i1~i2 f
*
* Draw an impulse vector from the span of the last two columns
*
compute v=%ransphere(2)
compute i3=%xsubmat(f,1,4,3,4)*v
*
* This will give a complete factorization of f
*
@ForcedFactor(force=column) s i1~i2~i3 f

[TomDoan]

There's an article in the forthcoming newsletter on this very topic. See

http://www.estima.com/newslett/March2009RATSLetter.pdf