*
* @ICSS  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.
*
* Parameters:
*   x = series to analyze
*   start end = range to analyze. By default, the defined range of x.
*
* Revision Schedule:
*   03/2005 Written by Tom Doan, Estima
*   04/2008 Stripped out of program file.
*
*******************************************************************************************
*
* 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 asymptotic .95 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
local series ck dk adk
local real nobs m
*
set ck start end = x**2
acc ck start end
compute nobs=(end-start+1)
set dk start end = ck/ck(end)-(t-start+1)/nobs
*graph
*# dk start end
set adk start end = abs(dk)
ext(noprint) adk start end
compute m=sqrt(.5*nobs)*%maximum
if m>1.358
   ;
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
type vect[int] ICSS1and2
type series x
type integer start end
*
local integer t1 t2 kstar kfirst klast
*
compute t1=start,t2=end
@IncTiaoBreak x start end
if %maxent==0 {
   dim ICSS1and2(0)
   return
}
compute t2=kstar=%maxent
loop
   @IncTiaoBreak x t1 t2
   if %maxent>0
      compute t2=%maxent
   else
      break
end loop
compute kfirst=t2
compute t1=kstar+1,t2=end
loop
   @IncTiaoBreak 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. Note that this does not
* do the "pruning" step (step 3), just the complete enumeration from steps 1 and
* 2. 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
*
local vect[int] breaks
local vect[int] newbreaks
local vect[int] t1 t2
local vect temp
local integer i k
local series xc
local integer startl endl
*
inquire(series=x) startl<<start endl<<end
*
diff(center) x startl endl xc
*
dim t1(10) t2(10)
compute k=1
dim breaks(0)
compute t1(k)=startl,t2(k)=endl
until k==0 {
   comp newbreaks=ICSS1and2(xc,t1(k),t2(k))
   compute breaks=%rlconcat(breaks,newbreaks)
   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(breaks))
ewise temp(i)=breaks(i)
compute temp=%sort(temp)
ewise breaks(i)=fix(temp(i))
disp 'Breakpoints found at entries'
disp breaks
end


