Identifying VARs with sign restrictions

[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 Cholesky 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 Cholesky 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||)

[MC128]

Dear Tom,

Thank you so much for correcting my codes!....I am now trying to adjust the code to calculate the impulse responses of the model subject to a restricted 1% shock to revenue. I then move on to calculate the cumulative response and the multiplier....the shape of the responses seem very similar to those in Mountford and Uhlig (2008). but somehow the mean response of multiplier is a bit weird....Is this normal? (for your information, I follow closely the GUASS codes of Mountford and Uhlig (2008), available at the website of journal of applied econometrics)

Again, many thanks!

MC

[TomDoan]

There's a discrepancy between the results I get with RATS and the ones in the paper regarding one of the six impulse vectors (the delayed revenue shock). We're trying to get that sorted out.

[TL]

Dear Tom,

For an impulse response, I would like to normalize one variable (interest rate) to increase by 1 percentage point initially, as the effect of a contractionary monetary policy shock.

Could you please give me suggestions on how to perform such normalization?

Thank you very much for your help and kindness.

Sincerely,

Tim

[TomDoan]

Just scale the shock - if x is your impulse vector (in data space, not orthogonalized response space) and you want it to be scaled to hit component i with a value 1, just replace it with x/x(i).

[TL]

Dear Tom,

Thank you very much for your reply. However, I still have one more question I would like to ask you. Sorry for not making the question clear previously.

With the following code, I would like to normalize the third variable (INT) to increase by 1 initially, as the effect of the first shock.

For the second shock, I do not want to normalize anything.

Could you please give me suggestions on how to perform this?

Thank you very much for your help and kindness.

Sincerely,

Tim

Code: Select all

     OPEN DATA All_4_Thailand.xls
    CALENDAR 1989 2 12
    compute missc=1.0e+32
    data(format=xlS,org=columns) 1989:2 2008:11 CPI Y INT SS FF PP

    set CPI = CPI*100.0
    set Y   = Y*100.0
    set SS  = SS*100.0
    set FF  = FF*100.0
    set PP  = PP*100.0
    *
    system(model=varmodel)
    variables Y CPI INT PP FF SS
    lags 1 to 6
    deterministic constant
    end(system)
    estimate(noprint,resid=resids)

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

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

    *************************************************************************************
    compute hstart = 1995:1
    compute hend   = 1999: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(2)>0.or.ik(1)>0.or.ik(3)<0.or.ik(4)>0
                  branch 105
           end do k
           *
           * Meets the first restriction
           * Draw from the orthogonal complement of i1 (last five columns of the
           * factor "f").
           *
           @forcedfactor(force=column) sigmad i1 f
           compute v2=%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(2)<0.or.ik(1)<0.or.ik(3)<0.or.ik(4)<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)
        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
    *************************************************************************************

    declare series upperH lowerH respH
    *
    spgraph(vfields=3,hfields=2,hlabel="Figure 4. 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)
             display respH(k)
          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 5. 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)
             display respH(k)
          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)