*
* @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  Controls printing of the final regression and breakpoints
*   TESTS/[NOTESTS]  Controls printing of the table of BIC, LWZ and F tests
*                    for choosing the number of breaks.
*
* Supplementary Card:
*   list of regressors with the fixed regressors (if any) listed first
*
* Variables Defined
*   %RSS  - sum of squared residuals
*   %BETA - final 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
*  12/2009 Add TESTS option and improve output from PRINT
*  06/2010 Add calculation of confidence bands for break points
*  04/2011 Correct calculation of BIC.
*  05/2011 Correct output when NFIX is 1. Suppress output from intermediate
*          LINREG's if NFIX>0.
*  06/2011 Added PRINT=PRINT option to base LINREG
*  12/2011 Created version without TITLE options on REPORT for use with 7.2 and older
*****************************************************************************
*
* Run the regression with the input set of breaks
*
procedure BPRegression depvar startr endr breaks eqnfix eqnswitch
type series        depvar
type integer       startr  endr
type vect[integer] breaks
type equation      eqnfix eqnswitch
*
option switch		 print  1
*
local integer i j m nshift
local rect[series]  dz
local series        di
*
compute m=%rows(breaks)
compute nshift=%eqnsize(eqnswitch)
dim dz(nshift,m+1)
do i=1,m+1
   if i==1
      set di = t<=breaks(i)+startr-1
   else
   if i==m+1
      set di = t>breaks(i-1)+startr-1
   else
      set di = t>breaks(i-1)+startr-1.and.t<=breaks(i)+startr-1
   do j=1,nshift
      set dz(j,i) = di*%eqnxvector(eqnswitch,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(print=print,title="Bai-Perron Break Analysis") depvar startr endr
# %rlfromtable(%eqntable(eqnfix)) dz
end
***************************************
*
* Format up the table of statistics on the number of breaks
*
procedure BPFinishBreakTable breakreport
type report *breakreport
*
report(use=breakreport,action=format,atcol=2,width=10)
report(use=breakreport,action=format,atcol=3,tocol=4,picture="*.##")
report(use=breakreport,action=format,atcol=5,tocol=6,picture="*.##")
report(use=breakreport,action=format,atcol=3,tocol=3,tag=minimum,special=onestar,align=decimal)
report(use=breakreport,action=format,atcol=4,tocol=4,tag=minimum,special=onestar,align=decimal)
report(use=breakreport,action=show)
end
***************************************
*
* 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
option integer   	report      0
option report     *use
*
local vector 		x xv
local symm     	intervalrss
local integer 		s e time test lowerlim minbreak nobs
local integer 		nsegments nbreaks totalregs
local real 			norm baserss lastrss total minrss
local real 			u supf thisf sigsq
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>> segments
* (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
compute baserss=lastrss=bestrss(1,nobs)
*
* 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)
   compute totalregs=%nreg*(nbreaks+1)
   compute sigsq =minrss/(nobs-totalregs)
   compute supf=(baserss-minrss)/(%nreg*nbreaks*sigsq)
   compute thisf=(lastrss-minrss)/(%nreg*sigsq)
   compute lastrss=minrss
   if report {
      report(use=use,atrow=nbreaks+2,atcol=1) nbreaks minrss $
         log(minrss/nobs)+totalregs*log(nobs)/nobs $
         log(minrss/(nobs-totalregs))+totalregs*.299*log(nobs)^(2.1)/nobs $
         supf thisf
   }
end do nbreaks
compute %rss=minrss
end
*********************************************************************
*
* BPFuncG computes the CDF for the distribution of the break fraction
* This is a direct translation of Pierre Perron's Gauss code
*
function BPFuncG x bet alph b deld gam
type real BPFuncG x bet alph b deld gam
local real aa xb g
*
if x <= 0 {
   compute xb=bet*sqrt(abs(x))
   if (abs(xb) <= 30)
      compute g=-sqrt(-x/(2*%pi))*exp(x/8)-$
        (bet/alph)*exp(-alph*x)*%cdf(-bet*sqrt(abs(x)))+$
       ((2*bet*bet/alph)-2-x/2)*%cdf(-sqrt(abs(x))/2)
   else {
      compute aa=log(bet/alph)-alph*x-xb^2/2-log(sqrt(2*%pi))-log(xb);
      compute g=-sqrt(-x/(2*%pi))*exp(x/8)-exp(aa)*%cdf(-sqrt(abs(x))/2)+$
        ((2*bet*bet/alph)-2-x/2)*%cdf(-sqrt(abs(x))/2)
   }
}
else {
   compute xb=deld*sqrt(x)
   if (abs(xb) <= 30)
      compute g=1+(b/sqrt(2*%pi))*sqrt(x)*exp(-b*b*x/8)+$
        (b*deld/gam)*exp(gam*x)*%cdf(-deld*sqrt(x))+$
        (2-b*b*x/2-2*deld*deld/gam)*%cdf(-b*sqrt(x)/2)
   else {
      compute aa=log((b*deld/gam))+gam*x-xb^2/2-log(sqrt(2*%pi))-log(xb)
      compute g=1+(b/sqrt(2*%pi))*sqrt(x)*exp(-b*b*x/8)+exp(aa)+$
        (2-b*b*x/2-2*deld*deld/gam)*%cdf(-b*sqrt(x)/2)
	}
}
compute BPFuncG=g
end BPFuncG
***********************************************************************
*
* BPCVG solves for the critical value for a given significance level by
* bisection.
*
function BPCVG eta phi1s phi2s sig
type real BPCVG
type real eta phi1s phi2s
*
local real a gam b deld alph bet
local real sig upb lwb crit xx pval
local integer cct
*
compute a=phi1s/phi2s
compute gam=((phi2s/phi1s)+1)*eta/2
compute b=sqrt(phi1s*eta/phi2s)
compute deld=sqrt(phi2s*eta/phi1s)+b/2
compute alph=a*(1+a)/2
compute bet=(1+2*a)/2

compute upb=2000.0
compute lwb=-2000.0
compute crit=999999.0
do cct=1,100
   if abs(crit)<=.000001
      break;
   compute xx=lwb+(upb-lwb)/2
   compute pval=BPFuncG(xx,bet,alph,b,deld,gam)
   compute crit=pval-sig
   if crit <= 0
      compute lwb=xx
   else
      compute upb=xx
end do cct
compute BPCVG=xx
end
**************************************************************************
procedure BPBreakRanges startr endr eqnshift limits
type integer 		startr endr
type equation 		eqnshift
type rect[int]    *limits
*
option vect[int]  breaks
option integer    maxbreaks
option integer    nfix
option integer    nshift
option switch     robust       0
option switch     qhet         0
option switch	   omegahet     0
*
local rect	      betashift
local integer     i k
local vect[int]   ranges
local vect 	      delv
local symm        qmat qmat1
local symm        omega omega1
local real        phi1s phi2s
local real        eta lower upper a
*
compute betashift=%vectorect(%xsubvec(%beta,nfix+1,%nreg),nshift)
dim ranges(maxbreaks+2) limits(maxbreaks,2)
ewise ranges(i)=breaks(i-1)
compute ranges(1)=startr-1,ranges(maxbreaks+2)=endr
cmom startr endr
# %rlfromtable(%eqntable(eqnshift))
compute qmat=%cmom/%nobs
compute qmat1=qmat
if robust {
   mcov(lags=.75*(endr-startr+1)^.333,lwindow=quadratic) startr endr %resids
   # %rlfromtable(%eqntable(eqnshift))
   compute omega=%cmom/%nobs
   compute omega1=%cmom/%nobs
}

compute phi1s=%sigmasq
compute phi2s=%sigmasq
*
do k=1,maxbreaks
   compute delv =%xcol(betashift,k)-%xcol(betashift,k+1)
   if omegahet {
      if .not.robust {
      	sstats(mean) ranges(k)+1   ranges(k+1) %resids^2>>phi1s
      	sstats(mean) ranges(k+1)+1 ranges(k+2) %resids^2>>phi2s
      }
      else  {
         mcov(lags=.75*(ranges(k+1)-ranges(k))^.333,lwindow=quadratic) ranges(k)+1 ranges(k+1) %resids
         # %rlfromtable(%eqntable(eqnshift))
         compute omega=%cmom/%nobs
         mcov(lags=.75*(ranges(k+2)-ranges(k+1))^.333,lwindow=quadratic) ranges(k+1)+1 ranges(k+2) %resids
         # %rlfromtable(%eqntable(eqnshift))
         compute omega1=%cmom/%nobs
		}
	}
   if qhet {
      cmom ranges(k)+1 ranges(k+1) %resids
      # %rlfromtable(%eqntable(eqnshift))
      compute qmat=%cmom/%nobs
      cmom ranges(k+1)+1 ranges(k+2) %resids
      # %rlfromtable(%eqntable(eqnshift))
      compute qmat1=%cmom/%nobs
	}

   if robust {
      compute phi1s=%qform(omega,delv)/%qform(qmat,delv)
      compute phi2s=%qform(omega1,delv)/%qform(qmat1,delv)
   }
   compute a    =%qform(qmat,delv)/phi1s
   compute eta  =%qform(qmat1,delv)/%qform(qmat,delv)
   compute lower=BPCVG(eta,phi1s,phi2s,.025)
   compute upper=BPCVG(eta,phi1s,phi2s,.975)
   *
   * Round intervals away from the break point. Have to reverse the
   * signs since we're doing confidence intervals.
   *
   compute limits(k,1)=fix(breaks(k)-upper/a)
   compute limits(k,2)=fix(breaks(k)-lower/a+.999999)
end do k
end
**************************************************************************
*
* Main controlling procedure
*
**************************************************************************
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
option switch  tests     0
*
local integer 	it
local integer 	minsize
local integer  i j k m
local integer 	startr endr nobs totalregs
local integer  lower upper
local integer 	slots
local real 		baserss lastrss minrss oldrss sigsq supf thisf
local series 	ytilde
local integer 	nshift
local vec 		beta
local equation eqnfull eqnswitch eqnfix
local report   breakreport regreport
local vect[integer] breaks
local rect[integer] etable ptable
local rect[integer] limits
local integer  baserow
*
* Run the base regression (this gets the supplementary card in the
* process).
*
linreg(define=eqnfull,print=print) depvar start end
compute startr =%regstart()
compute endr   =%regend()
compute nobs   =endr-startr+1
compute baserss=lastrss=%rss
*
if tests {
	report(use=breakreport,action=define)
   report(use=breakreport,atrow=1,atcol=1,align=center) $
      "Breaks" "RSS" "BIC" "LWZ" "F(m)" "F(m|m-1)"
   report(use=breakreport,atrow=2,atcol=1) 0 %rss $
      log(%rss/nobs)+%nreg*log(nobs)/nobs $
      log(%rss/(nobs-%nreg))+%nreg*.299*log(nobs)^2.1/nobs
}
*
* Split the regressors into the shifting and fixed ones.
*
compute nshift=%nreg-nfix
compute etable=%eqntable(eqnfull)
dim ptable(2,nshift)
ewise ptable(i,j)=etable(i,j+nfix)
equation eqnswitch depvar
# %rlfromtable(ptable)
dim ptable(2,nfix)
ewise ptable(i,j)=etable(i,j)
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 slots =nobs
if nfix==0 {
	*
	* Compute the optimum allowing all regressors to break.
	*
   @BPDynProg(breaks=breaks,slots=slots,minsize=minsize,$
      maxbreaks=maxbreaks,equation=eqnfull,report=tests,use=breakreport) $
        depvar startr endr
}
else {
   if tests
      compute lower=1,upper=maxbreaks
   else
      compute lower=maxbreaks,upper=maxbreaks
   *
   * Loop over range of breaks needed
   *
   do m=lower,upper
      compute oldrss=%na
      *
      * Compute the optimum allowing all regressors to break
      *
      @BPDynProg(breaks=breaks,slots=slots,minsize=minsize,$
        maxbreaks=m,equation=eqnfull,report=0) depvar startr endr
      do it=1,iters
         *
         * Run the regression with the current set of breaks
         *
   		@BPRegression(noprint) depvar startr endr breaks eqnfix eqnswitch
         compute totalregs=%nreg
         compute minrss   =%rss
         *
         * 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=m,equation=eqnswitch,report=0) ytilde startr endr
         *
         * Keep doing this until there is no change.
         *
         if it>1.and.oldrss==%rss
            break
         compute oldrss=%rss
      end do it
      if tests {
         compute sigsq =minrss/(nobs-totalregs)
         compute supf  =(baserss-minrss)/(nshift*m*sigsq)
         compute thisf =(lastrss-minrss)/(nshift*sigsq)
         report(use=breakreport,atrow=m+2,atcol=1) m minrss $
           log(minrss/nobs)+totalregs*log(nobs)/nobs $
           log(minrss/(nobs-totalregs))+totalregs*.299*log(nobs)^(2.1)/nobs $
           supf thisf
      }
   end do m
}
if tests
   @BPFinishBreakTable breakreport

compute %nreg=(maxbreaks+1)*nshift+nfix
compute %%bic=log(%rss/nobs)+log(nobs)*%nreg/nobs
*
if nfix>0
   compute %beta=beta
if .not.print
   return
*
* Run the final regression
*
@BPRegression(print) depvar startr endr breaks eqnfix eqnswitch
*
* Normalize breaks to standard entry numbers
*
ewise breaks(i)=breaks(i)+startr-1
report(use=regreport,action=define)
report(use=regreport,atrow=1,atcol=1,span) "Bai-Perron Break Point Analysis"
report(use=regreport,atrow=2,atcol=1,span) "Dependent Variable "+%l(depvar)
compute baserow=4
*
* Display the fixed regressors
*
if nfix>0 {
   report(use=regreport,atrow=4,atcol=1,span) "Fixed Regressors"
   do i=1,nfix
      report(use=regreport,atrow=4+i,atcol=1) %eqnreglabels(eqnfull)(i)
      report(use=regreport,atrow=4+i,atcol=maxbreaks+2) %beta(i)
   end do i
   compute baserow=6+nfix
}
*
* Display the shifting regressors, one per line, with one column per partition
*
report(use=regreport,atrow=baserow,atcol=1,span) "Shifting Regressors"
do i=1,nshift
   report(use=regreport,atrow=baserow+i+1,atcol=1) %eqnreglabels(eqnfull)(nfix+i)
end do i
do k=1,maxbreaks+1
   if k<=maxbreaks
      report(use=regreport,atcol=k+1,atrow=baserow+1,align=center) "<="+%datelabel(breaks(k))
   do i=1,nshift
      report(use=regreport,atcol=k+1,atrow=baserow+i+1) %beta(nfix+(k-1)*nshift+i)
   end do i
end do k
compute baserow=baserow+nshift+3
*
@BPBreakRanges(robust,qhet,omegahet,maxbreaks=maxbreaks,nfix=nfix,nshift=nshift,breaks=breaks) startr endr eqnswitch limits
*
* Display the break points with 95% confidence bands
*
report(use=regreport,atrow=baserow,atcol=1) "Breakpoint" "Lower 95%" "Upper 95%"
do k=1,maxbreaks
   report(use=regreport,atrow=baserow+k,atcol=1) %datelabel(breaks(k)) %datelabel(limits(k,1)) %datelabel(limits(k,2))
end do k
report(use=regreport,action=show)
end BaiPerron


