Identifying VARs with sign restrictions

[TomDoan]

Dear Tom,

I use code for two shocks from the forum. In addition to the impulse responses, I would like to see how much of the variation is accounted for by the first shock and by the second shock (variance decomposition).

I am not quite sure whether fraction explained (variance decomposition) 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.

Tim

That won't work without changing the way the impulse responses are computed. In the original Uhlig code, the impulse response and squares are computed by:

Code: Select all

   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

The shocks are computed using a factorization of the covariance matrix (decomp=swish). Because the shocks are orthogonal (across time by assumption, across variables by construction), the sums of squared impulse responses (in IRFSQUARED) provide all the information needed for computing the decomposition of variance. To get the variance produced by impulse vector swish * a (where ||a||=1) you just need the sum over j of irfsquared(i,j)*a(j)^2, which is what (irfsquared(i)*(a.^2)) is computing in:

Code: Select all

      ewise goodfevd(accept)(i,j)=ik=(irfsquared(i)*(a.^2))./(irfsquared(i)*ones),ik(j)

The irfsquared(i)*ones is the vector of total variances (sums across rows), so the ./ gives the fractions.

The code we wrote for doing multiple sign-restricted shocks uses:

Code: Select all

   impulse(noprint,model=varmodel,decomp=%identity(nvar),results=impulses,steps=nstep)

Because the shocks aren't orthogonal contemporaneously:

Code: Select all

   gset irfsquared 1 1     = %xt(impulses,t).^2
   gset irfsquared 2 nstep = irfsquared{1}+%xt(impulses,t).^2

[/code]

should have a big red X through it - the numbers produced are of almost no value; there are a whole bunch of non-zero covariances which this doesn't compute.

The solution to this is to return to the original base set of impulse responses, then back out the weights based upon those. That's really easy for the first one, since it just restores the original code:

Code: Select all

       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

The second (and any others) just requires pre-multiplying the non-orthogonal shock by inv(p).

Code: Select all

       @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 - Demand Shock
       *****************************************************
       do k=1,KMAX+1
          compute ik=%xt(impulses,k)*v2
          if ik(1)<0.or.ik(2)<0.or.ik(3)<0.or.ik(6)<0
             branch 105
       end do k

v1 and v2 are the two weight vectors on the "SWISH" factored draws. The FEVD's are computed with:

Code: Select all

       ewise goodfevd(accept)(i,j)=(ik=(irfsquared(i)*(v1.^2))./(irfsquared(i)*ones)),ik(j)
       ewise goodfevda(accept)(i,j)=(ik=(irfsquared(i)*(v2.^2))./(irfsquared(i)*ones)),ik(j)

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 12
 end(system)
 estimate(noprint,resid=resids)

 dec vect[strings] vl(6)
 compute vl=||'real GDP','GDP price defl','Comm. Price Ind.',$
      'Fed Funds Rate','Nonborr. Reserv.','Total Reserves'||

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

 * 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) goodfevd(1000)
 declare vect[rect] goodrespa(1000) goodfevda(1000)

 declare vector ik a(nvar) ones(nvar)
 declare series[rect] irfsquared
 compute ones=%const(1.0)
 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=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 - Demand Shock
       *****************************************************
       do k=1,KMAX+1
          compute ik=%xt(impulses,k)*v2
          if ik(1)<0.or.ik(2)<0.or.ik(3)<0.or.ik(6)<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)
       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)
       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='Subhea')
 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)
 smpl
 *
 clear upper lower resp
 *
 spgraph(vfields=3,hfields=2,hlabel='SECOND Impulse Responses with Pure-Sign Approach',subhea='Subhea')
 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='Impulse Responses 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 = 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(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='SECOND 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)
    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

[TL]

Dear Tom,

Thank you very much. I am now carefully reading your reply. I also have tried using the new code with my data set; it works.

Thank you very much for your help and kindness Tom.

Tim

[lhlee0506]

Dear Tom,

One problem which I am always confused about is what is the rule to set up sign restriction.
say, if we would like to disscuss positive government spending shock and it is put as the variable 1,
then we restric ik(1)<0 in the loop.
Is it correct?

[TomDoan]

The "IF" inside the loop lists the conditions which would cause you to reject a draw. So if you want only draws which have a positive reponse from variable 1, you would indeed test for ik(1)<0

[lhlee0506]

Dear Tom,

Thank you for your kind instruction.
I tried to apply your code into year delayed postivative shock for 4 periods, and it works.
I also tried to extend it to 3shocks case, but I think I could get wrong with something.
If I could have your guidance, it will ber very helpful.

@forcedfactor(force=column) sigmad i1~i2 f
compute step0=%xt(impulses,1)
compute step1=%xt(impulses,2)
compute step2=%xt(impulses,3)
compute step3=%xt(impulses,4)
dec rect v3(4,6)
compute %psubmat(v3,1,1,||0.0,1.0,0.0,0.0,0.0,0.0||*step0)
compute %psubmat(v3,2,1,||0.0,1.0,0.0,0.0,0.0,0.0||*step1)
compute %psubmat(v3,3,1,||0.0,1.0,0.0,0.0,0.0,0.0||*step2)
compute %psubmat(v3,4,1,||0.0,1.0,0.0,0.0,0.0,0.0||*step3)
compute rfperp=%perp(V3)
compute v3=rfperp*%ranmat(%cols(rfperp),1)
compute v3=v3/sqrt(%normsqr(v3))
compute v=i3=f*v3

[TomDoan]

What is it that you're trying to do? You've already created two shocks which must be satisfying sign restrictions. If you're trying to create a third shock which is orthogonal to the first two and satisfies your four equality constraints, it won't work. That's six constraints plus the length constraint. Even if you didn't have that problem, you need to do the equality constrained shock first since the calculation is based upon the shock being freely chosen from the original impulse vectors.

[lhlee0506]

Hi~

I tried to identify 3 shocks_Business cycle shock, monetary shock and government spending shock.
Government spending shock is taken as one year delayed shock.
I think I messed it up.
when I ran pure 3 shocks code modified from the previous post, I found one question that almost all error band of impulse response except those with sign restrictions include zero which imply "insignificant".
If I cut the sample periods, I have more variales which do not include zero, but I have the sawtooth shape of impulse responses.
Is there any tricks in dealing with data or sign restriction setting which I should notice?
Thanks for your great helps.

Code: Select all

open data data0508.xls
calendar 1970 1 4
data(format=xls,org=columns) 1970:01 2008:03 lrgspen1 lrgspen lrntax lrtax lrgexp lrgrev lrgcon lrginv1 lrginv $
lrgdp lrhgdp1 lrhgdp pgdpch cpich rint int lrpcon lrpinv1 lrpinv lrim lres lpgdp lm2 lpgdpch lrint lint vlpgdpch vlrint $
lrwmf lahmf lepmf depo ffr lrbc lrbcl lrtexp lrtexpca lrtexpcu lrtexpsp lrtrev lrtrevca lrtrevcu lrtrevsp
*
system(model=varmodel)
variables lrgdp lrpcon lrgspen lrntax lrpinv1 lint lm2 lpgdp
lags 1 to 4
end(system)
estimate(noprint)
*
dec vect[strings] vl(8)
compute vl=||'Real GDP','Private Consumption','Government Spending','Net Tax','Private Investment','Interest Rate','M2','GDP deflator'||
*
* 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 nkeep=1000
compute nvar=8
compute nstep=40
compute KMAX=5

* This is the standard setup for MC integration of an OLS VAR
*
dec symm s(8,8)
dec vect v1(8)      ;* For the unit vector on the 1st draw
dec vect v2(7)      ;* For the unit vector on the 2nd draw
dec vect v3(6)      ;* For the unit vector on the 3nd draw
dec vect v(8)     ;* 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(nkeep) goodfevd(nkeep)
declare vect[rect] goodrespa(nkeep) goodfevda(nkeep)
declare vect[rect] goodrespb(nkeep) goodfevdb(nkeep)

declare vector ik a(nvar) ones(nvar)
declare series[rect] irfsquared
compute ones=%const(1.0)
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(8),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=%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(1)<0 .or. ik(2)<0 .or. ik(5)<0  .or. ik(4)<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=%ranmat(8,1)            ;* Draw a random 8-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

      *****************************************************
      *  Second Set of Restrictions - monetary Shock
      *****************************************************
      do k=1,KMAX+1
         compute ik=%xt(impulses,k)*v
         if ik(6)<0 .or. ik(8)>0 .or.ik(7)>0
            branch 105
      end do k

      *
      * Meets the second restriction
      * Draw from the orthogonal complement of i2 (last four columns of the
      * factor "g").
      *
      @forcedfactor(force=column) sigmad i1~i2 f
      compute v3=%ranmat(8,1)            ;* Draw a random 8-vector
      compute v3(1)=0.0                  ;* Zero the 1st component
      compute v3(2)=0.0                  ;* Zero the 2nd component
      compute v3=v3/sqrt(%normsqr(v3))   ;* Normalize to unit length
      compute v=i3=f*v3                  ;* Transform to impulse vector

      *****************************************************
      *  Third Set of Restrictions - Govt Spending Shock
      *****************************************************
      do k=1,KMAX+1
         compute ik=%xt(impulses,k)*v
         if ik(3)<0
            branch 105
      end do k

      *
      * Meets Three restrictions
      *
      compute accept=accept+1
      dim goodrespb(accept)(nstep,nvar) goodfevdb(accept)(nstep,nvar)
      dim goodrespa(accept)(nstep,nvar) goodfevda(accept)(nstep,nvar)
      dim goodresp(accept)(nstep,nvar) goodfevd(accept)(nstep,nvar)
      ewise goodresp(accept)(i,j)=(ik=%xt(impulses,i)*i1),ik(j)
      ewise goodfevd(accept)(i,j)=ik=(irfsquared(i)*(i1.^2))./(irfsquared(i)*ones),ik(j)
      ewise goodrespa(accept)(i,j)=(ik=%xt(impulses,i)*i2),ik(j)
      ewise goodfevda(accept)(i,j)=ik=(irfsquared(i)*(i2.^2))./(irfsquared(i)*ones),ik(j)
      ewise goodrespb(accept)(i,j)=(ik=%xt(impulses,i)*i3),ik(j)
      ewise goodfevdb(accept)(i,j)=ik=(irfsquared(i)*(i3.^2))./(irfsquared(i)*ones),ik(j)
   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=4,hfields=2,hlabel='Impulse Responses with Pure-Sign Approach',subhea='Business Cycle Shock')
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)*100
      compute upper(k)=frac(2)*100
      compute resp(k)=%avg(work)*100
   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)

spgraph(vfields=4,hfields=2,hlabel='Impulse Responses with Pure-Sign Approach',subhea='Monetary Shock')
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)*100
      compute upper(k)=frac(2)*100
      compute resp(k)=%avg(work)*100
   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)

spgraph(vfields=4,hfields=2,hlabel='Impulse Responses with Pure-Sign Approach',subhea='Government Spending Shock')
do i=1,nvar
   compute minlower=maxupper=0.0
   smpl 1 accept
   do k=1,nstep
      set work = goodrespb(t)(k,i)
      compute frac=%fractiles(work,||.16,.84||)
      compute lower(k)=frac(1)*100
      compute upper(k)=frac(2)*100
      compute resp(k)=%avg(work)*100
   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)

spgraph(vfields=4,hfields=2,hlabel='Fraction of Variance Explained with Pure-Sign Approach',subhea='Business cycle shock')
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)*10000
      compute upper(k)=frac(3)*10000
      compute resp(k)=frac(2)*10000
   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)

spgraph(vfields=4,hfields=2,hlabel='Fraction of Variance Explained with Pure-Sign Approach',subhea='Monetary shock')
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)*10000
      compute upper(k)=frac(3)*10000
      compute resp(k)=frac(2)*10000
   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)

spgraph(vfields=4,hfields=2,hlabel='Fraction of Variance Explained with Pure-Sign Approach',subhea='Government Spending shock')
do i=1,nvar
   compute minlower=maxupper=0.0
   smpl 1 accept
   do k=1,nstep
      set work = goodfevdb(t)(k,i)
      compute frac=%fractiles(work,||.16,.50,.84||)
      compute lower(k)=frac(1)*10000
      compute upper(k)=frac(3)*10000
      compute resp(k)=frac(2)*10000
   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)



















































































[TomDoan]

Hi~

I tried to identify 3 shocks_Business cycle shock, monetary shock and government spending shock.
Government spending shock is taken as one year delayed shock.
I think I messed it up.
when I ran pure 3 shocks code modified from the previous post, I found one question that almost all error band of impulse response except those with sign restrictions include zero which imply "insignificant".
If I cut the sample periods, I have more variales which do not include zero, but I have the sawtooth shape of impulse responses.
Is there any tricks in dealing with data or sign restriction setting which I should notice?
Thanks for your great helps.

I'm not sure what you mean by "I found one question that almost all error band of impulse response except those with sign restrictions include zero". Are you saying that once you're past the end of the sign restrictions the responses no longer have the signs that you imposed on the shorter lags? That seems odd for a model with so many very persistent variables. You're not doing very many draws - you should probably see what happens when you increase that. Also, have you checked how easily the sub-draw procedure is getting a set of sign-restricted responses, that is, how many subdraws are typically needed in order to get a set which passes. If that's a high number, it may be an indication that the set of constraints is too rigid.

[lhlee0506]

Hi~

I am sorry for my imprecise expression.
I mean after meeting three shocks and producing all impulse responses of every shock,
their error band (16% and 84% percentils) will include zero, especially those variables without sign restriction.
for example, the impulse response graph of real GDP, private consumption, private investment...etc of government spending shock , thier error band will include zero.
According to Sims and Zha, error band includes zero seems to imply "stastiscally insignificant".
I worry my estimation finally produce almost insigificant results.
Perhaps I miss some points at dealing with data or setting sign restriction or something else.
By the way, how do I check if I have enough sub-draws to pass?
Thank you for your kindness.
I am really appreciated.

Best

[TomDoan]

Hi~

I am sorry for my imprecise expression.
I mean after meeting three shocks and producing all impulse responses of every shock,
their error band (16% and 84% percentils) will include zero, especially those variables without sign restriction.
for example, the impulse response graph of real GDP, private consumption, private investment...etc of government spending shock , thier error band will include zero.
According to Sims and Zha, error band includes zero seems to imply "stastiscally insignificant".
I worry my estimation finally produce almost insigificant results.
Perhaps I miss some points at dealing with data or setting sign restriction or something else.
By the way, how do I check if I have enough sub-draws to pass?
Thank you for your kindness.
I am really appreciated.

Best

You have to remember that your government spending shock is constructed to be orthogonal to your other two shocks, one of which was chosen to be positive on most of the real variables. It's certainly possible that isolated government spending shocks have no predictable effect on the other real variables. The lack of a "result" is a result.

Regarding the sub-draws, check your value of "draws" at the end of the loop. If that's fairly large (close to "n1"), you're not getting many good sub-draws.

[lhlee0506]

But why even though I identified only one government spending shock, the impulse responses of almost all variables also include zero?

[TomDoan]

But why even though I identified only one government spending shock, the impulse responses of almost all variables also include zero?

Do the following to verify that you're getting the programming correct.

1. Make sure you're getting enough acceptances. Do a DISP ACCEPT after the infobox(action=remove). If that isn't equal to your value of NKEEP, you need to increase the n1 and maybe n2 values.

2. If you're getting the right number of acceptances, do the following (again, after the infobox(action=remove))

Code: Select all

@forcedfactor(force=column) sigmad i1~i2~i3 f                                                     
impulse(factor=f,steps=20,model=varmodel,print)

See if the first three sets of responses satisfy the sign restrictions that you want. (The first block in the output are the responses to your first shock, the second block to your second shock, the third to your third. The last five are just fillers). They should. I don't see any programming errors, but it doesn't hurt to make sure.

If 1 is OK and 2 is OK, then your calculations are correct. If the responses of other variables to your government spending shock are insignificant, then, well, you have a topic for a paper.

[dennis0125hk]

Here's the code that I modified from this post. I also downloaded Forcefactor.src and installed it in WINRAT in order to run the program file.

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 v=i1=swish*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='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)
end(reset)


However, when I attempt to run it by WINRATS 7.1, below error message keeps appearing
## CP17. PROCEDURE/FUNCTION Must be Initial Statement in a Compiled Section
>>>>procedure <<<<
## SX11. Identifier F is Not Recognizable. Incorrect Option Field or Parameter Order?
>>>>%xsubmat(f,1,6,2,6)<<<<

How could I solved the problem?

Code: Select all

*
* ForcedFactor computes a factorization of an input NxN covariance matrix SIGMA
* which controls either
*
*   (a) a column or block of columns in the factorization itself
*   (b) a row or block of rows in the inverse of the factorization
*
* In the case (a), you provide an RxN matrix A. An RxR matrix PI is computed so
* that PI is upper triangular and A x PI makes up the first R columns in a matrix
* F that factors SIGMA (that is FF'=SIGMA). In particular, if R=1, then the first
* column of F is a scale multiple of A. Because of the structure of PI, if R>1,
* the first column of F is a scale multiple of the first column in A, the second
* column in F will be a linear combination of the first two columns of A, etc.
*
* In the case (b), you provide an NxR matrix A. An RxR matrix PI is computed so
* that PI is lower triangular and PI x A makes up the first R rows of the inverse
* of a matrix F that factors SIGMA. In particular, if R=1, the first row in F**-1
* is a scale multiple of A.
*
* This allows you to either force an orthogonalized component to hit the variables
* in a specific pattern (done by setting a column of the factorization), or to
* force that an orthogonalized component be formed from a particular linear
* combination of innovations (forces a row in the inverse). You set the one
* component, and ForcedFactor builds a factorization around it. With multiple
* input components, you can force the factorization to include linear
* combinations of your inputs.
*
* Syntax:
*
* @ForcedFactor( options ) sigma a f
*
*   Parameters:
*      sigma   SYMMETRIC covariance matrix (input)
*      a       RECTANGULAR (or VECTOR) - one dimension should be the same as sigma (input)
*      f       RECTANGULAR - output factor of sigma
*
*   Options:
*      force=[row]/column indicates whether a scale of "a" is to be a row in the
*        inverse or a column in the factorization itself.
*
* Examples:
*
* In a four variable system where the first two variables are two interest rates,
*
* @ForcedFactor sigma ||1.0,-1.0,0.0,0.0|| f1
* @ForcedFactor(force=column) sigma ||1.0,1.0,0.0,0.0|| f2
*
* f1 will be a factorization where the first orthogonal component is the
* (non-orthogonal) innovation in the difference between the rates. f2 will be a
* factorization where the first orthogonal component loads equally onto the two
* interest rates, and hits none of the other variables contemporaneously.
*
*  Revision Schedule
*   09/02 Written by Estima
*   11/05 Revised to allow multiple rows or columns
*
* %SVDZeros picks out the "r" smallest elements in the vector v, returning a
* VECT[INT] with the subscripts. This is used to locate the zero singular values
* in a singular value decomposition.
*
*
function %SVDZeros v r
type vect[int] %SVDZeros
type vector v
type integer r
local vector ranks
local integer i j
*
compute ranks=%ranks(v)
dim %SVDZeros(r)
compute j=0
do i=1,%rows(v)
   if ranks(i)<=r {
      compute j=j+1
      compute %SVDZeros(j)=i
   }
end do i
end
*
procedure forcedFactor sigma a f
type symm sigma
type rect a
type rect *f
*
option choice force 1 row column
*
local integer n r
local rect vv w pi p ff fill
local symm pp
local vect[rect] ss
local vect[int]  zeros
*
compute n=%rows(sigma)
compute r=%imin(%rows(a),%cols(a))
if %imax(%rows(a),%cols(a))<>n {
   disp "@ForcedFactor sigma a f"
   disp "A is dimension" %rows(a) "x" %cols(a)
   disp "SIGMA " n "x" n
   disp "One dimension of A must match SIGMA"
   return
}
if force==2 {
   compute p=inv(%decomp(sigma))
   if %rows(a)==n
      compute vv=p*a
   else
      compute vv=p*tr(a)
   compute w=%blockglue(||vv,%zeros(n,n-r)||)
   compute ss=%svdecomp(w)
   compute zeros=%SVDZeros(ss(2),n-r)
   compute pp=inv(tr(vv)*vv)
   compute pi=%psdfactor(pp,%seq(r,1))
   dim fill(n,n-r)
   ewise fill(i,j)=ss(1)(i,zeros(j))
   compute f=inv(p)*%blockglue(||vv*pi,fill||)
}
else {
   compute p=%decomp(sigma)
   if %cols(a)==n
      compute vv=a*p
   else
      compute vv=tr(a)*p
   compute w=%blockglue(||vv|%zeros(n-r,n)||)
   compute ss=%svdecomp(w)
   compute zeros=%SVDZeros(ss(2),n-r)
   compute pi=inv(%decomp(%outerxx(vv)))
   dim fill(n-r,n)
   ewise fill(i,j)=ss(3)(j,zeros(i))
   compute ff=%blockglue(||pi*vv|fill||)
   compute f=p*inv(ff)
}
end


[MC128]

Dear all,

I am thinking of including some exogenous variables (e.g. constant, dummies and more importantly, some stochastic variables) into Uhlig's sign restriction framwork. The problem that I have is how would the posterior distribution of Beta (the regression coefficients) and Sigma (the variance-covariance matrix) changed as a result. I know that Kadiyala and Karlsson (1997, p104) shows the general case of posterior distribution of normal-wishart priors in a model with exogenous variables but I am more concerned with the specific case of a weak (or flat) prior as in Uhlig (2005). So in my case, what should I input for the posterior distribution of Beta and Sigma?

Thanks for any comment!

MC128

[TomDoan]

Here's the code that I modified from this post. I also downloaded Forcefactor.src and installed it in WINRAT in order to run the program file.

<<snipped>>

However, when I attempt to run it by WINRATS 7.1, below error message keeps appearing
## CP17. PROCEDURE/FUNCTION Must be Initial Statement in a Compiled Section
>>>>procedure <<<<
## SX11. Identifier F is Not Recognizable. Incorrect Option Field or Parameter Order?
>>>>%xsubmat(f,1,6,2,6)<<<<

How could I solved the problem?

You missed a few things in copying the earlier code. First, you need to do a source forcedfactor.src instruction outside the loop. That's what's causing your immediate problem. You're also missing the declarations for v1 and v2. You need (again outside the loop), dec vect v1(6) v2(5). After correcting those, you'll find that you get a very low yield on acceptances, so you'll probably have to greatly increase n1 and n2.