*
* @ICSS(options)  x  start end
*
* Searches for breaks in variance using
*
* Inclan and Tiao, "Use of Cumulative Sums of Squares for
* Retrospective Detection of Changes in Variance", JASA 1994, vol
* 89, pp 913-923
*
* with refinement to allow for non-Gaussian distribution from
*
* Sanso, Arago and Carrion-i-Silvestre(2004), "Testing for changes in the
* unconditional variance of financial time series", Spanish Review
* of Financial Economics, vol 4, pp 32-53.
*
* Parameters:
*   x = series to analyze
*   start end = range to analyze. By default, the defined range of x.
*
* Options:
*   SIGNIF=.01/[.05]/.10
*   CV=critical value for determining variance change [not used]
*     This determines the critical value used for determining whether a
*     subsample split should be accepted. If CV isn't used, the SIGNIF
*     option controls (with a default of .05), which uses the asymptotic
*     critical values (1.358 for .05, 1.628 for .01 and 1.224 for .10). CV
*     can be used for inputting a specific alternative value.
*
*   BREAKS=(output) VECT[INT] of entries with breaks located [not used]
*   SUBSAMPLES=(output) RECT[INT] with subsamples [not used]
*     For the array created by SUBSAMPLES, ss(1,i) is the start entry of the
*     ith subsample and ss(2,i) is the end entry. Note that there if there
*     are K breaks, there will be K+1 subsamples. A breakpoint marks the end
*     of a subsample.
*
*   [PRINT]/NOPRINT
*   TITLE="title for report" ["ICSS Analysis of Series x"]
*
*   ROBUST/[NOROBUST]
*     If ROBUST, corrects the statistics for a non-Gaussian distribution
*     (using sample estimates of the 4th moment).
*
*   TRACE/[NOTRACE]
*     If TRACE, information showing the decisions taken in coming up with
*     sample breaks.
*
* Revision Schedule:
*   03/2005 Written by Tom Doan, Estima
*   04/2008 Stripped out of program file.
*   09/2010 Use REPORT
*   05/2011 Add BREAKS, PRINT and TRACE options. Add "Step 3" for pruning
*           and adjusting large set of breaks identified in the direct analysis.
*   10/2017 Add SIGNIF, CV and SUBSAMPLES options
*   04/2018 Add ROBUST option
*
*******************************************************************************************
*
* Procedure IncTiaoBreak is called from ICSS1and2 to do most of the
* work. It looks for a single variance break in x in the interval
* [start,end]. If a break which exceeds the critical value is found,
* %maxent is set equal to the entry at which that occurs. If no such
* break is found, %maxent will be 0 on return.
*
procedure IncTiaoBreak x start end
type     series   x
type     integer  start end
*
option   real     cv       1.358
option   switch   robust   0
*
local series   ck
local real     nobs ct sigsq eta4 kappa1
*
* This uses the notation from the Sanso, et al, paper
*
sstats(mean) start end x^2>>sigsq x^4>>eta4
set ck start end = x^2
acc ck start end
compute ct=ck(end)
*
* If not robust, use 3 sig^4 for Gaussian distribution
*
if .not.robust
   compute eta4=3*sigsq^2
compute nobs=(end-start+1)
sstats(max) start end abs((ck-ct*(t-start+1)/nobs))>>kappa1
compute kappa1=sqrt(1.0/(nobs*(eta4-sigsq^2)))*kappa1
if kappa1>cv
   ;
else
   compute %maxent=0
end IncTiaoBreak
*
*******************************************************************************************
*
* ICSS1and2 uses repeated calls to IncTiaoBreak to do steps 1 and 2 to
* locate possible breaks between start and end. It returns a VECT[INT]
* which will have dimension 0 (if no significant break occurs), 1 (if
* only a single break occurs), or 2 (if there is more than one
* significant break, in which case, it returns the most extreme break
* values).
*
function ICSS1and2 x start end cv robust
type vect[int] ICSS1and2
type series    x
type integer   start end
type real      cv
type integer   robust
*
local integer t1 t2 kstar kfirst klast
*
compute t1=start,t2=end
@IncTiaoBreak(cv=cv,robust=robust) x start end
if %maxent==0 {
   dim ICSS1and2(0)
   return
}
compute t2=kstar=%maxent
loop
   @IncTiaoBreak(cv=cv,robust=robust) x t1 t2
   if %maxent>0
      compute t2=%maxent
   else
      break
end loop
compute kfirst=t2
compute t1=kstar+1,t2=end
loop
   @IncTiaoBreak(cv=cv,robust=robust) x t1 t2
   if %maxent>0
      compute t1=%maxent+1
   else
      break
end loop
compute klast=t1-1
*
if kfirst==klast
   compute ICSS1and2=||kfirst||
else
   compute ICSS1and2=||kfirst,klast||
end
*
*********************************************************************************
*
* This is the controlling routine for the ICSS algorithm. Because the
* ICSS algorithm is recursive, and RATS doesn't allow recursive
* procedures, the recursion is unwound by using a stack of VEC[INT] for
* the starting and ending points being checked. The stack is popped each
* time a call to ICSS1and2 returns either 0 or 1 breakpoint and pushed
* when it returns 2.
*
procedure ICSS x start end
type series 		x
type integer 		start end
*
option vect[int]  *breaks
option rect[int]  *subsamples
option switch     print    1
option string     title
option switch     trace    0
option real       signif   .05
option real       cv
option switch     robust   0
*
local vect[int]   lbreaks
local vect[int]   newbreaks
local vect[int]   t1 t2
local vect[int]   leftinterval rightinterval
local vect        temp
local integer     i k
local series      xc
local integer     startl endl
local report      ireport
local string      ltitle
local integer     notdone count newcount ifill
local real        cvlocal
*
if .not.%defined(x) {
   disp "Syntax: @icss(options)  X start end"
   return
}
if %defined(cv)
   compute cvlocal=cv
else
if signif==.05
   compute cvlocal=1.358
else
if signif==.01
   compute cvlocal=1.628
else
if signif==.10
   compute cvlocal=1.224
else {
   disp "@icss(SIGNIF=.01/.05/.10). If you want a different critical value, use CV option."
   return
}

inquire(series=x) startl<<start endl<<end
*
diff(center) x startl endl xc
*
dim t1(50) t2(50)
compute k=1
dim lbreaks(0)
compute t1(k)=startl,t2(k)=endl
until k==0 {
   comp newbreaks=ICSS1and2(xc,t1(k),t2(k),cvlocal,robust)
   compute lbreaks=%rlconcat(lbreaks,newbreaks)
   if trace.and.%rows(newbreaks)>0
      disp "ICSS-Adding" newbreaks "between" t1(k) "and" t2(k)
   *
   * If we find two breaks in the interval, add to the recursion stack
   * and search between them. Otherwise, pop the recursion stack.
   *
   if %rows(newbreaks)==2
      compute t1(k)=newbreaks(1)+1,t2(k)=newbreaks(2)
   else
      compute k=k-1
}
*
* This transfers to a real vector so we can sort the breakpoints
*
dim temp(%rows(lbreaks))
ewise temp(i)=lbreaks(i)
compute temp=%sort(temp)
ewise lbreaks(i)=fix(temp(i))
if trace
   disp "ICSS-Step 2 breaks" lbreaks
*
* Pruning step
*
loop
   compute notdone=0
   compute newbreaks=lbreaks
   compute count   =%rows(lbreaks)
   compute newcount=count
   dim leftinterval(count) rightinterval(count)
   ewise leftinterval(k)=%if(k==1,startl,lbreaks(k-1)+1)
   ewise rightinterval(k)=%if(k==count,endl,lbreaks(k+1))
   do k=1,count
      @IncTiaoBreak(cv=cvlocal) x leftinterval(k) rightinterval(k)
      if trace {
         if %maxent==0
            disp "ICSS-Step 3. Drop " lbreaks(k)
         else
            disp "ICSS-Step 3. Find " %maxent "between" leftinterval(k) $
                      "and" rightinterval(k)
      }
      if %maxent==0.or.%iabs(%maxent-lbreaks(k))>2
         compute notdone=1
      compute newbreaks(k)=%maxent
      if %maxent==0
         compute newcount=newcount-1
		if k<count
         compute leftinterval(k+1)=%if(%maxent==0,leftinterval(k),%maxent+1)
   end do k
   if notdone==0
      break
   dim lbreaks(newcount)
   compute ifill=0
   do k=1,count
      if newbreaks(k)<>0
         compute ifill=ifill+1,lbreaks(ifill)=newbreaks(k)
   end do k
end loop
*
if %defined(title)
   compute ltitle=title
else
   compute ltitle="ICSS Analysis of Series "+%l(x)
*
report(use=ireport,action=define,title=ltitle)
report(use=ireport,atrow=1,atcol=1,span) ltitle
*
if %defined(cv)
   report(use=ireport,atrow=2,atcol=1,span) "Using critical value "+%strval(cv,"*.###")
else
   report(use=ireport,atrow=2,atcol=1,span) "Using significance level "+%strval(signif,"*.##")
*
report(use=ireport,atrow=3,atcol=1,span) "Breakpoints found at entries"
do i=1,%rows(lbreaks)
   report(use=ireport,atcol=1,row=new) %datelabel(lbreaks(i))
end do i
if print
   report(use=ireport,action=show)
if %defined(breaks)
   compute breaks=lbreaks
if %defined(subsamples) {
   dim subsamples(2,%rows(lbreaks)+1)
   ewise subsamples(i,j)=%if(i==1,%if(j==1,startl,lbreaks(j-1)+1),$
                                  %if(j>%rows(lbreaks),endl,lbreaks(j)))
}
end

