*
* Replication file for Kilian & Vigfusson(2011), "Are the responses of
* the U.S. economy asymmetric in energy price increases and decreases?",
* Quantitative Economics, vol. 2, no 3, 419-453.
*
* Generation of figure 6 (non-linear IRF's with different sized shocks)
*
* Number of lags in VAR
*
compute p=6
*
* Number of periods for impulse responses
*
compute h=8
*
* Number of outer draws for initial conditions
*
compute nhist=300
*
* Number of inner bootstrap draws. (The paper does 10000 at this point.
* However, the difference due to the extra simulations isn't noticeable).
*
compute ndraws=100
*
* Scale-ups of the shock
*
dec vect scales(4)
compute scales=||1.0,2.0,4.0,10.0||
dec vect[series] yirf(%size(scales))
dec vect[strings] irflabels(%size(scales))
compute irflabels=||"1 s.d.","2 s.d.","4 s.d.","10 s.d."||
*
* This uses quarterly data. The monthly series (POIL and CPI) use
* last-month-of-quarter values in compacting to quarterly.
*
open data "kvqedata.rat"
calendar(q) 1973:1
data(format=rats,compact=last) 1973:01 2007:4 urate cpi pgas cgas poil gdp
*
set rpoil  = log(poil/cpi)
*
set drpoil = 100.0*(rpoil-rpoil{1})
set dgdp   = 100.0*log(gdp/gdp{1})
set dcpi = 100*(CPI-CPI{1})
set durate = 100*(urate-urate{1})

*
* This uses X and Y (as in the paper) for copies of the data series used
* in the VAR.
*
set x     = drpoil
set y     = dgdp
set z     = dcpi
set m     = durate

*
set xplus = %max(0.0,x)
*
* Estimate the two equations. The "X" equation uses only lags of X and
* Y. The Y equation uses both the X and X+ series and also includes the
* current values of each.
*
* Save the standard deviation of the residuals for the X series for use
* in sizing shocks.
*
linreg(define=xeq) x / ux
# constant x{1 to p} y{1 to p} z{1 to p} m{1 to p}
compute sigmax=sqrt(%sigmasq)
*
linreg(define=yeq) y / uy
# constant x{0 to p} y{1 to p} z{1 to p} m{1 to p} xplus{0 to p}
*

linreg(define=zeq) z / uz
# constant x{0 to p} y{1 to p}  z{1 to p} m{1 to p} xplus{0 to p}

linreg(define=meq) m / um
# constant x{0 to p} y{1 to p}  z{1 to p} m{1 to p} xplus{0 to p}



* Range of residuals as source for bootstrapping
*
compute bstart=%regstart(),bend=%regend()
*
* Range of simulated forecasts
*
compute hstart=%regend()+1,hend=%regend()+h
*
* Range of data for blocks of pre-sample values
*
compute dstart=%regstart()-p,dend=%regend()
* Definitional identities for XPLUS and for the level (rather than
* difference) of the real activity variable.
*
frml xplusdef xplus = %max(0.0,x)
frml(identity) lydef ly = ly{1}+y

*
* The values of this don't matter since everything is done as deviations
* from the base value.
*
set ly dstart hend = 0.0
set lz dstart hend = 0.0

*
* Because the Y equation includes the current value of XPLUS, it's
* necessary to arrange the model so XPLUS is solved before Y.
*
group asymmvar xeq xplusdef yeq zeq meq lydef
*
* Make copies of data for bootstrapping
*
set xraw  dstart dend = x
set yraw  dstart dend = y
set zraw  dstart dend = z
set mraw  dstart dend = m
*
do i=1,%size(scales)
   compute delta=scales(i)*sigmax
   set yirf(i) hstart hend = 0.0

   set zeros   hstart hend = 0.0
   *
   do hist=1,nhist
      *
      * Draw block of length p for initial conditions
      *
      boot(block=p) ientries hstart-p hstart-1 dstart dend
      set x hstart-p hstart-1 = xraw(ientries(t))
      set y hstart-p hstart-1 = yraw(ientries(t))
      set z hstart-p hstart-1 = zraw(ientries(t))
      set m hstart-p hstart-1 = mraw(ientries(t))

      set xplus hstart-p hstart-1 = %max(x,0)
      *
      do draw=1,ndraws
         *
         * Draw random entries for residuals for forecast period
         *
         boot entries hstart hend bstart bend
         set xpath hstart hend = ux(entries(t))
         set ypath hstart hend = uy(entries(t))
         set zpath hstart hend = uz(entries(t))
         set mpath hstart hend = um(entries(t))
         *
         * Do FORECAST with simulated paths throughout
         *
         forecast(model=asymmvar,from=hstart,to=hend,results=base,paths,noprint)
         # xpath zeros ypath zpath mpath
         *
         * Do FORECAST with first period X shock replaced by delta
         *
         set xpath hstart hend = %if(t==hstart,delta,ux(entries(t)))
         forecast(model=asymmvar,from=hstart,to=hend,results=wshock,paths,noprint)
         # xpath zeros ypath zpath mpath
         *
         * Add the gap for the accumulated Y response (4th variable) to IRF(i)
         *
         set yirf(i) hstart hend = yirf(i)+(wshock(4)-base(4))
      end do draw
   end do hist
   *
   set yirf(i) hstart hend = (yirf(i)/(nhist*ndraws))/scales(i)
end do i
*
graph(series=yirf,key=upright,klabels=irflabels,number=0,footer=$
   "Figure 6: The Response of URATE to a Positive Energy Price Shock by Shock Size")