*
* @MultipleBreaks( options ) depvar start end
* # list of regressors
*
* This does multiple structural change analysis as described in Bai and Perron,
* "Computation and Analysis of Multiple Structural Change Models", Journal of
* Applied Econometrics, 2003, pp 1-22, but using a threshold variable other than
* time.
*
* Syntax:
*   depvar     dependent variable
*   start end  range to estimate (by default, the common range of all variables)
*
* Options:
*   MINSPAN   = smallest number of data points in a partition [number of regressors]
*   MAXBREAKS = maximum breaks allowed [2]
*   THRESH=threshold series. If omitted, time is used. If you're using a lagged
*    value, you need to create the lag as a separate series.
*
* Revision Schedule:
*  07/2005 Written by Tom Doan, Estima based upon a combination of the procedures
*          @BaiPerron and @Threshold
*
procedure MultipleBreaks depvar start end
type series depvar
type integer start end
*
option integer minspan
option integer maxbreaks 2
option series  thresh
*
local integer minsize slots
local integer startr endr nobs n
local symm intervalrss
local real u
local integer s e eadd lowerlim
local integer time test nbreaks
local vector x xv
local real norm baserss minrss total
local rect bestrss
local rect[integer] lastbreak
local integer minbreak
local vect[integer] breaks
local series threshold copy ix smpls
local label breakvar
*
* Run the base regression (this gets the supplementary card
* in the process).
*
linreg depvar start end
compute baserss=%rss
*
* Get the minimum span between breaks. Make sure it's
* at least as big as %nreg
*
if .not.%defined(minspan).or.minspan<%nreg
   compute minsize=%nreg
else
   compute minsize=minspan

compute startr=%regstart(),endr=%regend()
compute nobs  =endr-startr+1
compute slots =nobs
*
* If user provided a threshold variable, use it; otherwise
* use time.
*
if %defined(thresh) {
   set threshold startr endr = thresh
   compute breakvar=%l(thresh)
}
else {
   set threshold  startr endr = t
   compute breakvar='Time'
}
*
*  Make a copy of the threshold variable and sort it, along with a series of
*  entry numbers. (We can't just use RANK on ORDER because it will map ties
*  to an average rank).
*
set copy startr endr = threshold
set ix   startr endr = t
*
order copy startr endr ix
*
* <<intervalrss>> is the lookup table of the residual sums of squares for the
* different segments
*
dim intervalrss(slots,slots)
compute intervalrss=%const(%na)
*
* This uses a Kalman filter (written out as matrices) to fill in <<intervalrss>>
*
do s=1,nobs-minsize+1
*
*     Do smallest regression beginning <<s>> positions into the data (in
*     threshold variable order)
*
   set smpls startr endr = 0.0
   do e=s,s+minsize-1
      compute smpls(fix(ix(startr+e-1)))=1.0
   end do e
   linreg(lastreg,noprint,smpl=smpls) depvar startr endr
   compute intervalrss(s,s+minsize-1)=%rss
   do e=s+minsize,nobs
      compute eadd = fix(ix(startr+e-1))
      compute x = %eqnxvector(0,eadd)
      compute xv = %xx*x
      compute norm=1.0+%dot(xv,x)
      compute %xx=%xx-(1.0/norm)*%outerxx(xv)
      compute u=depvar(eadd)-%dot(x,%beta)
      compute %rss=%rss+u**2/norm
      compute %beta=%beta+(u/norm)*xv
      compute intervalrss(s,e)=%rss
   end do e
end do s
*
* bestrss(i,j) shows the best rss that can be created from i-1 breaks
* in observations through j. We initialize this with the straightforward
* case of i=1 (or 0 breaks).
*
dim bestrss(maxbreaks+1,nobs) lastbreak(maxbreaks+1,nobs)
compute bestrss=%const(%na)
do time=minsize,nobs
   compute bestrss(1,time)=intervalrss(1,time)
end do time
*
* Now run through the number of breaks from 1 to the maximum. At each stage,
* figure out the best sum of squared residuals that can be obtained with *nbreaks*
* break points between 1 and time, for all values of time for which that is
* sensible. This is done using the dynamic programming method described in Bai and
* Perron. For each combination of # of breaks and ending time period, we save the
* sum of squared residuals and the final break point. The latter is used to
* rebuild the optimal set of break points.
*
do nbreaks=1,maxbreaks
*
*     When we're at maxbreaks, we don't need to compute bestrss for
*     anything but the full sample. Otherwise, we need to do all
*     points from nbreaks*minsize to nobs
*
   if nbreaks==maxbreaks
      compute lowerlim=nobs
   else
      compute lowerlim=nbreaks*minsize
   do time=nobs,lowerlim,-1
      compute minrss=baserss
      do test=time-minsize,nbreaks*minsize,-1
         compute total=bestrss(nbreaks,test)+intervalrss(test+1,time)
         if total<minrss
            compute minrss=total,minbreak=test
      end do test
      compute bestrss(nbreaks+1,time)=minrss,lastbreak(nbreaks+1,time)=minbreak
   end do time

   dim breaks(nbreaks)
   compute test=nobs
   do s=nbreaks+1,2,-1
      compute breaks(s-1)=test=lastbreak(s,test)
   end do recurse
   compute minrss=bestrss(nbreaks+1,nobs)
   disp 'Best' nbreaks 'breakpoints'
   do n=1,nbreaks
     dis copy(breaks(n)+startr-1)
   end do n
   disp
   disp 'Sum of Squared Residuals' @27 bestrss(nbreaks+1,nobs)
   disp 'BIC' @27 log(minrss/nobs)+log(nobs)*(nbreaks+1)*%nreg/nobs
   disp
   disp
end do nbreaks
end MultipleBreaks


