*
* @RegHBreak( options )
*
* Regression post-processor which tests the regression just run for a
* full structural break with fixed regressor bootstrapping of the
* significance level. The test statistics are the Andrews-Quandt
* (maximum) and Andrews-Ploberger (geometric mean) of the F-statistics
* for the break.
*
* Options:
*   REPS=# of residual bootstrap replications [100]
*   PI1=fraction of early entries not examined as break points [.15]
*   PI2=fraction of late entries not examined as break points [PI1]
*   [PRINT]/NOPRINT
*
* Variables Defined:
*   %%BREAKPOINT 	Entry with largest break test
*   %%AQTEST 		Andrews-Quandt test statistic
*   %%APTEST 		Andrews-Ploberger test statistic
*   %%AQSIGNIF  	Bootstrapped significance value for AQ test
*   %%APSIGNIF		Bootstrapped significance value for AP test
*
* References:
*   Bruce Hansen(2000), "Testing for Structural Change in Conditional
*   Models", Journal of Econometrics, vol. 97, no. 1, 93-115, July.
*
*   Andrews and Ploberger(1994), "Optimal Tests When a Nuisance Parameter
*   is Present Only Under the Alternative", Econometrica, vol 62, no 6,
*   1383-1414.
*
* Revision Schedule:
*   08/2010 Written by Tom Doan, Estima
*   05/2015 %%BREAKPOINT changed from REAL to INTEGER
*
procedure RegHBreak depvar start end
type series 	depvar
type integer 	start end
*
option real   	pi1    .15
option real   	pi2    pi1
option integer reps	 100
option switch  print  1
*
local vect[int] reghold
local integer 	piStart piEnd rnobs rep time bestbreak
local series  	yb wstat u
local vect    	x
local symm    	dxx dxeex ss
local vect    	dxe
local real    	rssdiff fstat bestdiff apcount aqcount aqvalue apvalue
local string   ltitle
decl  real     %%aqtest %%aptest %%aqsignif %%apsignif
decl  integer  %%breakpoint
*
compute reghold=%regsave()
*
compute rnobs  =%regend()-%regstart()+1
compute piStart=%regstart()+fix((rnobs-1)*pi1)
compute piEnd  =%regstart()+fix((rnobs-1)*(1-pi2))
set u = %resids
*
dim dxx(%nreg,%nreg) dxeex(%nreg,%nreg) dxe(%nreg)
*
compute bestdiff=0.0
compute apcount=aqcount=0.0
if reps>0
   infobox(action=define,progress,lower=0,upper=reps) $
     "Bootstrapping P-Values"
do rep=0,reps
   infobox(current=rep)
   if rep==0
      set yb = u
   else
      set yb = %ran(1.0)
   *
   * Run the regression on the residuals or bootstrapped residuals
   *
   linreg(lastreg,noprint) yb
   *
   compute dxx  =%zeros(%nreg,%nreg)
   compute dxeex=%zeros(%nreg,%nreg)
   compute dxe  =%zeros(%nreg,1)
   set wstat piStart piEnd = 0.0
   do time=%regstart(),%regend()
      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.or.time>piEnd
         next
      compute ss=(dxx-dxx*%xx*dxx)
      compute rssdiff=%qform(inv(ss),dxe)
      *
      * This follows Hansen's example code. It's actually F/number of
      * restrictions.
      *
      compute fstat=(rssdiff)/((%rss-rssdiff)/(%ndf-%nreg))
      compute wstat(time)=fstat
      if rep==0.and.rssdiff>bestdiff
         compute bestbreak=time,bestdiff=rssdiff
   end do time
   *
   * Compute sup F and overflow-safe exp F statistic
   *
   sstats(max)  piStart piEnd wstat>>aqvalue
   sstats(mean) piStart piEnd exp(.5*(wstat-aqvalue))>>apvalue
   compute apvalue=log(apvalue)+.5*aqvalue
   *
   * Save the test statistics for rep==0; count the number above the
   * saved statistics for the bootstrapped data.
   *
   if rep==0
      compute %%aqtest=aqvalue,%%aptest=apvalue
   else
      compute aqcount=aqcount+(aqvalue>%%aqtest),apcount=apcount+(apvalue>%%aptest)
end do rep
infobox(action=remove)
*
compute %%apsignif=%if(reps>0,apcount/reps,%na)
compute %%aqsignif=%if(reps>0,aqcount/reps,%na)
compute %%breakpoint=bestbreak
*
compute ltitle="AQ-AP Break Tests/Bootstrapped P-Values"
report(use=hreport,action=define,title=ltitle)
report(use=hreport,atrow=1,atcol=1,span) ltitle
report(use=hreport,atrow=2,atcol=1) "Estimated Break" 			%datelabel(bestbreak)
report(use=hreport,atrow=3,atcol=2)                            "Statistic" "Signif"
report(use=hreport,atrow=4,atcol=1) "Andrews-Quandt Test"      %%aqtest    %%aqsignif
report(use=hreport,atrow=5,atcol=1) "Andrews-Ploberger Test"   %%aptest    %%apsignif
report(use=hreport,atcol=2,tocol=2,action=format,picture="#.###")
report(use=hreport,atcol=3,tocol=3,action=format,picture="#.#####")
if print
   report(use=hreport,action=show)

compute %regload(reghold)

end HBreakTest
