*
* @ThreshTest( options ) depvar start end
* # list of explanatory variables
*
* This checks for a threshold break in a linear regression, assuming
* homoscedastic errors.
*
* Options:
*   SMPL=SMPL series for entries to skip
*   THRESH=threshold series. If omitted, time is used
*   D=lag delay for threshold [0]
*   TRIM=fraction of values to skip at both ends [.15]
*   GRAPH/[NOGRAPH] If chosen, graphs the F-tests against the
*     threshold variables
*   NREPS=number of bootstrap replications [0]
*     If you provide a positive value, this computes an approximate
*     significance level by bootstrapping with random Normal draws.
*
*   TITLE=title for report ["Test of Threshold Effect"]
*   [PRINT]/NOPRINT
*
* Variables Defined:
*   %CDSTAT      Maximum test statistic
*   %SIGNIF      Bootstrapped p-value
*   %%BREAKVALUE Threshold value with highest F
*
* Revision Schedule:
* 07/2003 Written by Tom Doan, Estima. Based somewhat on Bruce Hansen's
*   Gauss procedure thr_test. The main difference on output is that this
*   computes F tests (the same as you could get on doing a Chow test
*   breaking at the indicated value). Internally, it works off an index
*   into the naturally ordered data set, so we don't have to make a copy
*   of the data set if we're using lags.
* 08/2010 Added PRINT and TITLE options. Switched to use of REPORT.
* 06/2012 Added definition of %%BREAKVALUE. Put it on the graph.
* 11/2014 Add check for TRIM too small for size of model/data set
* 11/2017 Add D option
*
procedure ThreshTest depvar start end
type series 	depvar
type integer 	start end
*
option series     thresh
option series     smpl
option integer    nreps    0
option real       trim     .15
option switch     graph    0
option switch     print    1
option string     title
option integer    d        0
*
local integer 	startl endl nobs nb ne time retime reps i
local vect 		fboots signif
local series 	resids smpll copy ix ssq threshold ftests e u flatspot
local vect 		xu x
local symm 		xx xxxx
local real 		ssqmax ssqtest sigma fmax breakvalue
local label 	breakvar
local report 	treport
local string	ltitle
local vect[symm] xxs
*
* Run the regression, which gets the regressors at the same time.
*
linreg(smpl=smpl,noprint) depvar start end resids
compute startl=%regstart(),endl=%regend()
set smpll startl endl = %valid(resids)

if trim<0.or.trim>=.5 {
   disp "### @THRESHTEST. Value of TRIM is " trim ". Must be between 0 and .5"
   return
}
if trim*%nobs<%nreg {
  disp "### @THRESHTEST. Value of TRIM is" trim
  disp "Too small for model with " %nreg "regressors and " %nobs " observations"
  return
}
*
* If user provided a threshold variable, use it; otherwise use time.
*
if %defined(thresh) {
   set threshold startl endl = thresh{d}
   if d
      compute breakvar=%l(thresh)+"{"+d+"}"
   else
      compute breakvar=%l(thresh)
}
else {
   set threshold  startl endl = t
   compute breakvar="Time"
}
*
* Make a copy of the threshold variable, patch it to take missing values
* in the regression into account 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 startl endl = %if(smpll,threshold,%na)
set ix   startl endl = t
*
order copy startl endl ix
*
* flatspot is 1 at locations where we don't want to allow a break.
* Initially it's at all positions where there are duplicate threshold
* values - it will be one until we hit the end of the span. flatspot
* will also be zero for the trimmed part of the sample.
*
set flatspot = (t<endl.and.copy==copy{-1})
compute nobs = %nobs
*
* Figure out the starting and ending entries at which we're going to
* allow breaks.
*
compute nb = startl+fix(trim*nobs),ne=startl+fix((1-trim)*nobs)
set flatspot startl nb = 1
*
* We only need to store the xx matrices if we're doing bootstrapping.
*
if nreps
   dim xxs(ne-nb+1)
*
* This computes sequentially the difference in the sum of squares
* between OLS and OLS with a sample split at "time" (according to the
* sort order on the threshold variable). This difference can be computed
* as u'x (x'x-x'x (X'X)**-1 x'x)**-1 x'u where X'X is the full sample
* X'X, x'x is the same for the partial sample, u are the full sample
* residuals and x'u is computed over the partial sample.
*
dim xu(%nreg) xx(%nreg,%nreg)
compute xx=%const(0.0),xu=%const(0.0)
compute ssqmax=0.0
set ssq nb ne = %na
do time=startl,ne
   compute retime=fix(ix(time))
   compute x=%eqnxvector(0,retime)
   compute xu=xu+x*resids(retime)
   compute xx=xx+%outerxx(x)
   if flatspot(time)
      next
   compute xxxx=inv(xx-xx*%xx*xx)
   *
   * xxxx is an array which will be needed for bootstrapping and which
   * is fairly time-consuming to compute.
   *
   if nreps>0
   	compute xxs(time-nb+1)=xxxx
   set ssq time time = %qform(xxxx,xu)
   *
   * See if this makes the largest break.
   *
   if ssq(time)>ssqmax
      compute ssqmax=ssq(time),breakvalue=threshold(retime)
end do time
compute %%breakvalue=breakvalue
*
*  Convert everything to F statistics
*
compute sigma = (%rss-ssqmax)/(%nobs-2*%nreg)
compute fmax  = (ssqmax/%nreg)/sigma
set ftests nb ne = (ssq/%nreg)/sigma
*
* If requested, do a graph
*
if graph {
   stats(noprint,fract) copy nb ne
   spgraph
   scatter(style=lines,hlabel=breakvar,vlabel="F-Tests",header="Threshold Tests",$
     hgrid=breakvalue)
   # copy ftests nb ne
   if breakvalue>%median
      grtext(align=right,x=breakvalue,y=ssqmax,valign=bottom) %strval(breakvalue,"*.##")
   else
      grtext(align=left,x=breakvalue,y=ssqmax,valign=bottom) %strval(breakvalue,"*.##")
   spgraph(done)
}
*
dim fboots(nreps) signif(nreps)
do reps=1,nreps
   *
   * Convert independent draws to the same correlation pattern as the
   * sample residuals.
   *
   set e startl endl = %ran(1.0)
   linreg(lastreg,smpl=smpll,noprint) e / u
   *
   * Given that, same calculation as above, but with xxxx taken from the
   * xxs cache.
   *
   compute xu=%const(0.0)
   compute ssqmax=0.0
   do time=startl,ne
      compute retime=fix(ix(time))
      compute x=%eqnxvector(0,retime)
      compute xu=xu+x*u(retime)
      if flatspot(time)
         next
      compute ssqtest = %qform(xxs(time-nb+1),xu)
      if ssqtest>ssqmax
         compute ssqmax=ssqtest
   end do time
   compute sigma = (%rss-ssqmax)/(%nobs-2*%nreg)
   compute fboots(reps)=ssqmax/(%nreg*sigma)
end do reps
*
* Figure out how many of the bootstrapped F's exceed the value of fmax.
*
ewise signif(i)=(fboots(i)>=fmax)
compute %signif=%sum(signif)/nreps

compute %cdstat=fmax

if %defined(title)
   compute ltitle=title
else
   compute ltitle="Test of Threshold Effect"

report(use=treport,action=define,title=ltitle)
report(use=treport,atrow=1,atcol=1,span) "Test of Null of No Threshold Against Alternative of Threshold"
report(use=treport,atrow=2,atcol=1,span) "Under Maintained Assumption of Homoskedastic Errors"
report(use=treport,atrow=4,atcol=1)  "Dependent Variable" 	%l(depvar)
report(use=treport,atrow=5,atcol=1)  "Threshold Variable" 	breakvar
report(use=treport,atrow=6,atcol=1)  "Observations"			%nobs
report(use=treport,atrot=7,atcol=1)  "From"				   	%datelabel(startl)
report(use=treport,atrow=8,atcol=1)  "To"				   		%datelabel(endl)
report(use=treport,atrow=9,atcol=1)  "Maximum F-Test" 		fmax
report(use=treport,atrow=10,atcol=1) "Achieved at" 			breakvalue
if nreps>0 {
   report(use=treport,atrow=11,atcol=1) "Bootstrap Replications " nreps
   report(use=treport,atrow=12,atcol=1) "Bootstrap P-Value"       %signif
}
report(use=treport,action=format,atcol=2,picture="*.####")
if print
   report(use=treport,action=show)
end
