*
* @APBreakTest(options)  depvar start end
* # list of regressors
*
* APBreakTest performs Andrews-Ploberger and Andrews-Quandt structural break tests
* for a linear regression, with p-values using Hansen's approximations. These
* are used to test for a single structural break at an unknown point within the
* sample. A series of LM statistics are generated for breaks at each of the points
* in the middle range of the data set. The Andrews-Quandt test uses as the test
* statistic the maximum of the LM statistics, while Andrew-Ploberger uses an
* exponentially weighted average. These both have highly non-standard distributions.
* Asymptotic p-values are computed.
*
* As currently written, this examines break points in the 15%-85% range of the
* entries of the data set. A different "pi" value than .15 would require a
* different table of entries for SupWTheta and ExpWTheta.
*
* References:
*   Andrews, Donald W. K. and Werner Ploberger, "Optimal Tests When a Nuisance
*    Parameter is Present Only Under the Alternative", Econometrica, 1994, pp 1383-1414
*
*   Hansen, Bruce E., "Approximate Asymptotic P-Values for Structural Change Tests",
*    Journal of Business and Economic Statistics, 1997, pp 60-67.
*
* Option:
*   ROBUST/[NOROBUST]  If ROBUST, heteroscedasticity consistent LM tests are used
*
* Revision Schedule:
*   05/2003 Written by Tom Doan, Estima, based upon a program written by Todd Clark,
*     Federal Reserve Bank, Kansas City, and a Gauss program written by Bruce Hansen,
*     University of Wisconsin.
*   08/2003 Revised to allow for missing data and add SMPL option
*   01/2004 Updated to use v6 features
*
procedure APBreakTest depvar start end
type series depvar
type integer start end
*
option real   pi     .15
option switch robust 0
option series smpl
*
local series resids resids2 dummy tempser smpll
local integer startl endl time piStart piEnd i
local integer count tstart tend
local vect x
*
* To hold full sample OLS information
*
local integer olsdf
local real olsrss olsvar
local symm olsxx olscv
*
* For single coefficient case
*
local real dss des dsees sss rssdiff lm
local vect dxs dxes dxees
*
* For full vector case
*
local symm dxx dxeex ss
local vect dxe
*
* For test bookkeeping
*
local vect[labels] brkslabels
local vect[int] dfs
local vect[series] wstat rss
local vector supwtests supw_asypval expwtests expw_asypval
local vect[int] brkdates
*
local integer nx ntests
local real zed
*
*** pars for computing Hansen's asymptotic p-values
*** note:  values are specific to pi0=.15
*** coefficient values taken from Hansen's Gauss code (Pv_sup.prg, Pv_exp.prg)
*** this allows for up to 25 coefficients to be tested for a break
*
fixed rect SupWtheta(25,3)
  -0.98828874    1.01604712    2.96752918
  -1.65273586    1.05932633    4.69340428
  -2.04820499    1.13491785    6.75208304
  -2.52186324    1.10727065    7.95561422
  -3.46464220    1.06781143    8.27206632
  -4.04745714    1.08339951    9.51466537
  -4.41913790    1.10217153   11.01503763
  -5.36063931    1.07646817   11.25922919
  -5.43335672    1.10431373   13.09302596
  -6.46819929    1.06017100   12.81597962
  -6.78917650    1.03551668   13.48939975
  -7.79706943    1.02343400   13.63889444
  -7.93078073    1.06995710   15.86478742
  -8.54188927    1.04729101   16.11546709
  -9.05194333    1.05477016   17.22576198
  -9.13026570    1.08679108   19.32503948
 -10.45110576    1.04570811   18.25837821
 -10.62998450    1.05028481   19.53064308
 -12.13704880    0.90008482   14.94170023
 -12.13536344    0.97206863   18.31406331
 -12.57665664    1.03475703   21.09789423
 -12.76806289    1.07068440   23.38510717
 -13.05498712    1.09274429   25.19280476
 -14.51359909    0.97639956   20.94992826
 -14.16209122    1.04603251   25.04274427
fixed rect ExpWtheta(25,5)
  -0.41546194    3.75135495   -0.65483135    0.10063927    1.68830619
  -0.76804658    2.66074495   -0.06635619    0.00000000    2.57852729
  -1.16062736    2.91102163   -0.08409545    0.00000000    4.16250418
  -1.46819315    2.25606670    0.00000000    0.00000000    4.41017364
  -1.94762736    2.19579911    0.00000000    0.00000000    5.17227554
  -2.45289663    2.12836984    0.00000000    0.00000000    5.82436481
  -2.73487377    2.19317991    0.00000000    0.00000000    7.12126224
  -2.91474255    2.22180407    0.00000000    0.00000000    8.46556736
  -3.58424345    2.19826545    0.00000000    0.00000000    9.06067942
  -4.03467984    2.11158946    0.00000000    0.00000000    9.39679396
  -4.66068165    2.14833682    0.00000000    0.00000000   10.30022638
  -5.20162379    2.07125608    0.00000000    0.00000000   10.49479632
  -5.49275767    2.04442714    0.00000000    0.00000000   11.28869335
  -5.55515654    2.16577246    0.00000000    0.00000000   13.45347527
  -6.66123413    2.10863845    0.00000000    0.00000000   13.12802116
  -6.92233210    2.10938011    0.00000000    0.00000000   14.18991387
  -7.58056409    2.07988244    0.00000000    0.00000000   14.43289540
  -8.15325008    1.93614845    0.00000000    0.00000000   13.50372176
  -8.77892414    1.92603583    0.00000000    0.00000000   13.97185104
  -8.96368910    2.02143495    0.00000000    0.00000000   16.05934545
  -8.90257536    2.13592724    0.00000000    0.00000000   18.77330812
  -9.91956050    2.06310660    0.00000000    0.00000000   18.04892040
  -9.62431215    2.19525793    0.00000000    0.00000000   21.41582378
 -10.95103780    1.95541117    0.00000000    0.00000000   17.89956485
 -12.17923782    1.88294675    0.00000000    0.00000000   16.70485644
*
* Do the base linear regression, save the residuals. Fetch the
* range actually used.
*
linreg(print=.not.robust,smpl=smpl) depvar start end resids
inquire(lastreg) startl endl
set smpll startl endl = %valid(resids)
*
* Figure out the start and end entries for break tests allowing for
* possible missing values
*
compute count=0,tstart=fix(pi*(%nobs-1)),tend=%nobs-1-tstart
do time=startl,endl
   if count==tstart
      compute piStart=time
   if count==tend {
      compute piEnd  =time
      break
   }
   compute count=count+fix(smpll(time))
end do time
*
compute olsrss=%rss,olsvar=%seesq,olsdf=%ndf,olsxx=%xx
compute nx=%nreg,ntests=%nreg+2
if robust {
   linreg(lastreg,print,robusterrors,smpl=smpll) depvar start end
   compute olscv=%xx
}
dim wstat(ntests) rss(ntests)
dim dfs(ntests) brkslabels(ntests) brkdates(ntests)
dim supwtests(ntests) supw_asypval(ntests) expwtests(ntests) expw_asypval(ntests)
*
* Individual coefficients
*
dim dxs(nx) dxes(nx) dxees(nx)
do i=1,nx
   compute brkslabels(i)=%eqnreglabels(0)(i)
   compute dfs(i)=1
   set wstat(i) piStart piEnd = 0.0
   set rss(i)   piStart piEnd = 0.0
   compute dxs=%const(0.0),dxes=%const(0.0),dxees=%const(0.0)
   compute dss=0.0,des=0.0,dsees=0.0
   do time=startl,piEnd
      if .not.smpll(time)
         next
      compute x    =%eqnxvector(0,time)
      compute dss  =dss  +x(i)**2
      compute dxs  =dxs  +x(i)*x
      compute des  =des  +x(i)*resids(time)
      compute dxes =dxes +x(i)*resids(time)*x
      compute dsees=dsees+(x(i)*resids(time))**2
      compute dxees=dxees+x(i)*resids(time)**2*x
      if time<piStart
         next
      compute sss=dss-%qform(olsxx,dxs)
      compute rssdiff=(1.0/sss)*des**2
      compute rss(i)(time)=olsrss-rssdiff
      if robust {
         compute sss=dsees-2.0*%dot(dxs,olsxx*dxees)+%qform(olscv,dxs)
         compute lm =(1.0/sss)*des**2
      }
      else
         compute lm=rssdiff/olsvar

      compute wstat(i)(time)=lm
   end do time
end do i
*
* Full coefficient vector
*
dim dxx(nx,nx) dxeex(nx,nx) dxe(nx)
compute dxx=%const(0.0),dxeex=%const(0.0),dxe=%const(0.0)
compute brkslabels(nx+1)='All Coeffs'
compute dfs(nx+1)=nx
set wstat(nx+1) piStart piEnd = 0.0
set rss(nx+1)   piStart piEnd = 0.0
do time=startl,piEnd
   if .not.smpll(time)
      next
   compute x=%eqnxvector(0,time)
   compute dxx  =dxx  +%outerxx(x)
   compute dxeex=dxeex+%outerxx(x*resids(time))
   compute dxe  =dxe  +x*resids(time)
   if time<piStart
      next
   compute ss=(dxx-dxx*olsxx*dxx)
   compute rssdiff=%qform(inv(ss),dxe)
   compute rss(nx+1)(time)=olsrss-rssdiff
   if robust {
      compute ss=dxeex-dxx*olsxx*dxeex-dxeex*olsxx*dxx+dxx*olscv*dxx
      compute lm=%qform(inv(ss),dxe)
   }
   else
      compute lm=rssdiff/olsvar
   compute wstat(nx+1)(time)=lm
end do time
*
* Variance
*
set resids2 startl endl = resids**2
compute brkslabels(nx+2)='Resid Var'
compute dfs(nx+2)=1
set wstat(nx+2) piStart piEnd = 0.0
set rss(nx+2)   piStart piEnd = 0.0
*
do time=piStart,piEnd
   if .not.smpll(time)
      next
   set dummy = t<=time
   linreg(noprint) resids2 startl endl
   # constant dummy
   compute wstat(nx+2)(time)=%tstats(2)
   compute rss(nx+2)(time)=%rss
end do time
*
do i=1,ntests
   *** Andrews
   extremum(noprint) wstat(i) piStart piEnd
   comp supwtests(i) = %maximum
   comp zed = %max(0.0,SupWtheta(dfs(i),1) + SupWtheta(dfs(i),2)*supwtests(i))
   comp supw_asypval(i) = %chisqr(zed,SupWtheta(dfs(i),3))
   *** Andrews date
   extremum(noprint) rss(i) piStart piEnd
   comp brkdates(i) = %minent
   *** Andrews-Ploberger
   set tempser piStart piEnd = exp(.5*wstat(i))
   stats(noprint) tempser piStart piEnd
   comp expwtests(i) = log(%mean)
   comp zed = %max(0.0,ExpWtheta(dfs(i),1) + ExpWtheta(dfs(i),2)*expwtests(i) + ExpWtheta(dfs(i),3)*expwtests(i)**2. + ExpWtheta(dfs(i),4)*expwtests(i)**3.)
   comp expw_asypval(i) = %chisqr(zed,ExpWtheta(dfs(i),5))
end do i
******** displaying results

report(action=define)
report(atrow=1,atcol=2,span) 'Andrews-Quandt'
report(atrow=1,atcol=5,span) 'Andrews-Ploberger'
report(atrow=2,atcol=2) 'Test' 'P-Val' 'Date' 'Test' 'P-Val'
do i=1,ntests
   report(row=new,atcol=1) brkslabels(i) supwtests(i) supw_asypval(i) %datelabel(brkdates(i)) expwtests(i) expw_asypval(i)
end do i
report(atcol=2,tocol=6,width=9)
report(atcol=3,tocol=3,picture='#.###')
report(atcol=6,tocol=6,picture='#.###')
report(action=show)
end APBreakTest



