*  MODE-HUNTING PROCEDURE  *
*
*  SOURCE: Silverman B.W.: "Density Estimation for Statistics
*     and Data Analysis", Monographs on Statistics and Applied
*     Probability 26",  Chapman & Hall, London 1986.
*
*   Written by Guido Travaglini.
*   Utilizes modified MBKERNEL.SRC by Tom Doan (1994).
*
*   This procedure finds the critical windows and the number of modes
*   of a given time series by KERNEL looping from an initially small
*   window up to the largest obtained through selected increments and
*   incremental values. Optional graphics is provided. Tests of the null
*   are obtainable from the companion procedure BOOTS.SRC. To this end,
*   the critical windows and modes are stored in a temporary file.
*
*
*   procedure(options) modes series start end maxm hini ninc inc 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
*    maxm         Maximum mode selected.
*    hini         Initial window.
*    ninc         Number of increments in window to loop starting from hini.
*    inc          Incremental value of window.
*    m            Number of bins (128 is standard)
*
procedure modes series start end maxm hini ninc inc m
*
type series     series
type integer    start end maxm ninc m
type real       hini inc
local real      hfin h0 ss mini maxh minh
local integer   i j k startl endl scs minninc sum maxmode minmode
local series    ly dyn xout yout x1 x2 xx2 imm1 hh1 mm1 mm2 summ2
local vector    h01 mod1 mm hh imm critwin critmode
local string    ghead gsubhead head subhead
option switch   scatter 0
option switch   spgraph 0
dim h01(ninc) mod1(ninc) mm(ninc) hh(ninc) imm(ninc) critwin(ninc) critmode(ninc)
*
inquire(series=series) startl>>start endl>>end
DISPLAY ' '
DISPLAY @2 '***********************************'
DISPLAY @2 '*      MODE-HUNTING PROCEDURE     *'
DISPLAY @2 '***********************************'
DISPLAY ' '
com hfin = hini+inc*ninc
disp @5 'You selected an initial window size of' @43 hini
disp @5 'Your window increment is' @33 #.#### inc
disp @5 'You selected a total number of' ninc @40 'increments'
disp @5 'Your final window will be' @34 #.#### hfin
*
*clear all
do i = 1,ninc
  com h0 = hini+i*inc
  @kernel(kernel=gaussian,window=h0) series start end xout yout
  set dyn = yout-yout{1}
  *Stage 1: binarize dyn by letting values gt. 0 be = 1.0 and others = 0.0
  set x1 2 m = %if(dyn.gt.0.0,1.0,0.0)
  *Stage 2: make 1st. diff of such values and square them to check for 0/1 changes
  set x2 2 m = (x1-x1{1})**2
  *Stage 3: compute total # of sign changes (0/1)
  acc x2 2 m xx2
  com scs = fix(xx2(m))              ;* Computing sign changes (scs)
  do j = 1,maxm+1
   com ss = 2*j-1
   if scs==ss
   com mod = j                       ;* Associating modes to sign changes
  end do j
  com h01(i) = h0, mod1(i) = mod
  com mm(i) = (mod1(i)-mod1(i-1))**2   ;* Checking for modes
  com hh(i) = h01(i)*mm(i)
  {
  if mm(i).eq.1.0
  com imm(i) = i*mm(i)
  }
end do i
*
smpl 1 ninc
set imm1 = imm(t)
ext(noprint) imm1
com mini = %minvalue(imm1), minninc = fix(mini) ;*Establish first valid ninc
smpl minninc ninc
set hh1 = h01(t)
set mm1 = mod1(t)
set mm2 = mm(t)
acc mm2 minninc ninc summ2
com sum = fix(summ2(ninc))
com maxh = %maxvalue(hh1), minh = %minvalue(hh1)
disp(unit=output) @5 'Windows found range from' @30 #.#### minh 'to' @40 #.#### maxh
com maxmode = fix(%maxvalue(mm1)), minmode = fix(%minvalue(mm1))
disp(unit=output) @5 'Modes found range from'  maxmode 'to'  minmode
{
if sum.gt.maxmode
disp(unit=output) @5 'Mode jump detected: you must tighten inc value or start from higher hini'
if maxmode.gt.maxm
disp @5 'Highest mode found exceeds maximum mode initially set'
}
*
open temp temp.dat      ;*Opens temp file where critical windows are stored.
do j = 1,1
   disp
   disp @1 'Critical window' ##.#### @20  'Number of modes'
   disp '***************************************'
   disp
 do i = minninc,ninc
   com critwin(i) =  %if(hh(i).gt.0.0,hh(i),0.0)
   com critmode(i) = mm1(i)
   if critwin(i).gt.0.0.and.mm(i).eq.1.0
   disp @4 ##.#### critwin(i) @26 fix(critmode(i))
   if critwin(i).gt.0.0.and.mm(i).eq.1.0
   disp(unit=temp) ##.#### critwin(i) fix(critmode(i))
 end do i
end do j
close temp
*
if scatter
scatter(header='Critical window width and number of modes',style=lines, $
hlabel='Critical window width',vlabel='Number of modes') 1
#hh1 mm1 minninc ninc
*
***Graphing alternative critical window kernels***
*

disp(store=ghead) 'Gaussian KDEs and critical windows of' %label(series)
disp(store=gsubhead) %label(series)
if spgraph
{
spgraph(hfield=3,vfield=2,head=ghead,vlabel='Frequency',hlabel=%label(series))
do i = minninc,ninc
   com critwin(i) =  %if(hh(i).gt.0.0,hh(i),0.0)
   com critmode(i) = mm1(i)
   disp(store=head) 'Critical window:' ##.#### critwin(i)
   disp(store=subhead) 'Modes:' # critmode(i)
   if critwin(i).gt.0.0.and.mm1(i).le.maxm
   {
   @kernel(kernel=gaussian,window=critwin(i)) series / xout yout
   graph(head=head,subhead=subhead,style=poly) 1
   #yout 1 m
   }
end do i
spgraph(done)
}
*
smpl
end modes



