*
* GRANGERBOOTSTRAP.RPF
* RATS Version 8.
* Bootstrap alternative to CAUSAL.RPF.
*
cal(m) 1959:1
allocate 2006:4
open data haversample.rat
data(format=rats) / gdph fm1
log gdph
log fm1
*
* Granger causality test with bootstrapped p-value.
*
* The bootstrap draws are conditional on the first 8 values of the two
* series. We do the equations in terms of "sample" series to make the
* bookkeeping easier later on.
*
set msample = fm1
set gsample = gdph
linreg(define=eqm) msample / rm
# constant msample{1 to 8}
linreg(define=eqg) gsample / rg
# constant gsample{1 to 8} msample{1 to 8}
*
group noncausal eqm eqg
*
compute ndraws=10000
set stats 1 ndraws = 0.0
do draw=1,ndraws
   *
   * Resample the residuals over the regression range
   *
   boot entries %regstart() %regend()
   set rmsample = rm(entries(t))
   set rgsample = rg(entries(t))
   forecast(paths,model=noncausal,from=%regstart(),to=%regend(),results=results)
   # rmsample rgsample
   set msample = results(1)
   set gsample = results(2)
   *
   * Do the causality test on the resampled data. Except for minor scale
   * factors, this will give identical results to a full system LR test
   * (under the assumptions of homoscedasticity).
   *
   linreg(noprint) msample
   # constant msample{1 to 8} gsample{1 to 8}
   exclude(noprint)
   # gsample{1 to 8}
   compute stats(draw)=%cdstat
end do draw
*
* Compute the percentiles of the test statistics
*
stats(fractiles) stats 1 ndraws
*
* Do the test with the real data
*
linreg(noprint) fm1
# constant fm1{1 to 8} gdph{1 to 8}
exclude(print)
# gdph{1 to 8}
*
* Figure out the bootstrapped p-value
*
sstats(mean) 1 ndraws (stats>%cdstat)>>pvalue
*
disp "Bootstrapped p-value" pvalue