*
* @BaiPerron( 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.
*
* Options:
*   MINSPAN = shortest distance between two breaks [number of regressors]
*   MAXBREAKS = maximum breaks allowed [2]
*   NFIX=number of regressors which are fixed over time[0]
*   ITERS=maximum number of iterations[20]. (Matters only if NFIX>0)
*   [PRINT]/NOPRINT
*
* Variables Defined
*   %RSS  - sum of squared residuals
*   %BETA - optimal fixed set of coefficients
*   %NREG - number of regressors
*   %NOBS - number of observations
*
* Revision Schedule:
*  08/2003 Written by Tom Doan, Estima
*  09/2003 Corrected problems when start<>entry 1. Thanks to Todd Clark, FRB-KC
*          for pointing this out.
*  09/2004 Improved output (now displays breakpoints in absolute terms, rather than
*          as number of entries from "start". Also, breakpoints are now displayed
*          as dates, rather than entry numbers whenever a CALENDAR is in use.
*  12/2006 NFIX and ITERS options added to allow for some coefficients to be fixed
*  12/2006 PRINT option added
*
* This is the dynamic programming routine for computing the best "maxbreaks" break
* points allowing all regressors to vary. It's not really designed to be used as
* a stand-alone procedure
*
procedure BPDynProg depvar startr endr
type series depvar
*
option vect[integer] *breaks
option integer slots
option integer minsize
option integer maxbreaks
option equation equation
*
local vector x xv
local symm intervalrss
local integer s e time test lowerlim minbreak nobs
local integer nsegments nbreaks
local real norm baserss total minrss
local real u
local rect bestrss
local rect[integer] lastbreak
*
dim intervalrss(slots,slots)
compute intervalrss=%const(%na)
*
* intervalrss is a lookup table of the residual sums of squares for the different
* segments. intervalrss(s,e) is the sum of squares on the interval [s,e].
*
* This uses a Kalman filter (written out as matrices) to fill in intervalrss.
*
compute nobs=endr-startr+1
do s=1,nobs-minsize+1
   linreg(equation=equation,noprint) depvar startr+(s-1) startr+(s-1)+(minsize-1)
   compute intervalrss(s,s+minsize-1)=%rss
   do e=s+minsize,nobs
      compute x = %eqnxvector(0,startr+e-1)
      compute xv = %xx*x
      compute norm=1.0+%dot(xv,x)
      compute %xx=%xx-(1.0/norm)*%outerxx(xv)
      compute u=depvar(startr+e-1)-%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 segment (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
   compute nsegments=nbreaks+1
*
*     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(nsegments-1,test)+intervalrss(test+1,time)
         if total<minrss
            compute minrss=total,minbreak=test
      end do test
      compute bestrss(nsegments,time)=minrss,lastbreak(nsegments,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 s
   compute minrss=bestrss(nsegments,nobs)
end do nbreaks
compute %rss=minrss
end
*
procedure BaiPerron depvar start end
type series depvar
type integer start end
*
option integer minspan
option integer maxbreaks 2
option integer nfix      0
option integer iters     20
option switch  print     1
*
local integer it
local integer minsize
local integer startr endr nobs n
local integer slots
local real baserss oldrss
local vect[integer] breaks
local series di ytilde
local rect[series] dz
local equation eqnfull eqnpart eqnfix
local rect[integer] etable ptable
local integer nshift
local vec beta
*
* Run the base regression (this gets the supplementary card
* in the process).
*
linreg(define=eqnfull) depvar start end
compute baserss=%rss
compute nshift=%nreg-nfix
if nfix>0 {
   compute etable=%eqntable(eqnfull)
   dim ptable(2,nshift)
   do i=1,nshift
      compute ptable(1,i)=etable(1,i+nfix),ptable(2,i)=etable(2,i+nfix)
   end do i
   equation eqnpart depvar
   # %rlfromtable(ptable)
   dim ptable(2,nfix)
   do i=1,nfix
      compute ptable(1,i)=etable(1,i),ptable(2,i)=etable(2,i)
   end do i
   equation eqnfix depvar
   # %rlfromtable(ptable)
}
*
* Get the minimum span between breaks. Make sure it's
* at least as big as the number of variables permitted to shift
*
if .not.%defined(minspan).or.minspan<nshift
   compute minsize=nshift
else
   compute minsize=minspan

compute startr=%regstart(),endr=%regend()
compute nobs  =endr-startr+1
compute slots =nobs
*
* Compute the optimum allowing all regressors to break.
*
@BPDynProg(breaks=breaks,slots=slots,minsize=minsize,maxbreaks=maxbreaks,equation=eqnfull) depvar startr endr
*
* If there are fixed regressors, iterate over the Bai-Perron procedure
*
if nfix>0 {
   dim dz(maxbreaks+1,nshift)
   compute oldrss=%na
   do it=1,iters
      do i=1,maxbreaks+1
         if i==1
            set di = t<=breaks(i)
         else
         if i==maxbreaks+1
            set di = t>breaks(i-1)
         else
            set di = t>breaks(i-1).and.t<=breaks(i)
         do j=1,nshift
            set dz(i,j) = di*%eqnxvector(eqnpart,t)(j)
         end do j
      end do i
      *
      * Run the regression on the fixed regressors + the segmented variables over the
      * current choices for the break points
      *
      linreg(noprint) depvar startr endr
      # %rlfromtable(ptable) dz
      *
      * Compute the deviations of y from the fixed regressors
      *
      compute beta=%xsubvec(%beta,1,nfix)
      set ytilde = depvar-%dot(beta,%eqnxvector(eqnfix,t))
      *
      * Find the best breaks over the varying regressors using the deviation as the
      * dependent variable
      *
      @BPDynProg(breaks=breaks,slots=slots,minsize=minsize,maxbreaks=maxbreaks,equation=eqnpart) ytilde startr endr
      disp %rss oldrss
      if it>1.and.oldrss==%rss
         break
      compute oldrss=%rss
   end do it
}
compute %nreg=(maxbreaks+1)*nshift+nfix
compute %%bic=log(%rss/nobs)+log(nobs)*%nreg/nobs
if nfix>0
   compute %beta=beta
if print {
   disp "Best" maxbreaks "break points"
   do n=1,maxbreaks
      dis %datelabel(breaks(n)+startr-1)
   end do n
   *
   disp "Sum of Squared Residuals" @27 %rss
   disp "BIC" @27 %%bic
}
end BaiPerron

