Identifying VARs with sign restrictions

[fabio fornari]

thanks a lot Tom.
will see how it works with three shocks before moving on ... needs really lots of draws to get acceptances ...
will let you know
thanks f

[mmk]

Hi

Did anyone try to extend Uhilg_3 file for the paper published (2005) with penalty function to more than 6 variables? can anyone tell me how the penalty is going to change in this case?

Thanks for help

[TomDoan]

Hi

Did anyone try to extend Uhilg_3 file for the paper published (2005) with penalty function to more than 6 variables? can anyone tell me how the penalty is going to change in this case?

Thanks for help

This uses an alternative and more flexible parameterization of the unit sphere than was used by Uhlig. (Note that it's really the parameterization that needs to change; the calculation of the penalty function is similar regardless). I didn't use it in the original uhlig3.prg because the %stereo function was added with version 6.35 of RATS. There are quite a few ways to parameterize the unit sphere. The stereoscopic projection used in the %stereo function takes the N-1 vector x and maps it to the N-vector formed by

y(i)=2x(i)/(||x||^2+1)
y(N)=(||x||^2-1)/(||x||^2+1)

All such mappings have at least one minor flaw. In this case, the only way to make y(N)=0 is to make the x's very large. (Uhlig's mapping has the problem of non-uniqueness, as well as being rather difficult to generalize easily).

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. Penalty function approach.
*
open data uhligdata.xls
calendar(m) 1965
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)
*
dec vect[strings] vl(6)
compute vl=||"real GDP","GDP price defl","Comm. Price Ind.",$
             "Fed Funds Rate","Nonborr. Reserv.","Total Reserves"||
*
* Get the scale factors to be used for impulse responses
*
dec vect scales(6)
dec real func
compute fill=0
dofor i = gdpc1 gdpdef cprindex fedfunds bognonbr totresns
   diff i / diff1
   stats(noprint) diff1
   compute fill=fill+1,scales(fill)=sqrt(%variance)
end dofor i
*
* ndraws is the desired number of draws from the posterior of 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 ndraws=1000
compute nvar=6
compute nstep=60
compute KMAX=5
declare vect[rect] goodresp(ndraws)
declare vector ik a(nvar)
*
dec vect g(nvar-1)
compute g=%fill(nvar-1,1,1.0)
nonlin g
*
* 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)
*
infobox(action=define,progress,lower=1,upper=ndraws) "Monte Carlo Integration"
*
* Plan for double the number of draws to allow for the ones which we
* reject.
*
compute accept=0
do draws=1,ndraws*2
   *
   * 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)
   *
   * Minimize the penalty function, starting from the last set of minimizers
   *
   find(noprint) min func
      compute a=%stereo(g)
      compute func=0.0
      do k=1,KMAX+1
         compute ik=(%xt(impulses,k)*a)./scales
         compute func=func+%if(ik(4)>0,-ik(4),-100*ik(4))+$
                           %if(ik(3)<0,ik(3),100*ik(3))+$
                           %if(ik(2)<0,ik(2),100*ik(2))+$
                           %if(ik(5)<0,ik(5),100*ik(5))
      end do k
   end find
   compute testfunc=func,testbeta=%beta
   *
   * Try the minimization again, starting from a standard set of values
   *
   compute g=%fill(nvar-1,1,1.0)
   find(noprint) min func
      compute a=%stereo(g)
      compute func=0.0
      do k=1,KMAX+1
         compute ik=(%xt(impulses,k)*a)./scales
         compute func=func+%if(ik(4)>0,-ik(4),-100*ik(4))+$
                           %if(ik(3)<0,ik(3),100*ik(3))+$
                           %if(ik(2)<0,ik(2),100*ik(2))+$
                           %if(ik(5)<0,ik(5),100*ik(5))
      end do k
   end find
   *
   * If the two estimates don't match, reject the draw. If they do, copy out the
   * impulse responses.
   *
   if abs(testfunc-func)<.01 {
      compute accept=accept+1
      dim goodresp(accept)(nstep,nvar)
      ewise goodresp(accept)(i,j)=ik=%xt(impulses,i)*a,ik(j)
   }
   *
   * If we've hit out desired number of accepted draws, break
   *
   if accept>=ndraws
      break
   infobox(current=accept)
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 14. Impulse Responses with Penalty Function 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)

[lhlee0506]

Hi~Tom,

Would please you kindly instruct me how to collect the value of the initial demand/monetary shock and how to show the variance of decomposition?

[lhlee0506]

Hi~

I tried to put the variance decomposition in the previous code with 2 shocks you posted few weeks ago.
and I got message as following,

## SX22. Expected Type REAL, Got MATRIX(REAL) Instead
>>>>%xt(impulses,t).^2<<<<

Would you kindly instruct me how to fix it?

Thank you very much!

[TomDoan]

You need a

Code: Select all

DEC SERIES[RECT] irfsquared

instruction. Otherwise GSET doesn't know what type it's expecting and defaults to just a series of reals.

[lhlee0506]

Hello,

I tried to follow up your instruction in the latest RATS letter which demonstrateds the impact responses on variables 1 and 2 to be equal. But I failed. Perhaps I misunderstood your instruction.
If you could give me some guidances, I will be very appreciated.

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 levelsa1.xls
calendar 1970 1 4
data(format=xls,org=columns) 1970:01 2008:03 lrgspen1 lrntax lrgdp pgdpch rint lrpcon lrpinv
*

system(model=varmodel)
variables lrgspen1 lrntax lrgdp pgdpch rint lrpcon
lags 1 to 4
end(system)
estimate(noprint)
*
dec vect[strings] vl(6)
compute vl=||'Govt spending','net tax','real gdp','GDP deflator','real interest rate','Private onsumption'||
*
* 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
      compute rperp=$%perp(||1.0,-1.0,0.0,0.0,0.0,0.0||*swish)
      compute a=rperp*%ranmat(5,1),a=a/sqrt(%normsqr(a))
      *
      * Check that the responses have the correct signs for steps 1 to KMAX+1
      * (+1 because in the paper, the steps used 0-based subscripts, rather than
      * the 1 based subscripts used by RATS).
      *
      do k=1,KMAX+1
         compute ik=%xt(impulses,k)*a
         if ik(1)<0
         branch 105
      end do k
      *
      * This is an accepted draw. Copy the information out. If we have our 1000
      * good ones, drop out of the overall loop.
      *
      compute accept=accept+1
      dim goodresp(accept)(nstep,nvar)
      ewise goodresp(accept)(i,j)=ik=%xt(impulses,i)*a,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',subhea='monetary 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)
      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]

Hello,

I tried to follow up your instruction in the latest RATS letter which demonstrateds the impact responses on variables 1 and 2 to be equal. But I failed. Perhaps I misunderstood your instruction.
If you could give me some guidances, I will be very appreciated.

Code: Select all

      compute rperp=$%perp(||1.0,-1.0,0.0,0.0,0.0,0.0||*swish)
 

The line I isolated has been garbled. (Probably didn't fit in one column in the newsletter). Get rid of the $ and it should be fine.

[lhlee0506]

Hi~

Thank you for your instruction.
After removing $, the code seems to wotk.
But I am still confused about how to set up the sign restriction after putting variable1,2 to be equal.
This modification seems only make variable 1,2 to be equal for one period, how to make it last for four periods?
Hope to have your kind responses.
Thanks a lot!

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 levelsa1.xls
calendar 1970 1 4
data(format=xls,org=columns) 1970:01 2008:03 lrgspen1 lrntax lrgdp pgdpch rint lrpcon lrpinv
*

system(model=varmodel)
variables lrgspen1 lrntax lrgdp pgdpch rint lrpcon
lags 1 to 4
end(system)
estimate(noprint)
*
dec vect[strings] vl(6)
compute vl=||'Govt spending','net tax','real gdp','GDP deflator','real interest rate','Private onsumption'||
*
* 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
      compute rperp=%perp(||1.0,-1.0,0.0,0.0,0.0,0.0||*swish)
      compute a=rperp*%ranmat(5,1),a=a/sqrt(%normsqr(a))
      *
      * Check that the responses have the correct signs for steps 1 to KMAX+1
      * (+1 because in the paper, the steps used 0-based subscripts, rather than
      * the 1 based subscripts used by RATS).
      *
      do k=1,KMAX+1
         compute ik=%xt(impulses,k)*a
         if ik(1)<0
         branch 105
      end do k
      *
      * This is an accepted draw. Copy the information out. If we have our 1000
      * good ones, drop out of the overall loop.
      *
      compute accept=accept+1
      dim goodresp(accept)(nstep,nvar)
      ewise goodresp(accept)(i,j)=ik=%xt(impulses,i)*a,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',subhea='monetary 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)
      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)

[TL]

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

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=4000
    compute n2=4000
    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=%identity(6),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(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=%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.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)*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)
          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)
          display resp(k)
       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)
          display resp(k)
       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)
          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='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)
          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

[TomDoan]

Hi~

Thank you for your instruction.
After removing $, the code seems to wotk.
But I am still confused about how to set up the sign restriction after putting variable1,2 to be equal.
This modification seems only make variable 1,2 to be equal for one period, how to make it last for four periods?
Hope to have your kind responses.
Thanks a lot!

This (the alpha weights on the "swish" matrix) will create responses for the first two variables which are equal to each other for the first four periods. Note that this severely restricts the possible shocks, so it may be difficult for the draws for shocks to satisfy any other constraints.

Code: Select all

   impulse(noprint,model=varmodel,decomp=swish,results=impulses,steps=nstep)
   compute step0=%xt(impulses,1)
   compute step1=%xt(impulses,2)
   compute step2=%xt(impulses,3)
   compute step3=%xt(impulses,4)
   dec rect r(4,6)
   compute %psubmat(r,1,1,||1.0,-1.0,0.0,0.0,0.0,0.0||*step0)
   compute %psubmat(r,2,1,||1.0,-1.0,0.0,0.0,0.0,0.0||*step1)
   compute %psubmat(r,3,1,||1.0,-1.0,0.0,0.0,0.0,0.0||*step2)
   compute %psubmat(r,4,1,||1.0,-1.0,0.0,0.0,0.0,0.0||*step3)
   compute rfperp=%perp(r)
   compute alpha=rfperp*%ranmat(%cols(rfperp),1)
   compute alpha=alpha/sqrt(%normsqr(alpha))

[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

[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