OPEN DATA shockdata.xls
    CALENDAR(M) 1986 1
    compute missc=1.0e+32
    data(format=xlS,org=columns) 1986:1 2010:12 Firon IT Piron Qiron ST YW
    
    SYSTEM(MODEL=VARMODEL)
    VARIABLES Qiron Piron YW IT Firon ST
    LAGS 1 TO 6
    END(SYSTEM)
    ESTIMATE(noprint)
    *
    dec vect[strings] vl(6)
    compute vl=||'iron production','price of iron','world economic activity',$
                 'iron inventories','future price of iron','future-spot spread'||
*
* 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=1000
    compute n2=1000
    compute nkeep=1000  
    compute nvar=6
    compute nstep=20
    compute KMAX=5

    *************************************************************************************
    compute hstart = 1986:1
    compute hend   = 2010: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)
    dec vect[series] upathsaa(nvar) baseplusaa(nvar)
    dec vect[series] upathsaaa(nvar) baseplusaaa(nvar)
    clear upaths
    clear upathsa
    clear upathsaa
    clear upathsaaa
    *************************************************************************************
    *
    * This is the standard setup for MC integration of an OLS VAR
    *
    dec symm s(nvar,nvar)
    dec vect v1(nvar)      ;* For the unit vector on the 1st draw
    dec vect v2(5)         ;* For the unit vector on the 2nd draw
    dec vect v3(4)         ;* For the unit vector on the 3nd draw
    dec vect v4(3)         ;* For the unit vector on the 4nd draw
    dec vect v(nvar)       ;* Working impulse vector

    compute sxx    =%decomp(%xx)
    compute svt    =%decomp(inv(%nobs*%sigma))
    compute betaols=%modelgetcoeffs(varmodel)
    compute ncoef  =%rows(sxx)
    compute wishdof=%nobs-ncoef
    dec rect ranc(ncoef,nvar)
    *
    * Most draws are going to get rejected. We allow for up to 1000
    * good ones. The variable accept will count the number of accepted
    * draws. GOODRESP will be a RECT(nsteps,nvar) at each accepted
    * draw.
    *
    declare vect[rect] goodresp(nkeep) goodfevd(nkeep)
    declare vect[rect] goodrespa(nkeep) goodfevda(nkeep)
    declare vect[rect] goodrespaa(nkeep) goodfevdaa(nkeep)
    declare vect[rect] goodrespaaa(nkeep) goodfevdaaa(nkeep)
    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) goodhistaa(1000)  goodhistaaa(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 - supply shock
        ************************************************
           compute v1=%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).
      *
           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(5)>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(nvar-1)
           compute i2=%xsubmat(f,1,nvar,2,nvar)*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(5)>0.or.ik(6)<0
                 branch 105
           end do k


           @forcedfactor(force=column) sigmad i1~i2 f
           compute v3=%ransphere(nvar-2)
           compute i3=%xsubmat(f,1,nvar,3,nvar)*v3
           compute v3=inv(p)*i3

           *****************************************************
           *  third Set of Restrictions - specific Shock
           *****************************************************
           do k=1,KMAX+1
              compute ik=%xt(impulses,k)*v3
              if ik(1)>0.or.ik(2)>0.or.ik(3)<0.or.ik(5)>0.or.ik(6)<0
                 branch 105
           end do k
           @forcedfactor(force=column) sigmad i1~i2~i3 f
           compute v4=%ransphere(nvar-3)
           compute i4=%xsubmat(f,1,nvar,4,nvar)*v4
           compute v4=inv(p)*i4
           
           *****************************************************
           *  Forth Set of Restrictions - speculation Shock
           *****************************************************
           do k=1,KMAX+1
              compute ik=%xt(impulses,k)*v4
              if ik(5)>0.or.ik(6)>0
                 branch 105
           end do k
           
           
           *
           * Meets all the restrictions
           *
           compute accept=accept+1
           dim goodresp(accept)(nstep,nvar) goodfevd(accept)(nstep,nvar)
           dim goodrespa(accept)(nstep,nvar) goodfevda(accept)(nstep,nvar)
           dim goodrespaa(accept)(nstep,nvar) goodfevdaa(accept)(nstep,nvar)
           dim goodrespaaa(accept)(nstep,nvar) goodfevdaaa(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)
           ewise goodrespaa(accept)(i,j)=(ik=%xt(impulses,i)*v3),ik(j)
           ewise goodfevdaa(accept)(i,j)=(ik=(irfsquared(i)*(v3.^2))./(irfsquared(i)*ones)),ik(j)
           ewise goodrespaaa(accept)(i,j)=(ik=%xt(impulses,i)*v4),ik(j)
           ewise goodfevdaaa(accept)(i,j)=(ik=(irfsquared(i)*(v4.^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) goodhistaa(accept)(nhor,nvar)  goodhistaaa(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
          compute remap=p*v3*tr(v3)*inv(p)
          do t=hstart,hend
             compute %pt(upathsaa,t,remap*%xt(u,t))
          end do t
          compute remap=p*v4*tr(v4)*inv(p)
          do t=hstart,hend
             compute %pt(upathsaaa,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
          forecast(model=varmodel,from=hstart,to=hend,results=baseplusaa,paths)
          # upathsaa
          forecast(model=varmodel,from=hstart,to=hend,results=baseplusaaa,paths)
          # upathsaaa
          *
          * 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)
          ewise goodhistaa(accept)(i,j)=baseplusaa(j)(i+hstart-1)-base(j)(i+hstart-1)
          ewise goodhistaaa(accept)(i,j)=baseplusaaa(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=6,hfields=4,hlabel='Impluse Response with Pure-Sign Approach',subhea="Supply shock")
    do i=1,nvar
        compute minlower=maxupper=0.0
        smpl 1 accept
        do k=1,nstep
           set work = goodresp(t)(k,i)
           compute frac=%fractiles(work,||.16,.84||)
           compute lower(k)=frac(1)
           compute upper(k)=frac(2)
           compute resp(k)=%avg(work)
        end do k
    *
        smpl 1 nstep
        graph(pattern,ticks,number=0,picture='##.##',header='Impluse Response'+vl(i)) 3
        # resp
        # upper / 2
        # lower / 2
    end do i
    *
    spgraph(vfields=6,hfields=4,hlabel='Impluse Response with Pure-Sign Approach',subhea="demand shock")
    do i=1,nvar
        compute minlower=maxupper=0.0
        smpl 1 accept
        do k=1,nstep
           set work = goodrespa(t)(k,i)
           compute frac=%fractiles(work,||.16,.84||)
           compute lower(k)=frac(1)
           compute upper(k)=frac(2)
           compute resp(k)=%avg(work)
        end do k
    *
        smpl 1 nstep
        graph(pattern,ticks,number=0,picture='##.##',header='Impluse Response'+vl(i)) 3
        # resp
        # upper / 2
        # lower / 2
    end do i
    *
    spgraph(vfields=6,hfields=4,hlabel='Impluse Response with Pure-Sign Approach',subhea="specific shock")
    do i=1,nvar
        compute minlower=maxupper=0.0
        smpl 1 accept
        do k=1,nstep
           set work = goodrespaa(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='Impluse Response'+vl(i)) 3
        # resp
        # upper / 2
        # lower / 2
    end do i
    *
    spgraph(vfields=6,hfields=4,hlabel='Impluse Response with Pure-Sign Approach',subhea="speculation shock")
    do i=1,nvar
        compute minlower=maxupper=0.0
        smpl 1 accept
        do k=1,nstep
           set work = goodrespaaa(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='Impluse Response'+vl(i)) 3
        # resp
        # upper / 2
        # lower / 2
    end do i
    spgraph(done)
    smpl
    *
    clear upper lower resp

*
    
    *
    clear upper lower resp
    *
*************************************************************************************
    *
    spgraph(vfields=3,hfields=2,hlabel='Fraction of Variance Explained with Pure-Sign Approach',subhea='supply shock')
    do i=1,nvar
        compute minlower=maxupper=0.0
        smpl 1 accept
        do k=1,nstep
           set work = goodfevd(t)(k,i)
           compute frac=%fractiles(work,||.16,.50,.84||)
           compute lower(k)=frac(1)
           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='Fraction of Variance Explained with Pure-Sign Approach',subhea='demand shock')
    do i=1,nvar
        compute minlower=maxupper=0.0
        smpl 1 accept
        do k=1,nstep
           set work = goodfevda(t)(k,i)
           compute frac=%fractiles(work,||.16,.50,.84||)
           compute lower(k)=frac(1)
           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='Fraction of Variance Explained with Pure-Sign Approach',subhea='specific shock')
    do i=1,nvar
        compute minlower=maxupper=0.0
        smpl 1 accept
        do k=1,nstep
           set work = goodfevdaa(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='Fraction of Variance Explained with Pure-Sign Approach',subhea='speculation shock')
    do i=1,nvar
        compute minlower=maxupper=0.0
        smpl 1 accept
        do k=1,nstep
           set work = goodfevdaaa(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 A. 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 supply 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 B. 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 demand shocks to "+vl(i)) 3
       # respH hstart hend
       # upperH / 2
       # lowerH / 2

    end do i
    *
    spgraph(done)

clear upperH lowerH respH


spgraph(vfields=3,hfields=2,hlabel="Figure C. 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 = goodhistaa(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 specific shocks to "+vl(i)) 3
       # respH hstart hend
       # upperH / 2
       # lowerH / 2

    end do i
    *
    spgraph(done)

clear upperH lowerH respH

spgraph(vfields=3,hfields=2,hlabel="Figure D. 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 = goodhistaaa(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 speculation shocks to "+vl(i)) 3
       # respH hstart hend
       # upperH / 2
       # lowerH / 2

    end do i
    *
    spgraph(done)