Identifying VARs with sign restrictions

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

[TomDoan]

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

If you want the 3rd shock to hit the 3rd variable with size one, add:

Code: Select all

compute i3=i3/i3(3)

If you don't want to rescale the 2nd, leave i2 alone.

[TL]

Dear Tom,

Thank you very much for your reply. I am sorry to ask you again.

Under the following line in the original code

Code: Select all

compute i1=p*v1

I have added the following code (for the 1st shock to hit the 3rd variable with size one)

Code: Select all

compute i1=i1/i1(3)

I got the following message

## SX20. Expected , Here
>>>>compute i1=i1/i1(3)<<<<

I am not sure whether I have added such code at the right place or I might do something wrong. Could you please give me suggestions?

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
           compute i1=i1/i1(3)
           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)

[TomDoan]

Since i1 is never declared, it's created as a n x 1 rectangular. If you add

Code: Select all

declare vector i1 i2 i3

before starting the loop, it should work.

[TL]

Dear Tom,

I am sorry to ask you again. I have added the following code

Code: Select all

declare vector i1 i2

and

Code: Select all

compute i1=i1/i1(3)

However, the impulse response acquired does not show that the 3rd variable is hit by the 1st shock with size one. In particular, the impulse response acquired is not different from one acquired from the original code not having those two lines.

I attach you the code added those two lines.

Could you please give me suggestions?

Thank you very much for your help and kindness.

Sincerely,

Tim

[TomDoan]

You need to pay careful attention to the difference between the v vectors, which are in rotated space, and the i vectors, which are in standard variable space. Rescaling the i vector gives you the correct loadings, but you're then using (later on) the v vector to get the impulse responses. If you change i, but not v, you get, not surprisingly, the same results as if you didn't change i. You need your weights on the impulse responses to be
inv(p)*i1, not v1. (You don't want to change v1 itself since you need the unscaled version to get the FEVD).

[TL]

Dear Tom,

Thank you very much for your help. I really appreciate your help.

I have changed the following code

Code: Select all

compute ik=%xt(impulses,k)*v1

to

Code: Select all

compute ik=%xt(impulses,k)*inv(p)*i1

and changed the following code

Code: Select all

ewise goodresp(accept)(i,j)=(ik=%xt(impulses,i)*v1),ik(j)

to

Code: Select all

ewise goodresp(accept)(i,j)=(ik=%xt(impulses,i)*inv(p)*i1),ik(j)

The impulse response acquired does show that the 3rd variable is hit by the 1st shock with size one. However, it seems that the dashed lines representing the 16 percent and 84 percent quantiles of the posterior distribution of impulse responses are unscaled. I am not sure if it has something to do with the spgraph section.

I attach you the code changed as above.

Could you please give me suggestions?

Thank you very much for your help and kindness.

Sincerely,

Tim

[TomDoan]

That looks correct. Why would you think the bounds are for unscaled responses?

[TL]

Dear Tom,

I am asking this because I see literature normalizing the impulse response to increase the variable by one and also showing probability bands corresponding to such normalized responses.

Nevertheless, it is fine if we would not have scaled bounds following scaled impulse responses.

I would like to thank you very much for all your help and suggestions. I really appreciate your kindness.

Thank you very much.

Sincerely,

Tim

[TomDoan]

I'm confused. I just ran your program and it gives perfectly reasonable-looking bands. Why do you think it's wrong?

[TL]

Dear Tom,

It is my mistake.

I have run my program again and, as you said, it gives perfectly reasonable-looking bands.

I am sorry for this.

Thank you very much for all your help and suggestions. I really appreciate your kindness.

Sincerely,

Tim