* MODE BOOTSTRAPPING PROCEDURE *
*
*  SOURCES: Silverman B.W.: "Density Estimation for Statistics
*     and Data Analysis", Monographs on Statistics and Applied
*     Probability 26",  Chapman & Hall, London 1986;
*     Efron B. & Tibshirani R.J., "An Introduction to the Bootstrap",
*     Monographs on Statistics and Applied Probability 57", Chapman
*     & Hall, London 1993.
*
*   Written by Guido Travaglini.
*
*   This procedure is companion to MODES.SRC. It retrieves its
*   critical windows stored in a temp file and KERNEL boots the
*   given time series for each of them. It finally tests for
*   the null of m versus the alternative of m+1 modes by supplying
*   the statistic ASL (Achieved Significance Level).
*
*
*   procedure boots series start end h1 h0modes nboots m
*
*   Parameters:
*
*   series       Series to analyze.
*   start end    Time range. If working with balanced panel series of N
*                individuals  and T observations, set start = 1, end = NT
*   h1           Window.
*   h0modes      Number m of modes to test fo the null.
*   nboots       Number of booted samples to KERNEL with h1.
*   m            Number of bins (128 is standard)
*
procedure boots series start end h1 h0modes nboots m
*
type series     series
type integer    start end nboots m h0modes
type real       h1
local real      ASL scs
local integer   i j mod ss startl endl
* Change for compatibility with RATS 6:
* local index     tek
local series[integer] tek
local series    ext yb z xout yout dynb x1 x2 xx2 hmod ax axs
local rec       bbb
*
dim bbb(end,1)
inquire(series=series) startl>>start endl>>end
DISPLAY ' '
DISPLAY @2 '***********************************'
DISPLAY @2 '*        BOOTSTRAPPING            *'
DISPLAY @2 '***********************************'
DISPLAY ' '
*
stat(noprint) series startl endl
com bbb = %ran(1.0)
set ext = bbb(t,1)
set hmod 1 nboots = %na
do i = 1,nboots
  boot tek /
  set yb = series(tek(t))
  stat(noprint) yb
  set z = %mean+( (yb-%mean+h1*ext)/(sqrt(1+(h1**2/%variance))) )
  @kernel(window=h1,kernel=gaussian) z / xout yout
  set dynb = yout-yout{1}
  set x1 2 m-1 = %if(dynb.gt.0.0,1.0,0.0)
  set x2 3 m-1 = (x1-x1{1})**2
  acc x2 2 m-1 xx2
  com scs = xx2(m-1)
   do j = 1,maxm
    com ss = 2*j-1
    if scs==ss
     com mod = j
   end do j
 compute hmod(i) = mod
end do i
*
* Create binary series 'axs' to test for H0; i.e. write 1.0 (0.0) if booted modes exceed (equal) H0 modes.
disp @5 'Null modes:' @17 # h0modes
disp @5 'Critical window:' @22 #.#### h1
set ax 1 nboots = %if(hmod.gt.h0modes,1.0,0.0)
acc ax 1 nboots axs
com ASL = (axs(nboots)/nboots)          ;*Percentage of booted modes exceeding H0 mode.
disp @5 'ASL' @12 #.### ASL @20
end boots

