*! ersur v1.0 JOtero 09dec2016
*!       v1.1 CFBaum 11mar2017 Mata logic added
*!       v1.2 CFBaum 01apr2017 onepanel added, returns modified
*!       v1.3 CFBaum 26jan2018 guard against rounding values in alphas
*!       v1.3 CFBaum 24sep2022 remove mata clear
capture program drop ersur
program ersur, rclass
version 13

// syntax varname(ts) [if] [in], CASE(integer) [, MAXLag(integer -1) noPRINT]

syntax varname(ts) [if] [in] [, TREND MAXLag(integer -1) noPRINT]

preserve
marksample touse
_ts tvar panelvar `if' `in', sort onepanel
markout `touse' `tvar'

quietly tsreport if `touse'
if r(N_gaps) {
   display in red "sample may not contain gaps"
   exit
}

global fileref "`c(sysdir_plus)'e/ersur.mtx"

tempname cbar rhobar
tempvar  glsy y trd glstrd glscons sspl

// Generate a time trend, which starts in 1, regardless of the start of the sample period
quietly gen `trd' = sum(`touse')

// Determine the number of observations "t" that is used in the response surfaces,
// which is equal to the number of observations of the variable of interest minus one

local lastobs  = `trd'[_N]

local t=`lastobs'-1

if `maxlag'==-1 {
   local maxlag = int(12*(`lastobs'/100)^0.25)
}

// compute regressors of the response surfaces

local t1 = 1/`t'
local t2 = 1/`t'^2
local t3 = 1/`t'^3
local t4 = 1/`t'^4

local p1t = (`maxlag')/`t'
local p2t = (`maxlag'^2)/`t'
local p3t = (`maxlag'^3)/`t'
local p4t = (`maxlag'^4)/`t'

local firstobs = `maxlag'+1

local case = cond("`trend'" == "", 1, 2)

// Case 1: GLS demeaned data
if `case'==1 {
   local cbar -7
   local rhobar = 1+`cbar'/`lastobs'

   quietly gen double `glscons' = 1                           if `trd'==1
   quietly replace    `glscons' = 1 -`rhobar'                 if `trd'>1

   quietly gen double `glsy' = `varlist'                      if `trd'==1
   quietly replace    `glsy' = `varlist'-`rhobar'*L.`varlist' if `trd'>1

   quietly regress    `glsy' `glscons' if `touse' , noconstant
   quietly gen double `y' = `varlist' - _b[`glscons']

   local treat "demeaned"
}

// Case 2: GLS detrended data
else if `case'==2 {
   local cbar -13.5
   local rhobar = 1+`cbar'/`lastobs'

   quietly gen double `glscons' = 1                           if `trd'==1
   quietly replace    `glscons' = 1 -`rhobar'                 if `trd'>1

   quietly gen double `glstrd' = `trd'                        if `trd'==1
   quietly replace    `glstrd' = `trd' -`rhobar'*L.`trd'      if `trd'>1

   quietly gen double `glsy' = `varlist'                      if `trd'==1
   quietly replace    `glsy' = `varlist'-`rhobar'*L.`varlist' if `trd'>1

   quietly regress    `glsy' `glscons' `glstrd' if `touse' , noconstant
   quietly gen double `y' = `varlist' - _b[`glscons'] - _b[`glstrd']*`trd'

   local treat "detrended"
}

qui gen byte `sspl' = 1 if `touse'
markout `touse' L(1/`firstobs').`sspl'
tempvar lag pval ers

tempvar aic    min_aic   lag_aic
tempvar sic    min_sic   lag_sic
tempvar gts05  max_gts05 lag_gts05
tempvar gts10  max_gts10 lag_gts10
tempname res

quietly gen `lag'   = .
quietly gen `aic'   = .
quietly gen `sic'   = .
quietly gen `pval'  = .
quietly gen `ers'   = .

qui tsset
loc tv `r(timevar)'
local rts `r(tsfmt)'
su `tv' if `touse', mean
loc fp `r(min)'
loc lp `r(max)'
loc minp = string(`fp',"`rts'")
loc maxp = string(`lp',"`rts'")

// Determine the optimal number of lags
// Run the ers regression with no lags (no deterministic components necessary)

quietly regress D.`y'  L1.`y' if `touse', noconstant
quietly replace `lag'  = 0                                                 in 1
quietly replace `pval' = 0                                                 in 1
quietly replace `aic'  = log(e(rss)/e(N)) + (2*(e(df_m)))/e(N)             in 1
quietly replace `sic'  = log(e(rss)/e(N)) + (log(e(N))*(e(df_m)))/e(N)     in 1
quietly replace `ers'  = _b[L1.`y']/_se[L1.`y']                            in 1

// Run the ers regression augmented with lags (no deterministic components necessary)

if `maxlag'>0 {

forvalues i=1/`maxlag' {
   quietly regress D.`y'  L1.`y' L(1/`i')D.`y' if `touse', noconstant
   local ii = `i' + 1
   quietly replace `lag'  = `i'                                                       in `ii'
   quietly replace `pval' = (2 * ttail(e(df_r), abs(_b[L`i'D.`y']/_se[L`i'D.`y'])))   in `ii'
   quietly replace `aic'  = log(e(rss)/e(N)) + (2*(e(df_m)))/e(N)                     in `ii'
   quietly replace `sic'  = log(e(rss)/e(N)) + (log(e(N))*(e(df_m)))/e(N)             in `ii'
   quietly replace `ers'  = _b[L1.`y']/_se[L1.`y']                                    in `ii'
   }
}

local store_fix = `maxlag'

// Determine the number of lags based on aic and sic

egen `min_aic' = min(`aic')
egen `min_sic' = min(`sic')

quietly gen `lag_aic' = `lag' if `min_aic' == `aic'
quietly gen `lag_sic' = `lag' if `min_sic' == `sic'

su `lag_aic', mean
local store_aic = r(mean)
su `lag_sic', mean
local store_sic = r(mean)

// Determines the number of lags based on gts05 and gts10

quietly gen     `gts05' = .
quietly gen     `gts10' = .
quietly replace `gts05' = .     if `lag' == .
quietly replace `gts10' = .     if `lag' == .
quietly replace `gts05' = `lag' if `pval'<=0.05
quietly replace `gts10' = `lag' if `pval'<=0.10
quietly replace `gts05' = `lag' if `lag' == 0
quietly replace `gts10' = `lag' if `lag' == 0

egen    `max_gts05' = max(`gts05')
egen    `max_gts10' = max(`gts10')

quietly gen `lag_gts05'  = `max_gts05' if `lag' !=.
quietly gen `lag_gts10'  = `max_gts10' if `lag' !=.

su `lag_gts05', mean
local store_gts05 = r(mean)
su `lag_gts10', mean
local store_gts10 = r(mean)

local ersm1 = `ers'[`store_fix'+1]
local ersm2 = `ers'[`store_aic'+1]
local ersm3 = `ers'[`store_sic'+1]
local ersm4 = `ers'[`store_gts05'+1]
local ersm5 = `ers'[`store_gts10'+1]

// Compute the 1, 5 and 10% critical values along with the associated p-value of the ers statistic.
// For alpha<=0.004 and alpha>=0.996, we use the actual quantile and the 14 observations closest to the
// desired quantile, as there will not be seven observations on either side.

// Response surface coefficients are retrieved from the Stata file "ersur.mtx"

   loc qhl 1 `t1' `t2' `t3' `t4' `p1t' `p2t' `p3t' `p4t'

   loc ersv `ersm1' `ersm2' `ersm3' `ersm4' `ersm5'

   mata: ersur1("`qhl'", "`ersv'")
   mat `res' = ___res
   loc lbl fix aic sic gts05 gts10
   loc i = 0
   foreach w of local lbl  {
		loc i = `i' + 1
		mat `res'[`i',1] = `store_`w''
		mat `res'[`i',2] = `ers'[`store_`w''+1]
	}

if "`print'" != "noprint" {

display as result _n "Elliott, Rothenberg & Stock (1996) test results for `minp' - `maxp'"
display as result "Variable name: `varlist'"
display as result "Ho: Unit root"
display as result "Ha: Stationarity"

if `case'==1 {
   display "GLS demeaned data" // (case 1)"
}
else if `case'==2 {
   display "GLS detrended data" // (case 2)"
}
display as text "{hline 74}"
display "Criteria{col 12}Lags{col 20}ERS stat.{col 33}p-value{col 45}1% cv{col 57}5% cv{col 66} 10% cv"
display as text "{hline 74}"
display as result "FIXED{col 13}" `res'[1,1] "{col 18}" %9.3fc `res'[1,2] "{col 31}" %9.3fc `res'[1,3] "{col 42}" %9.3fc `res'[1,4] "{col 54}" %9.3fc `res'[1,5] "{col 65}" %9.3fc `res'[1,6]
display as result "  AIC{col 13}" `res'[2,1] "{col 18}" %9.3fc `res'[2,2] "{col 31}" %9.3fc `res'[2,3] "{col 42}" %9.3fc `res'[2,4] "{col 54}" %9.3fc `res'[2,5] "{col 65}" %9.3fc `res'[2,6]
display as result "  SIC{col 13}" `res'[3,1] "{col 18}" %9.3fc `res'[3,2] "{col 31}" %9.3fc `res'[3,3] "{col 42}" %9.3fc `res'[3,4] "{col 54}" %9.3fc `res'[3,5] "{col 65}" %9.3fc `res'[3,6]
display as result "GTS05{col 13}" `res'[4,1] "{col 18}" %9.3fc `res'[4,2] "{col 31}" %9.3fc `res'[4,3] "{col 42}" %9.3fc `res'[4,4] "{col 54}" %9.3fc `res'[4,5] "{col 65}" %9.3fc `res'[4,6]
display as result "GTS10{col 13}" `res'[5,1] "{col 18}" %9.3fc `res'[5,2] "{col 31}" %9.3fc `res'[5,3] "{col 42}" %9.3fc `res'[5,4] "{col 54}" %9.3fc `res'[5,5] "{col 65}" %9.3fc `res'[5,6]
di as text "{hline 74}"
}

loc lblu = upper("`lbl'")
mat rownames `res' = `lblu'
mat colnames `res' = Lags ERS_stat p-value critval01 critval05 critval10
return matrix results = `res'
return scalar maxp = `lp'
return scalar minp = `fp'
return local tsfmt "`rts'"
return local treat "`treat'"
return scalar N = `t'
return local varname = "`varlist'"
restore
end

// DISABLE mata: mata clear
version 13
mata
void ersur1(string scalar consts,
            string scalar ersv)

{
fref = st_global("fileref")
fk = fopen(fref, "r")
ersur = fgetmatrix(fk)
// vl = fget(fk)
// vars = tokens(vl)'
fclose(fk)
alphas = ersur[.,1]
fix_inormal = invnormal(alphas)
kon = strtoreal(tokens(consts))'
ers = strtoreal(tokens(ersv))
trd = st_local("trend")
//  cases 1,2: 5 rows for fix, aic, sic, gts05, gts10
colmtx1 = J(10,2,.)
colmtx1 = (2, 10 \ 56, 64 \ 74, 82 \ 20, 28 \ 38, 46 \
          11, 19 \ 65, 73 \ 83, 91 \ 29, 37 \ 47, 55)
off = (trd != "" ? 5 : 0)
// result matrix, including 2 additional columns for labels
    res = J(5,6,.)

for(c1=1; c1<=5; c1++) {
	fc = colmtx1[c1+off,1]
	lc = colmtx1[c1+off,2]
	c1_fix_cb = ersur[.,(fc..lc)]
	fix_qhat = c1_fix_cb * kon
	fix_qhat2 = fix_qhat:^2
	fix_dist = abs(fix_qhat :- ers[c1])
	minindex(fix_dist, 1, fix_place, w)
// guard against rounding
//	if (alphas[fix_place] <= 0.004) {
	if (fix_place <= 7) {
		fix_start = 1
		fix_end   = 15
	}
//	else if (alphas[fix_place] >= 0.996) {
	else if (fix_place >= 215) {
		fix_start = 207
		fix_end   = 221
	}
	else {
		fix_start = fix_place - 7
		fix_end   = fix_place + 7
	}
	y = fix_inormal[fix_start..fix_end]
	iota = J(fix_end - fix_start + 1, 1, 1)
	x = fix_qhat[fix_start..fix_end],fix_qhat2[fix_start..fix_end], iota
	beta = invsym(quadcross(x,x)) * quadcross(x,y)
	kon2 = ers[c1] , ers[c1]^2 , 1
	fix_inormalhat = kon2 * beta
	pv_fix = normal(fix_inormalhat)
	cvfix_01 = fix_qhat[13]
    cvfix_05 = fix_qhat[21]
    cvfix_10 = fix_qhat[31]
	res[c1,3] = pv_fix
	res[c1,4] = cvfix_01
	res[c1,5] = cvfix_05
	res[c1,6] = cvfix_10
	}
	st_matrix("___res", res)
}
end