*
* @PanelThresh( options ) depvar start end
* # list of regressors (omit CONSTANT)
*
* Tests a fixed effects regression for one or two breaks in a single
* regressor, with breaks determined by the values of another variable.
*
* The regression (with two breaks) takes the form:
*
* y(i,t)=a(i)+x(t,i)b+z(t,i)*(q(t,i)<=gamma1)+z(t,i)*(q(t,i)<=gamma2)
*
* This requires a separate fixed effects regression for each test
* setting for the gammas. If q is a continuous variable, a complete
* search over all possible values in a panel data set could be very time
* consuming. As a result, you can examine a reduced number of quantiles
* of the values of q. It also only does an exhaustive search for the
* first break, then examines the second break given the first.
*
* Options:
*   THRESHVAR=series that determines the breaks [required]
*   BREAKVAR =series that is dummied out at the break
*   QN=number of quantiles of the threshold variable to examine
*      default [max(400,actual values)]
*   T1=trimming fraction on ends for first threshold [.01]
*   [PRINT]/NOPRINT
*   TITLE=title for report ["Panel Threshold Analysis, Threshold Variable xxx"]
*   SMPL=Dummy series with 0's in entries to skip [include all]
*
*   GRAPH/[NOGRAPH]
*   HEADER=header for graph [same as title]
*   FOOTER=footer for graph [not used]
*   SUBHEADER=subheader for graph [not used]
*     Draws an X-Y graph of the F-tests for a single break vs the break values
*
* Variables Defined:
*   %%rss1     Minimum sum of squared residuals for one break
*   %%ftest1_0 Maximal F-test for one break
*   %%gamma1_1 Best single break value
*   %%rss2     (Nearly) minimum sum of squared residuals for two breaks
*   %%ftest2_1 F-test for two breaks given one break
*   %%gamma1_2 Break point in two break model
*   %%gamma2_2 Other break point in two break model
*
* Reference:
* Bruce Hansen(1999), "Threshold effects in non dynamic panels:
* estimation, testing and inference", Journal of Econometrics, vol 93,
* pp 345-368.
*
* The original Gauss program for this had an error, (unnecessarily)
* dropping a data point from the end of each individual.
*
* Revision Schedule:
*  08/2011 Written by Tom Doan, Estima
*  05/2012 Revised to add PRINT and TITLE options
*  06/2012 Revised to add SMPL, GRAPH and related options
*
*********************************************************************
procedure MakeGammas gammas
type vector *gammas
option integer obs       %nobs
option integer quantiles obs
option real    trim      .01
option vector  sorted
*
local real    	range base
local integer 	i
*
compute range=obs*(1-2*trim),base=obs*trim
dim gammas(quantiles)
ewise gammas(i)=sorted(fix(range*(i-1)/(quantiles-1)+base))
end
*********************************************************************
procedure PanelThresh depvar start end
option integer qn 400
option real    t1 .01
option real    t2 .01
option real    t3 .05
option series  threshvar
option series  breakvar
option switch  print  	1
option integer nbreaks 	2
option string  title
option series  smpl
option switch  graph    0
option string  header
option string  subheader
option string  footer
*
local real     rss0 minrss test1 test2
local series   cf0 cf2
local vector   sorted
local index    baselist
local integer  threshnobs pass
local vector   gammas rssthr1 rssthr2
local integer  startl endl
local report   preport
local string   ltitle
local series   gammaser ftestser
*
declare real %%rss1 %%ftest1_0 %%gamma1_1
declare real %%rss2 %%ftest2_1 %%gamma1_2 %%gamma2_2
*
* Estimate the base model by FE
*
if .not.%defined(breakvar).or..not.%defined(threshvar) {
   disp "###Usage:@PanelThresh(THRESHVAR=threshold variable,BREAKVAR=breaking regressor,...)"
   return
}
*
* To reproduce the results in the paper, use the first PREG, not the
* second.
*
*preg(method=fixed,print=print,smpl=%period(t)<%panelobs()) depvar start end
preg(method=fixed,smpl=smpl,print=print) depvar start end
compute startl=%regstart(),endl=%regend()
*
* Save the sum of squared residuals and the regressor list
*
compute rss0=%rss
compute baselist=%rlfromtable(%eqntable(0))
*
* Sort unique values for the threshold
*
@uniquevalues(values=sorted,smpl=smpl) threshvar startl endl
compute threshnobs=%size(sorted)
*
* Single threshold
*
* Create vector of gammas
*
@MakeGammas(obs=threshnobs,quantiles=qn,trim=t1,sorted=sorted) gammas
dim rssthr1(qn)
compute minrss=rss0
do i=1,qn
   set cf0 = breakvar*(threshvar<=gammas(i))
   preg(method=fixed,smpl=smpl,noprint) depvar startl endl
   # baselist cf0
   compute rssthr1(i)=%rss
end do i
compute %%rss1    =%minvalue(rssthr1)
compute %%gamma1_1=gammas(%minindex(rssthr1))
compute %%ftest1_0=(rss0-%%rss1)/(%%rss1/%nobs)

if %defined(title)
   compute ltitle=title
else
   compute ltitle="Panel Threshold Analysis, Threshold Variable "+%l(threshvar)

if graph {
   set gammaser 1 qn = gammas(t)
   set ftestser 1 qn = (rss0-rssthr1(t))/(rssthr1(t)/%nobs)
   scatter(style=line,header=header,subheader=subheader,footer=footer)
   # gammaser ftestser 1 qn
}
*
* Double threshold, by sequential analysis. First, take the one break
* value as given and look for optimal second break. Then, take the
* second break as given and replace the first break.
*
compute %%gamma1_2=%%gamma1_1
dim rssthr2(qn)
do pass=1,2
   do i=1,qn
      compute test1=%min(gammas(i),%%gamma1_2)
      compute test2=%max(gammas(i),%%gamma1_2)
      set cf0 = breakvar*(threshvar<=test1)
      set cf2 = breakvar*(threshvar>test2)
      preg(method=fixed,smpl=smpl,noprint) depvar startl endl
      # baselist cf0 cf2
      compute rssthr2(i)=%rss
   end do i
   compute %%gamma2_2=gammas(%minindex(rssthr2))
   if pass==1
      compute %%gamma1_2=%%gamma2_2
   else
      break
end
compute %%rss2    =%minvalue(rssthr2)
compute %%ftest2_1=(%%rss1-%%rss2)/(%%rss2/%nobs)

if .not.print
   return

report(use=preport,action=define,title=ltitle)
report(use=preport,atrow=1,atcol=1,span) ltitle
report(use=preport,atrow=2,atcol=1) "Single Threshold" %%gamma1_1
report(use=preport,atrow=3,atcol=1) "Sum of Squared Residuals" %%rss1
report(use=preport,atrow=4,atcol=1) "F Test vs 0 Threshold" %%ftest1_0
report(use=preport,atrow=5,atcol=1) "Double Threshold" %%gamma1_2 %%gamma2_2
report(use=preport,atrow=6,atcol=1) "Sum of Squared Residuals" %%rss2
report(use=preport,atrow=7,atcol=1) "F Test vs 1 Threshold" %%ftest2_1
report(use=preport,atcol=2,tocol=3,action=format,width=12)
report(use=preport,action=show)
end
