Identifying VARs with sign restrictions

[TL]

Dear Tom,

Thank you very much for your suggestions.

Sincerely,

Tim

[TL]

Dear Tom,

I have used the following codes with the data set (Jan 1995 – Nov 2008). Impulse responses acquired are unusual. Please see pdf attached.

When I use longer data set (Feb 1989 – Nov 2008), results are fine.

I am wondering whether this is because the data set is too short. Could you please give me suggestions on this?

Thank you very much for your help and kindness.

Sincerely,

Tim

[TomDoan]

Dear Tom,

I have used the following codes with the following data set (Jan 1995 – Nov 2008). Impulse responses acquired are unusual. Please see pdf attached.

When I use longer data set (Feb 1989 – Nov 2008), results are fine.

I am wondering whether this is because the data set is too short. Could you please give me suggestions on this?

Thank you very much for your help and kindness.

Sincerely,

Tim

Offhand, the only thing that might be described as "unusual" is the rapid widening towards the end of the IRF. That's not, in fact, that unusual; if you take almost any model, even one estimated with quite a bit of data, and run the IRF out too far, you'll eventually see that. This is due to the fact that the impulse responses will eventually be dominated by the largest root, which will typically be just a bit larger than 1. With a fairly small data set, the largest root is less well-determined than with a larger one, hence you will see the exploding confidence bands at shorter horizons. If you cut the IRF off at 36 rather than 60, the graphs won't look as odd.

[TL]

Dear Tom,

Thank you very much for your suggestions.

Sincerely,

Tim

[MC128]

Dear Tom,

Do you know how to modify your replication codes for the paper Mountford and Uhlig (2009) such that the size of the basic government spending and revenue shocks are normalised to 1% in the initial period?

Thank you so much!

MC

[luching]

I am also trying to identify multiple shocks using sign restrictions. I think I understand the codes above. But to begin with, I was wondering why it may be useful/necessary to identify multiple shocks SIMULTANEOUSLY. As in, what value additions do we get by identifying the shocks simultaneously instead of identifying each shocks SEQUENTIALLY using the original (one shocks) Uhlig codes? That is, run the Uhlig codes with different sign restrictions for each shock each time.

Any thoughts?

[TomDoan]

The two are quite different. If you run the procedure twice with two different sign restrictions, you isolate a pair of shocks which are very likely to be correlated with each other---there's nothing preventing them from being correlated, and, in fact, there may be nothing preventing them from being the same (it depends upon what restrictions you use). On the other hand, if you use this procedure with multiple sign restrictions, the shocks are uncorrelated by construction, which is what you would like if you're trying to identify separate fundamental shocks.

[luching]

Thanks. This question was motivated from the variance decomposition numbers I got after running these codes, SEQUENTIALLY. I understand when shocks are not orthogonal, the variance decomposition is not strictly valid. However, what I get and surprises me is that the variance decompositions yield a more or less the same number (say 15%) for all variables due to a given shock. The impulse responses however look good. Any thoughts? I was looking for coding errors, but found none.

[TomDoan]

Not strictly valid is an understatement. There's no real justification for using the standard FEVD calculation on a set of non-orthogonal shocks. You can compute the fraction of variance for any one shock in isolation, using any generic method of completing a factorization. If that's what you think you're doing, and are getting odd results, you'll have to send the code and data to support@estima.com so we can look at it.

[luching]

Could you give me some reference on the factorization so that I can look it up? I think this wasn't what I did. I just applied the variance operator which is incorrect. Thanks again!!!

[TomDoan]

That's what the ForcedFactor procedure does---it completes a factorization holding one or more columns fixed (up to scale). It's used extensively in the various Uhlig programs.

[shugo]

Dear Tom,

I have modified your Uhlig(2005)'s sign restriction code with historical decomposition for 4shocks
It works but since I'm a beginner, I might make some mistakes.
If I could have your instruction or checking, I will be very grateful.

Thank you very much.

Code: Select all


OPEN DATA uhligdata.xls
CALENDAR 1965 1 12
compute missc=1.0e+32
data(format=xlS,org=columns) 1965:1 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 6
    deterministic constant
    end(system)
    estimate(noprint,resid=resids)

    dec vect[strings] vl(6)
    compute vl=||'Consumer Prices','Real GDP','Interest Rates',$
          'Market Index Real Stock Prices','Financial Sector Index Real Stock Prices','Real Exchange Rates'||

    compute n1=1000
    compute n2=1000
    compute nvar=6
    compute nstep=20
    compute KMAX=5

    *************************************************************************************
    compute hstart = 1995:1
    compute hend   = 2003:12
    compute nhor   = hend-hstart+1
    dec vect[series] upaths(nvar) u(nvar) onesteps(nvar) base(nvar) baseplus(nvar)
    dec vect[series] upathsa(nvar) baseplusa(nvar)
    dec vect[series] upathsaa(nvar) baseplusaa(nvar)
    dec vect[series] upathsaaa(nvar) baseplusaaa(nvar)
    clear upaths
    clear upathsa
    clear upathsaa
    clear upathsaaa
    *************************************************************************************
    *
    * 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 v3(4)
    dec vect v4(3)
    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) goodfevd(1000)
    declare vect[rect] goodrespa(1000) goodfevda(1000)
    declare vect[rect] goodrespaa(1000) goodfevdaa(1000)
    declare vect[rect] goodrespaaa(1000) goodfevdaaa(1000)
    declare vector ik a(nvar) ones(nvar)
    declare series[rect] irfsquared
    compute ones=%const(1.0)
    source forcedfactor.src

    *************************************************************************************
    declare vector hk(nvar)
    declare vect[rect] goodhist(1000) goodhista(1000) goodhistaa(1000)  goodhistaaa(1000)
    *************************************************************************************

    *
    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)
     *************************************************************************************
        *
        * 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=onesteps,errors=u,static)
        forecast(model=varmodel,from=hstart,to=hend,results=base)
     *************************************************************************************
        *
        * This is changed to unit shocks rather than orthogonalized shocks.
        *
        impulse(noprint,model=varmodel,decomp=p,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 Set of Restrictions - Unique Impulse Vector
        ************************************************
           compute v1=%ransphere(6)
           compute i1=p*v1
           do k=1,KMAX+1
               compute ik=%xt(impulses,k)*v1
               if ik(1)>0.or.ik(2)>0.or.ik(3)<0.or.ik(6)>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=%ransphere(5)
           compute i2=%xsubmat(f,1,6,2,6)*v2
           compute v2=inv(p)*i2
           *****************************************************
           *  Second Set of Restrictions - 2nd Shock
           *****************************************************
           do k=1,KMAX+1
              compute ik=%xt(impulses,k)*v2
              if ik(1)<0.or.ik(2)<0
                 branch 105
           end do k


           @forcedfactor(force=column) sigmad i1~i2 f
           compute v3=%ransphere(4)
           compute i3=%xsubmat(f,1,6,3,6)*v3
           compute v3=inv(p)*i3

           *****************************************************
           *  third Set of Restrictions - 3rd Shock
           *****************************************************
           do k=1,KMAX+1
              compute ik=%xt(impulses,k)*v3
              if ik(1)<0
                 branch 105
           end do k
           @forcedfactor(force=column) sigmad i1~i2~i3 f
           compute v4=%ransphere(3)
           compute i4=%xsubmat(f,1,6,4,6)*v4
           compute v4=inv(p)*i4
           
           *****************************************************
           *  Forth Set of Restrictions - 4th Shock
           *****************************************************
           do k=1,KMAX+1
              compute ik=%xt(impulses,k)*v4
              if ik(1)>0
                 branch 105
           end do k
           
           
           *
           * Meets both restrictions
           *
           compute accept=accept+1
           dim goodrespa(accept)(nstep,nvar) goodfevda(accept)(nstep,nvar)
           dim goodresp(accept)(nstep,nvar) goodfevd(accept)(nstep,nvar)
           dim goodrespaa(accept)(nstep,nvar) goodfevdaa(accept)(nstep,nvar)
           dim goodrespaaa(accept)(nstep,nvar) goodfevdaaa(accept)(nstep,nvar)


           ewise goodresp(accept)(i,j)=(ik=%xt(impulses,i)*v1),ik(j)
           ewise goodfevd(accept)(i,j)=(ik=(irfsquared(i)*(v1.^2))./(irfsquared(i)*ones)),ik(j)
           ewise goodrespa(accept)(i,j)=(ik=%xt(impulses,i)*v2),ik(j)
           ewise goodfevda(accept)(i,j)=(ik=(irfsquared(i)*(v2.^2))./(irfsquared(i)*ones)),ik(j)
           ewise goodrespaa(accept)(i,j)=(ik=%xt(impulses,i)*v3),ik(j)
           ewise goodfevdaa(accept)(i,j)=(ik=(irfsquared(i)*(v3.^2))./(irfsquared(i)*ones)),ik(j)
           ewise goodrespaaa(accept)(i,j)=(ik=%xt(impulses,i)*v4),ik(j)
           ewise goodfevdaaa(accept)(i,j)=(ik=(irfsquared(i)*(v4.^2))./(irfsquared(i)*ones)),ik(j)
     *************************************************************************************
          *
          * Compute the historical paths of shocks to the system corresponding to the
          * <<a>> weights on the Choleski factor.
          *
          dim goodhist(accept)(nhor,nvar) goodhista(accept)(nhor,nvar) goodhistaa(accept)(nhor,nvar)  goodhistaaa(accept)(nhor,nvar)
          compute remap=p*v1*tr(v1)*inv(p)
          do t=hstart,hend
             compute %pt(upaths,t,remap*%xt(u,t))
          end do t
          compute remap=p*v2*tr(v2)*inv(p)
          do t=hstart,hend
             compute %pt(upathsa,t,remap*%xt(u,t))
          end do t
          compute remap=p*v3*tr(v3)*inv(p)
          do t=hstart,hend
             compute %pt(upathsaa,t,remap*%xt(u,t))
          end do t
          compute remap=p*v4*tr(v4)*inv(p)
          do t=hstart,hend
             compute %pt(upathsaaa,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
          forecast(model=varmodel,from=hstart,to=hend,results=baseplusa,paths)
          # upathsa
          forecast(model=varmodel,from=hstart,to=hend,results=baseplusaa,paths)
          # upathsaa
          forecast(model=varmodel,from=hstart,to=hend,results=baseplusaaa,paths)
          # upathsaaa
          *
          * 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)(i,j)=baseplus(j)(i+hstart-1)-base(j)(i+hstart-1)
          ewise goodhista(accept)(i,j)=baseplusa(j)(i+hstart-1)-base(j)(i+hstart-1)
          ewise goodhistaa(accept)(i,j)=baseplusaa(j)(i+hstart-1)-base(j)(i+hstart-1)
          ewise goodhistaaa(accept)(i,j)=baseplusaaa(j)(i+hstart-1)-base(j)(i+hstart-1)
     *************************************************************************************
           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=6,hfields=4,subhea='Impluse Response')
    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=+vl(i)) 3
        # resp
        # upper / 2
        # lower / 2
    end do i
    *
    do i=1,nvar
        compute minlower=maxupper=0.0
        smpl 1 accept
        do k=1,nstep
           set work = goodrespa(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=+vl(i)) 3
        # resp
        # upper / 2
        # lower / 2
    end do i
    *
    do i=1,nvar
        compute minlower=maxupper=0.0
        smpl 1 accept
        do k=1,nstep
           set work = goodrespaa(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=+vl(i)) 3
        # resp
        # upper / 2
        # lower / 2
    end do i
    *
    do i=1,nvar
        compute minlower=maxupper=0.0
        smpl 1 accept
        do k=1,nstep
           set work = goodrespaaa(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=+vl(i)) 3
        # resp
        # upper / 2
        # lower / 2
    end do i
    spgraph(done)
    smpl
    *
    clear upper lower resp

 *
   
    *
    clear upper lower resp
    *
*************************************************************************************
    *
    spgraph(vfields=3,hfields=2,hlabel='Fraction of Variance Explained with Pure-Sign Approach',subhea='Subhea')
    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)
           display resp(k)
        end do k
    *
        smpl 1 nstep
        graph(pattern,ticks,number=0,min=0.0,picture="##.##",header="Fraction Explained for "+vl(i)) 3
        # resp
        # upper / 2
        # lower / 2
    end do i
    *
    spgraph(done)
    smpl
    *
    clear upper lower resp
    *
    spgraph(vfields=3,hfields=2,hlabel='Fraction of Variance Explained with Pure-Sign Approach',subhea='Subhea')
    do i=1,nvar
        compute minlower=maxupper=0.0
        smpl 1 accept
        do k=1,nstep
           set work = goodfevda(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)
           display resp(k)
        end do k
    *
        smpl 1 nstep
        graph(pattern,ticks,number=0,min=0.0,picture="##.##",header="Fraction Explained for "+vl(i)) 3
        # resp
        # upper / 2
        # lower / 2
    end do i
    *
    spgraph(done)
    smpl
    *
    clear upper lower resp
    *
    spgraph(vfields=3,hfields=2,hlabel='Fraction of Variance Explained with Pure-Sign Approach',subhea='Subhea')
    do i=1,nvar
        compute minlower=maxupper=0.0
        smpl 1 accept
        do k=1,nstep
           set work = goodfevdaa(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)
           display resp(k)
        end do k
    *
        smpl 1 nstep
        graph(pattern,ticks,number=0,min=0.0,picture="##.##",header="Fraction Explained for "+vl(i)) 3
        # resp
        # upper / 2
        # lower / 2
    end do i
    *
    spgraph(done)
    smpl

    clear upper lower resp
    *
    spgraph(vfields=3,hfields=2,hlabel='Fraction of Variance Explained with Pure-Sign Approach',subhea='Subhea')
    do i=1,nvar
        compute minlower=maxupper=0.0
        smpl 1 accept
        do k=1,nstep
           set work = goodfevdaaa(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)
           display resp(k)
        end do k
    *
        smpl 1 nstep
        graph(pattern,ticks,number=0,min=0.0,picture="##.##",header="Fraction Explained for "+vl(i)) 3
        # resp
        # upper / 2
        # lower / 2
    end do i
    *
    spgraph(done)
    smpl

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

    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 1 accept
          clear upperH lowerH respH
          do k=hstart,hend
             set work 1 accept = goodhist(t)(k-hstart+1,i)
             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
    *      compute maxupper=%max(maxupper,%maxvalue(upper))  ,min=minlower,max=maxupper
    *      compute minlower=%min(minlower,%minvalue(lower))

    *
    /*
       display 'contributions to variable' +vl(i)
       do h=1,nhor
          display respH(h)
       end do h
    */

       smpl hstart hend
       graph(pattern,ticks,dates,picture="##.##",header="Contribution of 1st shocks to "+vl(i)) 3
       # respH hstart hend
       # upperH / 2
       # lowerH / 2

    end do i
    *
    spgraph(done)
    smpl
    *
    clear 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 1 accept
          clear upperH lowerH respH
          do k=hstart,hend
             set work 1 accept = goodhista(t)(k-hstart+1,i)
             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
    *      compute maxupper=%max(maxupper,%maxvalue(upper))  ,min=minlower,max=maxupper
    *      compute minlower=%min(minlower,%minvalue(lower))

    *
    /*
       display 'contributions to variable' +vl(i)
       do h=1,nhor
          display respH(h)
       end do h
    */

       smpl hstart hend
       graph(pattern,ticks,dates,picture="##.##",header="Contribution of 2nd shocks to "+vl(i)) 3
       # respH hstart hend
       # upperH / 2
       # lowerH / 2

    end do i
    *
    spgraph(done)

clear 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 1 accept
          clear upperH lowerH respH
          do k=hstart,hend
             set work 1 accept = goodhistaa(t)(k-hstart+1,i)
             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
    *      compute maxupper=%max(maxupper,%maxvalue(upper))  ,min=minlower,max=maxupper
    *      compute minlower=%min(minlower,%minvalue(lower))

    *
    /*
       display 'contributions to variable' +vl(i)
       do h=1,nhor
          display respH(h)
       end do h
    */

       smpl hstart hend
       graph(pattern,ticks,dates,picture="##.##",header="Contribution of 3rd shocks to "+vl(i)) 3
       # respH hstart hend
       # upperH / 2
       # lowerH / 2

    end do i
    *
    spgraph(done)

clear upperH lowerH respH

spgraph(vfields=3,hfields=2,hlabel="Figure y. Historical Decomposition with Pure-Sign Approach")
    do i=1,nvar
       compute minlower=maxupper=0.0
       smpl 1 accept
          clear upperH lowerH respH
          do k=hstart,hend
             set work 1 accept = goodhistaaa(t)(k-hstart+1,i)
             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
    *      compute maxupper=%max(maxupper,%maxvalue(upper))  ,min=minlower,max=maxupper
    *      compute minlower=%min(minlower,%minvalue(lower))

    *
    /*
       display 'contributions to variable' +vl(i)
       do h=1,nhor
          display respH(h)
       end do h
    */

       smpl hstart hend
       graph(pattern,ticks,dates,picture="##.##",header="Contribution of 4th shocks to "+vl(i)) 3
       # respH hstart hend
       # upperH / 2
       # lowerH / 2

    end do i
    *
    spgraph(done)






[TomDoan]

Dear Tom,

I have modified your Uhlig(2005)'s sign restriction code with historical decomposition for 4shocks
It works but since I'm a beginner, I might make some mistakes.
If I could have your instruction or checking, I will be very grateful.

Thank you very much.

That looks correct.

[shugo]

Dear Tom,

Thank you very much!

[pei]

Dear Tom,

I use code for four shocks from the forum.I am not quite sure whether the code is correct. Could you please give me some suggestions and feedbacks whether my code is correct?

Thank you very much for your help and kindness.


pei