The RATS Software Forum

Code: Select all

procedure acKupiecBlockBoot trigger start end
*
   type series    trigger
   type integer   start end
*
   local integer           startl endl draw
   local integer           Nobs kobs kboot
   local integer           blocklen
   local real              violratio alphahat alphahatboot
   local real              pofobs pofboot countextreme
   local series[integer]   bootentries
   local series            H0DIST bootviol
   local real              acsum acsum_SS res_S_M res_S_MSS res_S_SE
*
   option real    alpha       .01
   option real    siglevel    .05
   option integer ndraws      10000
   option integer horizon     1
*
   local series   temp gPOFBoot fPOFBoot
   local series   pofbootdist
*
*
   inquire(reglist) startl<<start endl<<end
   # trigger
*
* --- Set block length ---
   compute blocklen = horizon
*
* --- Calculate observed statistics ---
   sstats(mean) startl endl trigger>>alphahat
   compute Nobs = %nobs
   compute kobs = fix(alphahat*Nobs)
   compute violratio = alphahat
*
* --- Observed POF statistic ---
   if (kobs>0.and.kobs<Nobs)
      compute pofobs=2*Nobs*((1-alphahat)*log((1-alphahat)/(1-alpha))+alphahat*log(alphahat/alpha))
   else if (kobs==0)
      compute pofobs=-2*Nobs*log(1-alpha)
   else
      compute pofobs=-2*Nobs*log(alpha)
*
*
* --- BLOCK-BOOTSTRAPPING ---
*
*  --- Calculate p-value ---
   compute countextreme = 0.0

   compute acsum = 0.0
   compute acsum_SS = 0.0

*
* --- Block Bootstrap series ---
   set pofbootdist 1 ndraws = 0.0

* --- Block Bootstrap loop ---
   do draw=1,ndraws
*
*     STEP 1: Generate fresh H0DIST for this replication
      set H0DIST startl endl = %if(%ranflip(alpha),1.0,0.0)
*
*     STEP 2: Resample H0DIST in blocks to create serial dependence
      BOOT(block=blocklen) bootentries startl endl startl endl
      set bootviol startl endl = H0DIST(bootentries(t))
*
*     STEP 3: Calculate bootstrap POF statistic
      sstats(mean) startl endl bootviol>>alphahatboot
      compute kboot = fix(alphahatboot*Nobs)
*
      if (kboot>0.and.kboot<Nobs)
         compute pofboot=2*Nobs*((1-alphahatboot)*log((1-alphahatboot)/(1-alpha))+alphahatboot*log(alphahatboot/alpha))
      else if (kboot==0)
         compute pofboot=-2*Nobs*log(1-alpha)
      else
         compute pofboot=-2*Nobs*log(alpha)
*
*     Store pof bootstrap to a vector pofbootdist for POF statistics and density plot
      compute pofbootdist(draw) = pofboot
*
*     STEP 4: Count if bootstrap statistic > observed
      if (pofboot>pofobs)
         compute countextreme = countextreme + 1.0
*
*     Alternatively
      compute acsum = acsum+pofboot
      compute acsum_SS = acsum_SS+(pofboot)^2
*
   end do draw
*
*
* --- Store results ---
   compute %%pofboot = countextreme/ndraws
*
*
* https://estima.com/RATSBetaHelp/topics/grangerbootstraprpf.html
* Compute the percentiles of the test statistics
   stats(fractiles) pofbootdist 1 ndraws
   compute %%pofbootmean = %mean
   compute %%pofbootse = sqrt(%variance)
*
* CONFIRMATION
* Figure out the bootstrapped p-value (SAME AS %%pofboot)
   sstats(mean) 1 ndraws (pofbootdist>pofobs)>>pvalue
   disp "Bootstrapped p-value" pvalue
*

* Alternatively
   compute res_S_M = acsum/ndraws
   compute res_S_MSS = acsum_SS/ndraws
   compute res_S_SE = sqrt(res_S_MSS-(res_S_M)^2)
   disp 'res_S_M' res_S_M
   disp 'res_S_SE' res_S_SE


* --- Plot density function for pofbootdist ---
*  convert VECTOR to a SERIES for plotting
   set temp 1 ndraws = pofbootdist(t)

   density(smoothing=2.0) temp 1 ndraws gpofboot fpofboot

   * 10 horizons, 10 colours
   compute col = %if(horizon==1,2,$
                     %if(horizon==2,4,$
                     %if(horizon==3,5,$
                     %if(horizon==4,6,$
                     %if(horizon==5,7,$
                     %if(horizon==6,8,$
                     %if(horizon==7,9,$
                     %if(horizon==8,10,$
                     %if(horizon==9,11,12)))))))))
   spgraph(hea="",vfi=1,hfi=1)
   scatter(key=upright,head="A non-parametric empirical approximation to the sampling distribution\\of the POF statistic under H0 with serial dependence",$
            footer="horizon="+%strval(horizon,"#")+",   "+$
                   "POF mean under H0 (bootstrap):"+%strval(%%pofbootmean,".####")+",   "+$
                   "p-value (bootstrap):"+%strval(%%pofboot,".####"),$
            hgrid=||%%pofbootmean||,style=line) 1
   # gpofboot fpofboot 1 ndraws col
   spgraph(done)
*
*
* --- Output Results ---
   display "======================================="
   display "--- Block Bootstrap Kupiec POF Test ---"
   display "======================================="
   display "From "+%datelabel(startl)+" to "+%datelabel(endl)
   display "Horizon (h):" horizon
   display "Block length:" blocklen
   display

   display "--- Sample Statistics ---"
   display "Sample size (N):" Nobs
   display "Observed exceedances (k):" kobs
   display "Violation ratio (k/N):" violratio
   display "Target alpha:" alpha
   display

   display "--- Test Configuration ---"
   display "Significance level:" siglevel
   display "Bootstrap replications:" ndraws
   display

   display "--- Test Results ---"
   display "POF statistic (observed):" pofobs
   display "POF mean under H0 (bootstrap) approximately 1.0:" %%pofbootmean
   display "POF se under H0 (bootstrap):" %%pofbootse
   display "p-value (bootstrap):" %%pofboot
   display "=========================================="
*
end acKupiecBlockBoot