*
* This does a parametric bootstrap (to get error bands for an IRF) of a
* VECM with known cointegrating vector. The model used is from King,
* Plosser, Stock and Watson(1991), "Stochastic Trends and Economic
* Fluctuations", AER, vol 81, pp 819-840.
*
compute ndraws=500
compute nvar  =3
compute nstep =25
*
open data kpswdata.rat
calendar(q) 1947
data(format=rats) 1947:1 1988:4 c in y mp dp r
*
* Create the equations describing the cointegrating relations (c-y and
* in-y). This means that there is only one stochastic trend linking the
* three variables.
*
equation(coeffs=||1.0,-1.0||) covery
# c y
equation(coeffs=||1.0,-1.0||) iovery
# in y
*
* Create the VECM on the original variables
*
system(model=vecmmodel)
variables y c in
lags 1 to 9
det constant
ect covery iovery
end(system)
*
* Estimate the model, saving the residuals and the estimation range.
*
estimate(resids=baseresids)
compute bstart=%regstart(),bend=%regend()
*
* Substitute out for the error correction terms
*
compute basevecm=%modelsubstect(vecmmodel)
*
source varbootsetup.src
source forcedfactor.src
*
* Create a parallel VECM for the bootstrapped series. This uses the
* variable names from @VARBootSetup so @VARBootDraw can be used to do
* the actual draws.
*
dec vect[series] %%VARResample(nvar)
do i=1,nvar
   set %%VARResample(i) = %modeldepvars(basevecm)(i){0}
end do i
equation(coeffs=||1.0,-1.0||) coveryboot
# %%VARResample(2) %%VARResample(1)
equation(coeffs=||1.0,-1.0||) ioveryboot
# %%VARResample(3) %%VARResample(1)
*
system(model=bootvecm)
variables %%VARResample
lags 1 to 9
det constant
ect coveryboot ioveryboot
end(system)
*
* To isolate the balanced growth shock
*
compute atilde=||1.0|1.0|1.0||
*
* Bookkeeping series for the IRF for the balanced growth shock and the
* error decomposition percentages.
*
dec rect[series] %%responses(nstep,nvar)
dec rect[series] errors(nstep,nvar)
do i=1,nstep
   do j=1,nvar
      set %%responses(i,j) 1 ndraws = 0.0
      set errors(i,j) 1 ndraws    = 0.0
   end do j
end do i
*
infobox(action=define,progress,lower=1,upper=ndraws) "Bootstrapping"
do draws = 1,ndraws
   *
   * Bootstrap the data. This enforces the parametrization that the
   * three series are cointegrated.
   *
   @VARBootDraw(model=basevecm,resids=baseresids) bstart bend
   *
   * Estimate the VECM on the bootstrapped data
   *
   estimate(noprint) bstart bend
   *
   * Compute the responses to unit shocks and compute an estimate of the
   * long-run response as the 500th element of the IRF.
   *
   impulse(model=bootvecm,factor=%identity(3),results=baseimp,noprint,steps=500)
   compute lrsum=%xt(baseimp,500)
   *
   * Compute a factor with the balanced growth shock as the first shock.
   * (We don't care about the rest).
   *
   compute d=inv(%innerxx(atilde))*tr(atilde)*lrsum
   @forcedfactor(force=row) %sigma d f
   compute lrfactor=lrsum*f
   impulse(noprint,model=bootvecm,factor=f/lrfactor(1,1),results=impulses,steps=25)
   *
   * Store the impulse responses. In this case, we're only interested in
   * the responses to the first shock.
   *
   do i=1,nstep
      do j=1,nvar
         set %%responses(i,j) draws draws = impulses(j,1)(i)
      end do j
   end do i
   *
   * Store the decomposition of variance. Again, we're only interested
   * in the fraction explained by the first shock.
   *
   errors(noprint,model=bootvecm,factor=f,results=decvar,steps=25)
   do i=1,nstep
      do j=1,nvar
         set errors(i,j) draws draws = decvar(j,1)(i)
      end do j
   end do i
   infobox(current=draws)
end do draws
infobox(action=remove)
*
report(action=define)
report(atrow=1,atcol=1,tocol=4,span) "Fraction of the forecast-error variance"
report(atrow=2,atcol=1,tocol=4,span) "attributed to the real permanent shock"
report(atrow=3,atcol=1,align=center) "Horizon" "y" "c" "i"
compute row=1
dofor horizon = 1 4 8 12 16 20 24
   compute row=row+3
   report(atrow=row,atcol=1) horizon
   do j=1,nvar
      stats(noprint) errors(horizon,j) 1 ndraws
      report(atrow=row,atcol=j+1) %mean
      report(atrow=row+1,atcol=j+1,special=parens) sqrt(%variance)
   end do j
end do horizon
report(action=format,picture="*.##",atrow=4,align=decimal)
report(action=show)
*
*
dec vect[series] mid(nvar) upper(nvar) lower(nvar)
do j=1,nvar
   do i=1,nstep
      sstats(mean) 1 ndraws %%responses(i,j)>>first %%responses(i,j)^2>>second
      compute stddev=sqrt(second-first^2)
      set mid(j)   i i = first
      set upper(j) i i = first+stddev
      set lower(j) i i = first-stddev
   end do i
end do j
*
table(noprint) 1 nstep upper lower
*
spgraph(vfields=nvar,ylabels=||"Output","Consumption","Investment"||,$
  footer="Figure 2 - Responses to Shock in Real Permanent Component")
do j=1,nvar
   graph(max=%maximum,min=%minimum,nodates) 3
   # mid(j)
   # upper(j) / 2
   # lower(j) / 2
end do j
spgraph(done)