Identifying VARs with sign restrictions

[TL]

Dear Tom,

I have been using the following code (two shocks) you kindly helped me.

In addition to impulse responses and variance decomposition, would it be possible to directly acquire the followings from this code?

1. The baseline series (in order to see the prediction errors (the deviation between the baseline and the actual series))

2. The contribution of each shock in explaining the deviation between the baseline and the actual series. That is, series of the baseline with the first shock and series of the baseline with the second shock.

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 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=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=6,hfields=2,subhea='Monetary Policy Shocks and Exchange Rate Shocks, Thailand')
    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
    *
    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]

This does the historical decomposition for the one shock model. It isn't that hard to modify it for two.

Code below corrected June 25, 2009

Code: Select all

*
* Based upon 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.
*
* This does an historical decomposition of the data showing the cumulative effects
* of the monetary policy shocks. Note that this is not part of the paper.
*
open data uhligdata.xls
calendar(m) 1965
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 nkeep=1000
compute nvar=6
compute nstep=60
compute KMAX=5

*************************************************************************************
compute hstart = 1998:1
compute hend   = 2002:12
compute nhor   = hend-hstart+1
dec vect[series] upaths(nvar) u(nvar) onesteps(nvar) base(nvar) baseplus(nvar)
clear upaths
*************************************************************************************
*
* 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
*
* Most draws are going to get rejected. We allow for up to nkeep 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 vector ik a ones(nvar)
declare series[rect] irfsquared
compute ones=%const(1.0)

*************************************************************************************
declare vector hk(nvar)
declare vect[rect] goodhist(nkeep)
*************************************************************************************

*
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 betau   =%ranmvkron(swish,sxx)
   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,errors=u,static)
   forecast(model=varmodel,from=hstart,to=hend,results=base)
*************************************************************************************
   *
   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
   *
   * Do the subdraws over the unit sphere. These give the weights on the
   * orthogonal components.
   *
   do subdraws=1,n2
      compute a=%ransphere(nvar)
      *
      * 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(4)<0.or.ik(3)>0.or.ik(2)>0.or.ik(5)>0
            branch 105
      end do k
      *
      * This is an accepted draw. Copy the information out. If we have enough
      * good ones, drop out of the overall loop.
      *
      compute accept=accept+1
      dim goodresp(accept)(nstep,nvar) goodfevd(accept)(nstep,nvar)
      ewise goodresp(accept)(i,j)=ik=%xt(impulses,i)*a,ik(j)
      ewise goodfevd(accept)(i,j)=ik=(irfsquared(i)*(a.^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)
      compute remap=swish*a*tr(a)*inv(swish)
      do t=hstart,hend
         compute %pt(upaths,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
      *
      * 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)
*************************************************************************************
      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=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)

spgraph(vfields=3,hfields=2,hlabel="Figure 8. Fraction of Variance Explained with Pure-Sign Approach")
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(ticks,number=0,min=0.0,picture="##.##",header="Fraction Explained for "+vl(i)) 3
   # resp
   # upper / 2
   # lower / 2
end do i
*
spgraph(done)

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

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(ticks,dates,picture="##.##",header="Contribution of Monetary Policy shock to "+vl(i)) 3
   # respH hstart hend
   # upperH / 2
   # lowerH / 2

end do i
*
spgraph(done)

[TL]

Dear Tom,

I have modified the code for two shocks. Nevertheless, results do not appear. I might do something wrong. Could you please give me some suggestions?

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 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=60
    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)
    clear upaths
    *************************************************************************************
    *
    * 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

    *************************************************************************************
    declare vector hk(nvar)
    declare vect[rect] goodhist(1000) goodhista(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 - 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)
     *************************************************************************************
          *
          * 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)
          do t=hstart,hend
             compute %pt(upaths,t,p*(a.*(inv(p)*%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
          *
          * 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)=baseplus(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=2,subhea='Monetary Policy Shocks and Exchange Rate Shocks, Thailand')
    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
    *
    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
    *************************************************************************************

    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(ticks,dates,picture="##.##",header="Contribution of Monetary Policy 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(ticks,dates,picture="##.##",header="Contribution of Exchange Rate shocks to "+vl(i)) 3
       # respH hstart hend
       # upperH / 2
       # lowerH / 2

    end do i
    *
    spgraph(done)

[TomDoan]

Code: Select all

          do t=hstart,hend
             compute %pt(upaths,t,p*(a.*(inv(p)*%xt(u,t))))
          end do t


With your coding, this isn't <<a>> in the middle line, it's <<i1>> for the first and <<i2>> for the second. You need to repeat the calculation for each of those (this + the forecast + ewise) for each of the two shocks.

[TL]

Dear Tom,

Thank you very much for your suggestions. I believe that I have successfully modified the code for two shocks. Could you please give me some 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 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=60
    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) ua(nvar) onestepsa(nvar) basea(nvar) baseplusa(nvar)
    clear upaths
    clear upathsa
    *************************************************************************************
    *
    * 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

    *************************************************************************************
    declare vector hk(nvar)
    declare vect[rect] goodhist(1000) goodhista(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)
        forecast(model=varmodel,from=hstart,to=hend,results=onestepsa,errors=ua,static)
        forecast(model=varmodel,from=hstart,to=hend,results=basea)
     *************************************************************************************
        *
        * 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)
     *************************************************************************************
          *
          * 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)
          do t=hstart,hend
             compute %pt(upaths,t,p*(i1.*(inv(p)*%xt(u,t))))
          end do t
          do t=hstart,hend
             compute %pt(upathsa,t,p*(i2.*(inv(p)*%xt(ua,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
          *
          * 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)-basea(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=2,subhea='Monetary Policy Shocks and Exchange Rate Shocks, Thailand')
    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
    *
    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
    *************************************************************************************

    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 Monetary Policy 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 Exchange Rate shocks to "+vl(i)) 3
       # respH hstart hend
       # upperH / 2
       # lowerH / 2

    end do i
    *
    spgraph(done)

[TomDoan]

Dear Tom,

Thank you very much for your suggestions. I believe that I have successfully modified the code for two shocks. Could you please give me some feedbacks whether my code is correct?

Thank you very much for your help and kindness.

Tim

Code: Select all

        forecast(model=varmodel,from=hstart,to=hend,results=onesteps,errors=u,static)
        forecast(model=varmodel,from=hstart,to=hend,results=base)
        forecast(model=varmodel,from=hstart,to=hend,results=onestepsa,errors=ua,static)
        forecast(model=varmodel,from=hstart,to=hend,results=basea)

This is overkill. Your ua and basea will be exactly the same as u and base.

Code: Select all

          do t=hstart,hend
             compute %pt(upaths,t,p*(i1.*(inv(p)*%xt(u,t))))
          end do t
          do t=hstart,hend
             compute %pt(upathsa,t,p*(i2.*(inv(p)*%xt(ua,t))))
          end do t


The correct calculation for this is described in the "Historical Decomposition (Technical Paper)" topic. First, with your definitions of the vectors, you want to be using the v1 and v2, not the i1 and i2. (v1 and v2 are your norm one weights on the Choleski factor). The following should replace yours above:

Code: Select all

          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


Code: Select all

          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)-basea(j)(i+hstart-1)

As pointed out above, there's no need to compute a

Code: Select all

base

and a

Code: Select all

basea

. Just use base on the right side of both.

[MC128]

Dear Tom,

I am trying to do the same analysis as Mountford and Uhlig (2005) "What are the effects of fiscal policy shocks?"

In particular, the following three shocks are identified:

1. business cycle shock, which causes output, revenue, consumption and investment to increase
2. monetary policy shock, which causes interest rate to rise, and prices and reserves to fall. This shock is also assumed to be orthogonal to business cycle shock
3. government spending shock, which cause government spending to increase. This shock is assumed to be orthogonal to the first two shocks.

All shocks were identified using penalty function approach.

But one difference from Mountford and Uhlig is that I tried to make one variable, the producer price of crude material exogenous to the rest. In doing so, I follow the analysis in Zha (1999) "Block recursion and Structural Vector Autoregressions", and use the posterior distribution of beta and sigma derived under an improper prior on sigma.

I have used the codes in this forum and take into of the above in my analysis. But I am very worried that there may be some mistakes. Please could you give me some comments on my code and approach?

The code and data have been attached for your reference. Thank you very much.....

MC

[TomDoan]

Don't use the same "G" for each of the minimizations. The first one uses a NVAR-1 vector, the second an NVAR-2 and the third an NVAR-3. For instance, for your second you would want:

Code: Select all

dec vect g2(nvar-2)

at the top, then

Code: Select all

find(noprint) min func
      compute v2=%stereo(g2)
      compute v=i2=f*v2
      ...


[TL]

Dear Tom,

Following your suggestions, I have corrected the code.

I hope the following code will prove useful to other users working on impulse responses, variance decomposition and historical contribution using sign restrictions for two shocks.

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 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=60
    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)
    clear upaths
    clear upathsa
    *************************************************************************************
    *
    * 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

    *************************************************************************************
    declare vector hk(nvar)
    declare vect[rect] goodhist(1000) goodhista(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 - 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)
     *************************************************************************************
          *
          * 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)
          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
          *
          * 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
          *
          * 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)
     *************************************************************************************
           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=2,subhea='Monetary Policy Shocks and Exchange Rate Shocks, Thailand')
    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
    *
    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
    *************************************************************************************

    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 Monetary Policy 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 Exchange Rate shocks to "+vl(i)) 3
       # respH hstart hend
       # upperH / 2
       # lowerH / 2

    end do i
    *
    spgraph(done)

[MC128]

Dear Tom,

Thank you so much for your reply. I have made some changes following your suggestions.. but one problem is that the resulting impulse responses no longer satisfy the sign-restriction. For instance, following the monetary policy shock, interest rate and prices move in the wrong direction........do you know why this happened?

Many thanks!

MC

[TomDoan]

Dear Tom,

Thank you so much for your reply. I have made some changes following your suggestions.. but one problem is that the resulting impulse responses no longer satisfy the sign-restriction. For instance, following the monetary policy shock, interest rate and prices move in the wrong direction........do you know why this happened?

Many thanks!

MC

First, you're not going to get the correct FEVD information when you start with the IMPULSE instruction with identity matrix shocks. You need to start with orthogonalized shocks in order to use the fairly simple calculations shown. If you switch to Choleski factor shocks, you'll have to make adjustments to everything else since the weights on the impulse responses will be the draws for the unit vector, rather than the factor times the draw.

You're having trouble because you're not switching non-linear parameter sets inside the loop. You need NONLIN G1 before solving the first subproblem, NONLIN G2 before the second and NONLIN G3 before the third.

[MC128]

Dear Tom,

Thank you so much for your helpful comment! The problem of having wrong sign is now solved....

I am now trying to come up with two set of codes for replicating the two analysis in Mountford and Uhlig (2008):

1. identification of business cycle, monetary and government spending and revenue shocks.
2. implementation lag of 4 quarters for government spending shock, which will then remain positive for 4 quarters

I find that:

For the first set of code: the shape of the impulse responses are very similar to that in Mountford and Uhlig (2008), although there are some minor quantitative differences.
For the second set of code: The shape of the impulse responses are a bit different

May I ask whether there are any mistake, particularly in the second set of codes?

Many thanks for your kind attention!

MC

[TomDoan]

Dear Tom,

Thank you so much for your helpful comment! The problem of having wrong sign is now solved....

I am now trying to come up with two set of codes for replicating the two analysis in Mountford and Uhlig (2008):

1. identification of business cycle, monetary and government spending and revenue shocks.
2. implementation lag of 4 quarters for government spending shock, which will then remain positive for 4 quarters

I find that:

For the first set of code: the shape of the impulse responses are very similar to that in Mountford and Uhlig (2008), although there are some minor quantitative differences.
For the second set of code: The shape of the impulse responses are a bit different

May I ask whether there are any mistake, particularly in the second set of codes?

Many thanks for your kind attention!

MC

Did you get the data set from the authors, or is it just your attempt to reconstruct their data?

[MC128]

Dear Tom,

Thanks for your reply. The dataset is downloaded from the website of journal of applied econometrics, where the paper of mountford and ublig is published.

Many thanks for your kind attention!

MC

[TomDoan]

I reorganized this to use a function to calculate the penalty function, which shortens the program considerably. (Also makes it much easier to read). This does everything in the paper that uses only the impact shocks. You need the MCGraphIRF procedure from the Procedures forum as well.

The zero restrictions are done with (for pegging variable 2):

dec rect r(KMAX+1,nvar)
ewise r(i,j)=ix=%xt(impulses,i),ix(2,j)
@forcedfactor(force=column) sigmad i1~i2~p*tr(r) f

The "R" matrix is the same as the one in the paper. The p*tr(r) input in the forcedfactor is to get it in the same form as the orthogonality conditions on the new shocks that forcedfactor is imposing. (forcedfactor is doing basically a "Gram-Schmidt" orthogonalization on a rotated basis).

Code: Select all

*
* Replication file for Mountford and Uhlig (2009), "What are the Effects of Fiscal
* Policy Shocks?", Journal of Applied Econometrics.
*
open data mudata.xls
calendar(q) 1955:1
data(format=xls,org=columns) 1955:1 2000:4 rgdp rbpexp rbprev ffrt ares ppic gdpdef rcon rnresinv wage
*
set rgdp     = rgdp*100
set rbpexp   = rbpexp*100
set rbprev   = rbprev*100
set ffrt     = ffrt
set ares     = ares*100
set ppic     = ppic*100
set gdpdef   = gdpdef*100
set rcon     = rcon*100
set rnresinv = rnresinv*100
set wage     = wage*100
*
system(model=varmodel)
variables rgdp rbpexp rbprev ffrt ares ppic gdpdef rcon rnresinv wage
lags 1 to 6
end(system)
estimate(noprint)
*
dec vect[strings] vl(10)
compute vl=||"GDP","Expenditure","Revenue","Federal Rate","Reserves","PPI","GDP deflator","Consumption","Investment","Wages"||
*
* Compute scale factors (standard deviations of first differences) for weighting
* the penalty function
*
dec vect scales(10)
dec real func
compute fill=0
dofor i = rgdp rbpexp rbprev ffrt ares ppic gdpdef rcon rnresinv wage
   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=250
compute nvar  =10
compute nstep =25
compute KMAX  =3
dec rect[series] impulses
*
* Parameter vectors for the nested optimization problems. They should have a size
* equal to nvar less 1 + the number of constraints imposed prior to their estimation.
* Those prior constraints can be orthogonality to previous shocks, or can be zero
* restrictions.
*
compute nshocks=6
dec vect[strings] sl(nshocks)
compute sl=||"Business Cycle","Monetary Policy","Revenue","Revenue (delayed)","Spending","Spending (delayed)"||
*
compute [vect] g1 =%fill(nvar-1,1,1.0)      ;* Parameters for the first identified shock
compute [vect] g2 =%fill(nvar-2,1,1.0)      ;* Parameters for the second identified shock
compute [vect] g3c=%fill(nvar-3,1,1.0)      ;* Parameters for the third identified shock (revenue)
compute [vect] g3d=%fill(nvar-7,1,1.0)      ;* Parameters for the third identified shock (revenue, delayed)
compute [vect] g3e=%fill(nvar-3,1,1.0)      ;* Parameters for the third identified shock (spending)
compute [vect] g3f=%fill(nvar-7,1,1.0)      ;* Parameters for the third identified shock (spending, delayed)
*******************************************************************************
*
* UhligPenalty(q,first,last,constrained) returns the penalty function for the
* weight vector <<q>> (applied to the impulse responses in the global VECT[SERIES]
* impulses). <<first>> and <<last>> are the range of responses constrained (1 =
* impact response). <<constrained>> is a VECT[INTEGER] with the variable positions
* being constrained. Use +slot for a positivity constraint and -slot for a
* negativity constraint. That is, if you want the first 4 responses to be positive
* on the 2nd variable and negative on the third, you would use the function
*
* UhligPenalty(q,1,4,||2,-3||)
*
function UhligPenalty q first last constrained
type real      UhligPenalty
type vector    q
type integer   first last
type vect[int] constrained
*
local integer i k
local real    func
local vector  ik
*
compute func=0.0
do k=first,last
   compute ik=(%xt(impulses,k)*q)./scales
   do i=1,%rows(constrained)
      if constrained(i)<0
         compute value= ik(-constrained(i))
      else
         compute value=-ik(constrained(i))
      compute func=func+%if(value<0.0,value,100*value)
   end do i
end do k
compute UhligPenalty=func
end
*******************************************************************************
*
* 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 vect[rect] %%responses
dim %%responses(ndraws)
*
declare vector ik a(nvar) ones(nvar)
declare series[rect] irfsquared
compute ones=%const(1.0)
source forcedfactor.src
*
infobox(action=define,progress,lower=1,upper=ndraws) 'Monte Carlo Integration'
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 p       =%decomp(sigmad)
   compute betadraw=betaols+%ranmvkron(p,sxx)
   compute %modelsetcoeffs(varmodel,betadraw)
   *
   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
   *
   * Minimize the penalty function, starting from the last set of minimizers
   *
   ************************************************
   * First in Order - Business Cycle Vector
   ************************************************
   nonlin g1
   find(noprint) min func
      compute v1=%stereo(g1)
      compute func=UhligPenalty(v1,1,KMAX+1,||+1,+3,+8,+9||)
   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 v1=%stereo(g1)
      compute func=UhligPenalty(v1,1,KMAX+1,||+1,+3,+8,+9||)
   end find
   *
   * If the two estimates don't match, reject the draw.
   *
   if abs(testfunc-func)>.01
      branch newdraw

   compute i1=p*v1
   *
   * Compute a factor with i1 forced as the first column.
   * Transform the orthogonal complement (remaining columns) back into weights on
   * the Choleski factor.
   *
   @forcedfactor(force=column) sigmad i1 f
   compute r2=inv(p)*%xsubmat(f,1,nvar,2,nvar)

   *****************************************************
   *  Second in Order - Monetary Shock
   *****************************************************
   *
   nonlin g2
   find(noprint) min func
      compute v2=r2*%stereo(g2)
      compute func=UhligPenalty(v2,1,KMAX+1,||+4,-5,-6,-7||)
   end find
   compute testfunc=func,testbeta=%beta
   *
   * Try the minimization again, starting from a standard set of values
   *
   compute g2=%fill(nvar-2,1,1.0)
   find(noprint) min func
      compute v2=r2*%stereo(g2)
      compute func=UhligPenalty(v2,1,KMAX+1,||+4,-5,-6,-7||)
   end find
   *
   * If the two estimates don't match, reject the draw.
   *
   if abs(testfunc-func)>.01
      branch newdraw

   compute i2=p*v2
   *
   * Compute a factor with i1~i2 forced as the first two columns.
   * Transform the orthogonal complement (remaining columns) back into weights on
   * the Choleski factor.
   *
   @forcedfactor(force=column) sigmad i1~i2 f
   compute r3=inv(p)*%xsubmat(f,1,nvar,3,nvar)

   *****************************************************
   *  Third in Order - Govt Revenue Shock (no delay)
   *****************************************************
   *
   * Draw from the orthogonal complement of i1 and i2 (final nvar-2 columns
   * in the factor).
   *
   nonlin g3c
   find(noprint) min func
      compute v3c=r3*%stereo(g3c)
      compute func=UhligPenalty(v3c,1,KMAX+1,||+3||)
   end find
   compute testfunc=func,testbeta=%beta
   *
   * Try the minimization again, starting from a standard set of values
   *
   compute g3c=%fill(nvar-3,1,1.0)
   find(noprint) min func
      compute v3c=r3*%stereo(g3c)
      compute func=UhligPenalty(v3c,1,KMAX+1,||+3||)
   end find
   *
   * If the two estimates don't match, reject the draw.
   *
   if abs(testfunc-func)>.01
      branch newdraw

   compute i3c=p*v3c

   *****************************************************
   *  Third in Order - Govt Revenue Shock (one year delay)
   *****************************************************

   dec rect r(KMAX+1,nvar)
   ewise r(i,j)=ix=%xt(impulses,i),ix(3,j)
   @forcedfactor(force=column) sigmad i1~i2~p*tr(r) f
   compute r3x=inv(p)*%xsubmat(f,1,nvar,7,nvar)

   nonlin g3d
   find(noprint) min func
      compute v3d=r3x*%stereo(g3d)
      compute func=UhligPenalty(v3d,KMAX+2,KMAX+5,||+3||)
   end find
   compute testfunc=func,testbeta=%beta
   *
   * Try the minimization again, starting from a standard set of values
   *
   compute g3d=%fill(nvar-7,1,1.0)
   find(noprint) min func
      compute v3d=r3x*%stereo(g3d)
      compute func=UhligPenalty(v3d,KMAX+2,KMAX+5,||+3||)
   end find
   *
   * If the two estimates don't match, reject the draw.
   *
   if abs(testfunc-func)>.01
      branch newdraw

   compute i3d=p*v3d

   *****************************************************
   *  Third in Order - Govt Spending Shock (no delay)
   *****************************************************

   nonlin g3e
   find(noprint) min func
      compute v3e=r3*%stereo(g3e)
      compute func=UhligPenalty(v3e,1,KMAX+1,||+2||)
   end find
   compute testfunc=func,testbeta=%beta
   *
   * Try the minimization again, starting from a standard set of values
   *
   compute g3e=%fill(nvar-3,1,1.0)
   find(noprint) min func
      compute v3e=r3*%stereo(g3e)
      compute func=UhligPenalty(v3e,1,KMAX+1,||+2||)
   end find
   *
   * If the two estimates don't match, reject the draw.
   *
   if abs(testfunc-func)>.01
      branch newdraw

   compute i3e=p*v3e

   *****************************************************
   *  Third in Order - Govt Revenue Shock (one year delay)
   *****************************************************

   dec rect r(KMAX+1,nvar)
   ewise r(i,j)=ix=%xt(impulses,i),ix(2,j)
   @forcedfactor(force=column) sigmad i1~i2~p*tr(r) f
   compute r3x=inv(p)*%xsubmat(f,1,nvar,7,nvar)

   nonlin g3f
   find(noprint) min func
      compute v3f=r3x*%stereo(g3f)
      compute func=UhligPenalty(v3f,KMAX+2,KMAX+5,||+2||)
   end find
   compute testfunc=func,testbeta=%beta
   *
   * Try the minimization again, starting from a standard set of values
   *
   compute g3f=%fill(nvar-7,1,1.0)
   find(noprint) min func
      compute v3f=r3x*%stereo(g3f)
      compute func=UhligPenalty(v3f,KMAX+2,KMAX+5,||+2||)
   end find
   *
   * If the two estimates don't match, reject the draw.
   *
   if abs(testfunc-func)>.01
      branch newdraw

   compute i3f=p*v3f

   *
   * Meets all restrictions
   *
   compute accept=accept+1
   *
   dim %%responses(accept)(nshocks*nvar,nstep)
   compute vweights=v1~v2~v3c~v3d~v3e~v3f
   ewise %%responses(accept)(i,j)=(ik=%vec(%xt(impulses,j)*vweights)),ik(i)

*   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)
*   ewise goodfevdb(accept)(i,j)=ik=(irfsquared(i)*(v3.^2))./(irfsquared(i)*ones),ik(j)
*   ewise goodfevdc(accept)(i,j)=ik=(irfsquared(i)*(v4.^2))./(irfsquared(i)*ones),ik(j)

   if accept>=ndraws
      break
   infobox(current=accept)
   :newdraw
end do draws
infobox(action=remove)
disp accept
*
@MCGraphIRF(model=varmodel,varlabels=vl,shocklabels=sl,$
  page=byshock,columns=2,include=||1,2,10,4,6,8,3,9,5,7||)