*
* @PerronRodriguez( options )  series start end
* does the Perron-Rodriguez(2003) test for unit roots using GLS
* detrending, allowing for a break point at an unknown date.
*
* Parameters:
*   series      Series to test for unit roots
*   start end   Range of series to use
*
* Options:
*   MAXLAGS=maximum number of lags in the auxiliary regression (depends on obs)
*   MINLAGS=minimum number of lags in the auxiliary regression [0]
*   CRITERION=BIC/[MAIC]/MBIC
*     Selection criterion of lag length
*
*   BREAK=TREND/[BOTH]
*     Determines whether the break is just to the trend rate, or to the
*     trend rate and level.
*   PI=fraction of entries on each end of data to exclude as break points [.10]
*   TEST=[ADF]/MZT/PT
*     Unit root test to use
*   SELECTBY=[TEST]/TSTAT
*     Determines whether to select the break point by the most extreme unit
*     root test [SELECTBY=TEST], or the largest absolute value of the
*     t-statistic for the trend break dummy [SELECTBY=TSTAT]. The authors
*     prefer the default.
*
*   [PRINT]/NOPRINT
*
*  Variables Defined:
*    %CDSTAT=test statistic (selected break point)
*    %%AUTOP=lag length chosen (at selected break point)
*    %%BREAKPOINT=entry chosen
*
*  Reference:
*     Perron and Rodriguez(2003), "GLS Detrending, Efficient Unit Root Tests
*     and Structural Change",  Journal of Econometrics, vol 115, pp 1-27.
*
*  Revision Schedule:
*   09/2009 Written by Tom Doan, Estima
*   05/2021 Revised to use newer @ADFAutoSelect and include report output
*
source adfautoselect.src
source perronngmtests.src
*
procedure PerronRodriguez series start end
type series series
type integer start end
*
option choice  break    2  trend  both
option integer maxlag
option integer minlag   0
option real    pi       .10
option switch  print    1
option choice  selectby 1  test   tstat
option choice  test     1  adf    mzt    pt
option string  title
option choice  crit     2  bic    maic   mbic
fixed vect[label] critlabel(3)
 BIC MAIC MBIC
*
local integer  startl endl pinobs lower upper
local integer  bestlag nobs lags tb glsbreak bestbreak
local series   dy ya
local real     selectstat bestselect beststat
local real     adf mzt pt
local integer  adfautocrit
local string   ltitle mtitle
local report   prreport
*
declare integer %%autop %%breakpoint
*
if .not.%defined(series) {
   disp "Syntax: @PerronRodriguez(options)  Y start end"
   return
}
*
inquire(series=series) startl<<start endl<<end
stats(noprint) series startl endl
if %nmiss>0 {
   disp "@PerronRodriguez can't be used on data with missing values"
   return
}
compute nobs=%nobs
*
if .not.%defined(maxlag)
   compute lags=fix(10*(.01*nobs)^.25)
else
   compute lags=maxlag
*
set dy    startl+1 endl = series-series{1}
*
* lower and upper bounds of regressions, allowing for lags
*
compute lower=startl+lags+1
compute upper=endl
*
* Shortest interval between the ends and the breaks
*
compute pinobs=fix(pi*(upper-lower+1))
*
* Remap option values to match those used by GLSDetrend
*
if break==1
   compute glsbreak=3
else
   compute glsbreak=4

if crit==2
   compute adfautocrit=4
else
if crit==3
   compute adfautocrit=5
else
   compute adfautocrit=2

compute bestselect=%na
do tb=lower+pinobs,upper-pinobs
   @GLSDetrend(tb=tb,break=glsbreak) series startl endl ya
   *
   * We judge most extreme by the most negative, so this is minus
   * the t-statistic on the trend shift.
   *
   compute selectstat=-abs(%tstats(%nreg))
   @ADFAutoSelect(minlag=minlag,maxlag=lags,det=none,crit=adfautocrit,noprint) ya startl endl
   @PerronNgMTests(lags=%%autop,mzt=mzt,adf=adf,pt=pt) ya startl endl
   compute %cdstat=%if(test==1,adf,%if(test==2,mzt,pt))
   *
   * If we are choosing based upon most extreme unit root test,
   * override the minus t-stat with the unit root test.
   *
   if selectby==1
      compute selectstat=%cdstat
   if .not.%valid(bestselect).or.selectstat<bestselect
      compute bestselect=selectstat,beststat=%cdstat,bestbreak=tb,bestlag=%%lags
end do tb
*
if %defined(title)
   compute ltitle=title
else
   compute ltitle="Perron-Rodriguez Test, Series "+%l(series)

if break==1
   compute mtitle="Model I (trend rate change only)"
else
   compute mtitle="Model II (trend and level break)"

report(use=prreport,action=define,title=ltitle)
report(use=prreport,atrow=1,atcol=1,span) ltitle
report(use=prreport,atrow=2,atcol=1,span) "Using data from "+%datelabel(startl)+" to "+%datelabel(endl)
report(use=prreport,atrow=3,atcol=1,span) mtitle
report(use=prreport,atrow=4,atcol=1,span) "Break location chosen by most extreme "+%if(selectby==1,"unit root test","shift dummy t")
report(use=prreport,atrow=5,atcol=1,span) "Using test "+%if(test==1,"ADF",%if(test==2,"MZT","PT"))
report(use=prreport,atrow=6,atcol=1,span) "Choosing lag length using" critlabel(crit)
report(use=prreport,atrow=7,atcol=1,span) "Null is unit root. Reject in left tail"
report(use=prreport,atrow=9,atcol=1) "Break date" %datelabel(bestbreak)
report(use=prreport,atrow=10,atcol=1) "Lags" bestlag
report(use=prreport,atrow=11,atcol=1) "Test" beststat
*
compute %%autop=bestlag
compute %cdstat=beststat
compute %%breakpoint=bestbreak
if print
   report(use=prreport,action=show)
end
