Identifying VARs with sign restrictions

[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.

[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

[TomDoan]

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

First off, do not start with a NOPRINT option on your ESTIMATE. Look at the output and see if there are problems. Once you're sure that the VAR itself is OK, you can put NOPRINT back on. In your case, it's not OK. ST is the difference between FIRON and PIRON, so you can't include all three in a VAR. You also aren't including a CONSTANT in the system, so the residuals appear to be much more correlated that they are. Offhand, it also appears that you're including in levels data series which would ordinarily be modeled in logs.

Your HSTART is also wrong; it needs to be the start of the estimation period (1986:7 with six lags), rather than the start of the data.

Finally, don't start with four shocks. Start with one. See if you can get that to work successfully. As it is, you're not even getting successful draws for the first shock.

[pei]

Dear Tom,
Thank you very much for your suggestions.
I am a beginner, so there will be many problems cannot understand. Thank you for your patience to answer my questions.
I modified the sample，it works partially. But I don’t know what is worry with it, and how can I correct the mistakes.
Please give me some suggestions.


Sincerely,
pei

[pei]

Dear Tom,
The above is my code, I have modified the code for three shocks and four variables. Nevertheless, results do not appear.
I might do something wrong. Could you please give me some suggestions?


pei

[FRegensburg]

Dear Mr Doan,

I need to calculate a statistic, presented by Mountford and Uhlig in their 2009 article, that allows one to arrive at the impact multipliers of the three policy experiments (as in pages 983-985 in the original paper)

the authors calculate, at different horizons, the impact multipliers, according to this formula:

Multiplier for GDP=(GDP response/initial -lag zero- fiscal shock)/(average fiscal variable share of GDP).

The last ratio is from the raw data, so that's not a problem. But what lines of code do I need add to the replication file code to extract the hard numbers to compute the first ratio?

I also need to calculate the maximum and minimum multipliers of the policy experiments, with confidence intervals (16th and 84th quantiles), with the same formula.

I look forward to hear from you, as I don't know how to do this.

[TomDoan]

Those are already being generated at the end of the mu2009x1.rpf and mu2009b1.rpf programs. They are being graphed by @MCGRAPHIRF. If you want the values to put into a table, use @MCPROCESSIRF instead (or in addition).

[cordura]

Dear Tom,

I am wondering is there any possibilities that in Uhlig code, realizing restriction like following:

Sign restrictions on impulse response functions



It is like for a specific shock, giving different restrictions on different lags. i tried with uhligaccept but seems it is not working.

thank you for your help! Also thanks for your advices about my historical decomposition problem.

cheers.

[TomDoan]

You would have to modify UhligAccept (it would probably make sense to call it something else) to take a VECT[VECT[INT]] where the outer VECT is over steps and the inner over variables. For instance, the loan supply shock would have an input of ||||+1,+2,-4||,||+1,+2,+3,-4||,||+1,+2,+3,-4,+5||||