*! radf.   v1.10 CFBaum 30jan2022
*! radf     v1.9 CFBaum 27oct2021
*! radf     v1.8 CFBaum 24aug2021
*! radf     v1.7 CFBaum 18may2021
*! radf     v1.6 CFBaum 14oct2020
*! radf     v1.5 CFBaum 23jul2020
*! radf     v1.4 CFBaum 10jul2020
*! radf     v1.3 CFBaum 07jul2020
*! radf     v1.2 CFBaum 27jun2020
*! radf     v1.1 CFBaum 01jun2020
*! radf     v1.0 JOtero 04apr2020

capture program drop radf
mata: mata clear
program radf, rclass
version 13

syntax varname(ts) [if] [in] , [ PREfix(string) MAXLag(integer -1) ///
                                 CRITerion(string)   ///
								 WINdow(integer -1) ///
								 SEED(integer -1) ///
								 BS BOOT(integer -1) ///
								 GRAPH PRINT]

marksample touse
_ts  tvar panelvar `if' `in', sort onepanel
markout `touse' `varlist'
quietly tsreport if `touse'
if r(N_gaps) {
   display in red "sample may not contain gaps"
   exit
}
loc tsfmt `r(tsfmt)'

global fradf "`c(sysdir_plus)'r/radf.mtx"
global frsadf "`c(sysdir_plus)'r/rsadf.mtx"
global frgsadf "`c(sysdir_plus)'r/rgsadf.mtx"
global frball_sadf "`c(sysdir_plus)'r/rall_bsadf.mtx"

// validate variables created by the routine
if "`prefix'" != "" {
		confirm name `prefix'
}
loc vvv rolADF SADF BSADF
loc vvvv SADF BSADF
loc ttt 90 95 99
loc newvars
loc newvars2
loc cvv
local cmdline : copy local 0
loc bootrepl 0

tempvar en y trd yboot resid
tempname tstat res allres b

// Generate a time trend, which starts in 1, regardless of the start of the sample period
quietly gen `trd' = sum(`touse')
qui count if `touse'
loc ntouse `r(N)'
local enobs  = `trd'[_N]
// prefix, graph not available if enobs>600
loc nok = `enobs'<=600
if !`nok' {
	di _n "BSADF critical values and graphs not available for `enobs' observations"
	loc prefix
	loc graph
}
   if "`prefix'" != "" {
		foreach v of local vvv {
			confirm new var `prefix'`v'
		}
// also must validate BSADF 95, 90, indicators for exceedance at 90 and 95
		confirm new var `prefix'BSADF_95
		confirm new var `prefix'BSADF_90
		confirm new var `prefix'Exceeding95
		confirm new var `prefix'Exceeding90
    }
   loc i 0
   foreach v of local vvv {
		tempvar `v'
		qui gen ``v'' = .
		loc newvars "`newvars' ``v''"
		loc i = `i' + 1
		if (`i'>1) {
			if "`prefix'" == "" {
				loc newvars2 "`newvars2' ``v''"
// disable graph if no prefix given
 				loc graph
			}
		else {
				qui generate  `prefix'`v' = .
				loc newvars2 "`newvars2' `prefix'`v'"
			}
		}
   }
   loc i 0
   foreach t of local ttt {
		tempvar cv`t'
		qui gen `cv`t'' = .
		loc cvv "`cvv' `cv`t''"
		label var `cv`t'' " right-tail CV[`t'%]"
	}
if `maxlag'==-1 {
// use more conservative estimate from Schwert (1989, 2002)
   local maxlag = int(4*(`enobs'/100)^0.25)
//   di as res _n "maxlag = `maxlag'"
}
// remove trend option
local case 1
// local case = cond("`trend'" == "", 1, 2)

// BS logic
if  ("`bs'" != "" | `boot'>0) {
   local bootrepl = cond(`boot'==-1, 199, `boot')
}

/*
//   di as res _n "Bootstrap replications = `bootrepl'"
}
else {
	local bootrepl = `boot'
}
*/
if `seed'!=-1 {
   local seednum = `seed'
   set seed `seednum'
// di as res _n "seed = `seed'"
}

* Window size for rolling estimation; see Table 3 in Caspi (2017)
	loc wwidauto = floor(`enobs'*(0.01 + 1.8/sqrt(`enobs')))
if `window'==-1 {
//	loc wauto 1
	loc wwid = `wwidauto'
}
else {
//	loc wauto 0
	loc wwid `window'
}
// ensure sufficient data for lags so that dof>0; reset if window set automatically
if `wwid' <=`maxlag'+ (`case' + 1) {
	di as err _n "Error: `maxlag' too large for window  = `wwid'"
	loc maxlag = `wwid' - `case'
	di as err "Resetting maxlag to `maxlag'"
}

if "`criterion'" == "" | "`criterion'" == "fix" {
   local ncrit = 1
}
else if "`criterion'" == "aic" {
   local ncrit = 2
}
else if "`criterion'" == "sic" {
   local ncrit = 3
}
else if "`criterion'" == "gts05" {
   local ncrit = 4
}
else if "`criterion'" == "gts10" {
   local ncrit = 5
}
else {
	di "Error in specifying criterion: `criterion'"
	error 198
}
//  need to be able to define lags; do not apply touse
qui gen double `y' = `varlist'
g `en' = _n
// do not reference the tsset var, but rather the trend
// must also deal with sum(touse) when it returns to 0
su `en' if `trd'>0 & !mi(`trd') & `touse', mean
// ensure sufficient observations to handle max lag order + differencing at begnning of sample
loc first = max(`r(min)',`maxlag'+2)
loc LAST = `r(max)'
// di as err "*** `lastobs' `LAST'"
// l `tvar' `en' `touse' `varlist', noobs sep(0)
// add one to pick up last possible window
loc last = `LAST' - `wwid' + 1
// ensure sufficient observations
if `first' > `last' {
	di as err _n "Insufficient observations for rolling analysis"
	error 198
}
loc crits FIX AIC SIC GTS05 GTS10
loc rets  fix aic sic gts05 gts10
loc crit : word `ncrit' of `crits'
loc retl : word `ncrit' of `rets'

// l `en' `tvar' `y' in `first'/`LAST', sep(0) noobs
//  di as err _n " first last  l-f+1   LAST   L-f+1  maxlag   criterion   wwid"
//  di as err "  `first' |  `last'  | `=`last'-`first'+1'  |  `LAST'  |  `=`LAST'-`first'+1' |    `maxlag'   |     `criteria'    |  `wwid'"

tempvar  dy ly
qui g `dy' = D.`y' if `touse'
qui g `ly' = L.`y' if `touse'
// di as err "dy, ly"
// l `dy' `ly' if `touse'
loc lagv
forv i=1/`maxlag' {
		tempvar ly`i'
		qui g `ly`i'' = DL`i'.`y'
		loc lagv "`lagv' `ly`i''"
}
// pass in entire series, ignoring touse
mata: lagadf(`first',`last',`wwid',`LAST',`case',`ncrit',`maxlag',"`touse'","`ly'","`dy'","`lagv'")
mata: psy(`enobs',`first',`last',`wwid',`LAST',`case',"`touse'","`ly'","`dy'","`lagv'","`newvars'")

loc tee = `LAST' - `first' + 1

di as res _n "Right-tail ADF statistics for `varlist' with first observations " `tsfmt' `tvar'[`first'] " - " `tsfmt' `tvar'[`last']
di as res _n "Number of obs = `tee'  lag selection[`crit']  maxlag = `maxlag'  window = `wwid' periods"

mat `tstat' = __adfstat \ __pwystat \ __psystat
mat `res' = __adfcv[1,2..4] \ __sadfcv[1,2..4] \  __gsadfcv[1,2..4]

if ("`bs'" == "" & `boot'==-1) {
mat `allres' = `tstat', `res'
mat rownames `allres' = ADF0 SADF GSADF
mat colnames `allres' = Test Tab90 Tab95 Tab99
matlist `allres', format(%8.4f) twidth(5)
di _n as res "Test: ADF0, SADF (PWY,2011), GSADF (PSY,2015)"
di    as res "Tab : right-tail tabulated critical values for 90, 95, 99 confidence levels"
di    as res "      from Vasilopoulos, Pavlidis, Spavound and Martínez-García (2020)"

/*
di as res _n "ADF " _col(18) "= "  %8.4f __adfstat " `adf_s' Right-tail CVs:  90% "               %8.4f __adfcv[1,2]   " 95% "  %8.3f __adfcv[1,3]   "  99% "   %8.4f __adfcv[1,4]

di as res _n "SADF (PWY,2011) " _col(18)  "= " %8.4f __pwystat " `sadf_s' Right-tail CVs:  90% "  %8.4f __sadfcv[1,2]  " 95% "  %8.4f __sadfcv[1,3]  "  99% "   %8.4f __sadfcv[1,4]

di as res _n "GSADF (PSY,2015) " _col(18) "= " %8.4f __psystat " `gsadf_s' Right-tail CVs:  90% " %8.4f __gsadfcv[1,2] " 95% "  %8.4f __gsadfcv[1,3] "  99% "   %8.4f __gsadfcv[1,4]
*/
if (`wwid' != `wwidauto') {
	loc winwarn "Note: critical values based on auto window width = `wwidauto', not specified width = `wwid'"
	di as res _n "`winwarn'"
	}
}
// BS logic
if ("`bs'" != "" | `boot'>0) {
	qui gen double `yboot' = .
	if `maxlag' == 0 {
	qui regress D.`y' if `touse'
	}
	else {
	qui regress D.`y' L(1/`maxlag')D.`y' if `touse'
	}
	matrix `b' = e(b)
	qui predict double `resid' if `touse', residuals
	local maxlag2 = `maxlag'+2
// call original radf routine (renamed to radf0)
	loc cmd "radf0 `yboot' in `maxlag2'/`enobs', maxlag(`maxlag') window(`wwid')"
	mata: bradf("`y'","`yboot'","`dy'","`b'",`maxlag',`enobs',`bootrepl',"`resid'","`touse'", "`lagv'","`cmd'")
	mat `allres' = `tstat', __bsres, `res'
	mat rownames `allres' = ADF0 SADF GSADF
	mat colnames `allres' = Test RTMC90 RTMC95 RTMC99 Tab90 Tab95 Tab99
	matlist `allres', format(%8.4f) twidth(5)
	di _n as res "Test: ADF0, SADF (PWY,2011), GSADF (PSY,2015)"
	di    as res "RTMC: right-tail Monte Carlo critical values for 90, 95, 99 percentiles"
	di    as res "      based on wild bootstrap with `bootrepl' replications"
	di    as res "Tab : right-tail tabulated critical values for 90, 95, 99 confidence levels"
	di    as res "      from Vasilopoulos, Pavlidis, Spavound and Martínez-García (2020)"
	if (`wwid' != `wwidauto') {
		loc winwarn "Note: tabulated critical values based on auto window width = `wwidauto', not specified width = `wwid'"
		di as res _n "`winwarn'"
	}
}
// end BS logic

// base on supadf not missing
loc v1: word 2 of `newvars'
su `tvar' if !mi(`v1'), mean
loc fc0 = `r(min)'
loc lc0 = `r(max)'
loc fd = `r(min)' - (`wwid' - 1)
loc ld = `r(max)' - (`wwid' - 1)
loc rstripe
loc cstripe
forv j=`fd'/`ld' {
	loc lbl = string(`j',"`tsfmt'")
	loc rstripe "`rstripe' `lbl'"
}
forv j=`fc0'/`lc0' {
	loc lbl = string(`j',"`tsfmt'")
	loc cstripe "`cstripe' `lbl'"
}
su `en' if !mi(`v1'), mean
loc fr = `r(min)'  - (`wwid' - 1)
loc lr = `r(max)'  - (`wwid' - 1)
loc fc = `r(min)'
loc lc = `r(max)'

// matlist __results

tempname res2 opt2 nobs2
mat `res2' = __results[`fr'..`lr', `fc'..`lc']
mat rownames `res2' = `rstripe'
mat colnames `res2' = `cstripe'
mata: st_matrix("__optlag",optlag)
mat `opt2' = __optlag[`fr'..`lr', `fc'..`lc']
mat rownames `opt2' = `rstripe'
mat colnames `opt2' = `cstripe'
mat `nobs2' =  __nobs[`fr'..`lr', `fc'..`lc']
mat rownames `nobs2' = `rstripe'
mat colnames `nobs2' = `cstripe'

if (`nok') {
if "`prefix'" != "" {
	mata: pwy(`tee',"`cvv'")
	loc frow = __to1
	loc lrow = __to2
//	if ("`print'" == "print") {
//		di _n "gensupADF CVs for 10, 5, 1"
//		l `cvv' in `frow'/`lrow', noobs sep(0)
//	}
	foreach v of local vvvv {
		qui replace `prefix'`v' = ``v''
	}
// add 5%, 10%^ CV for gensupadf to created variables
	qui g `prefix'BSADF_95 = `cv95' if !mi(`cv95')
	qui g `prefix'BSADF_90 = `cv90' if !mi(`cv90')
	qui g `prefix'Exceeding95 = (`prefix'BSADF > `cv95') if !mi(`cv95')
	qui g `prefix'Exceeding90 = (`prefix'BSADF > `cv90') if !mi(`cv90')
	if ("`print'" == "print") {
		di as res _n "Labeled by endpoint of estimation window"
		l `tvar' `newvars2' `prefix'BSADF_95 `prefix'Exceeding95 in `fc'/`lc', sep(0) noobs abb(20)
	}
}
if ("`print'" == "print") {
	matlist `res2', tit("Estimation window from rowdate to coldate, crit[`crit']")
	if (`ncrit'>1) {
		matlist `opt2', tit("Optimal lags, crit[`crit']") nohalf
	}
	else {
		di as res _n  "Computed with crit[FIX], `maxlag' lags"
	}
	matlist `nobs2', tit("Number of observations in test")
}
}
if (`nok') {
if "`graph'" != "" {
	lab var `tvar' " "
	loc i 0
	di _n
	foreach v of local newvars2 {
		loc i=`i'+1
		loc vn: word `i' of `newvars2'
		lab var `v' "`vn'"
		loc adfv
		loc note
		if (`i' == 1) {
	// CVs for sadf, bsadf
		tempvar adf90 adf95 adf99
		loc j 1
		foreach l in 90 95 { // 1  {
			loc j = `j' + 1
			g `adf`l'' = __adfcv[1,`j']
			label var `adf`l'' " right-tail CV[`l'%]"
			loc adfv "`adfv' `adf`l''"
		}
		loc vti "Date-stamping explosive behavior of `varlist', SADF test"
		}
		else {
			loc adfv `cv90' `cv95'
			if (`wwid' != `wwidauto') {
				loc note "note("`winwarn'",size(vsmall))"
			}
		loc vti "Date-stamping explosive behavior of `varlist', BSADF test"
		}
		tsline `v' `adfv' if !mi(`v'), ylab(,angle(0) labs(small)) scheme(s2mono) graphregion(color(white)) ti("`vti'", size(medium)) ///
		 xlab(#6, labs(small)) legend(row(1) size(small) symxsize(medium) region(lcolor(white))) saving(`v', replace) ///
		 name(`v',replace) `note'
	}
}
}
return local cmd = "radf"
return local cmdlne = "`cmdline'"
return local varname = "`varlist'"
return local first "`first'"
return local last "`last'"
return scalar nobs = `enobs'
return scalar N = `tee'
return scalar maxlag = `maxlag'
return scalar window = `wwid'
return local crit = "`criterion'"
return scalar adfstat = `allres'[1,1]
return scalar sadfstat = `allres'[2,1]
return scalar gsadfstat = `allres'[3,1]
return matrix radfstats = `allres'
return scalar ntests = __tests
if (`bootrepl' > 0) {
	return scalar nreps = `bootrepl'
}
if (`seed' != -1) {
	return scalar seed = `seed'
}

end
// -----------------------------------------------------------------------------
version 13
mata:
mata set matastrict on
void  lagadf(real scalar  first,    ///
             real scalar  last,     ///
		     real scalar  wwid,     ///
		     real scalar  LAST,     ///
		     real scalar  kase,     ///
			 real scalar  krit,     ///
		     real scalar  maxlag,   ///
		     string scalar touse,   ///
		     string scalar lyy,     ///
		     string scalar dyy,     ///
		     string scalar lagvv)

{
	external real matrix optlag

	st_view(ly=., ., tokens(lyy)) // touse)
	st_view(dy=., ., tokens(dyy)) // touse)
	st_view(lagd=., ., tokens(lagvv)) // touse)

// krit=1, just calc ADF for maxlag
	if (krit == 1) {
		optlag = J(LAST, LAST, maxlag)
		return
	}
// -----------------------------------------------------------------------------
// krit > 1
		iota = J(LAST,1,1)
		trd = J(LAST,1,.)
// kase 2: include trd *DISABLE*
/*
		if (kase==2) {
			for(t=1;t<=LAST;t++){
				trd[t]=t
			}
			iota = iota, trd
		}
*/
		// reverse results matrix to contain subscripts for start, end of estimation window
		optlag = J(LAST, LAST, .)
// -----------------------------------------------------------------------------
// krit=2 (aic), estimate each and keep track of min value
	if (krit==2) {

	for(t=first;t<=last;t++) {
		wo = t + wwid - 1
		for(tt=wo;tt<=LAST;tt++) {
			minval = 1e10
			minlag = -1
// loop over lags for each combination of start/end
// calc standard DF, 0 lags
			X = (ly[t..tt],iota[t..tt,.])
			wye = dy[t..tt]
			xxinv = invsym(X'*X)
			beta = xxinv*(X'*wye)
			s = sqrt(((wye - X*beta)'*(wye - X*beta))/(rows(wye)-cols(X)))
			dful = beta[1]/(s * sqrt(xxinv[1,1]))
			if (dful<.) {
			rss = (wye - X*beta)'*(wye - X*beta)
			aic = log(rss/rows(wye)) + 2*(cols(X)-1)/rows(wye)
				if (aic < minval) {
					minval = aic
					minlag = 0
				}
			}
		for(i=1;i<=maxlag;i++) {
			X = (ly[t..tt],lagd[t..tt,1..i],iota[t..tt,.])
			xxinv = invsym(X'*X)
			beta = xxinv*(X'*wye)
			s = sqrt(((wye - X*beta)'*(wye - X*beta))/(rows(wye)-cols(X)))
			dful = beta[1]/(s * sqrt(xxinv[1,1]))
			if (dful<.) {
// aic = log(e(rss)/e(N)) + (2*(e(df_m)))/e(N)
// rss: resid sum sq
// n: rows(wye)
// df_m: cols(X)-1
			rss = (wye - X*beta)'*(wye - X*beta)
			aic = log(rss/rows(wye)) + 2*(cols(X)-1)/rows(wye)
				if (aic < minval) {
					minval = aic
					minlag = i
				}
			}
// printf("t, tt %4.0f %4.0f lag = %4.0g dful= %10.3g aic = %10.5g minval, minlag = %10.5g %5.0f\n", t,tt, i, dful, aic, minval, minlag)
//		end loop over lags
			}
//      store optimal lag in results matrix
			optlag[t,tt] = minlag
//	end loop over start/end
			}
			}
// end for krit=2 (aic)
		}
// -----------------------------------------------------------------------------
// krit=3 (sic), estimate each and keep track of min value
	if (krit==3) {

	for(t=first;t<=last;t++) {
		wo = t + wwid - 1
		for(tt=wo;tt<=LAST;tt++) {
			minval = 1e10
			minlag = -1
// loop over lags for each combination of start/end
// add standard DF, 0 lags, unless crit==1
			X = (ly[t..tt],iota[t..tt,.])
			wye = dy[t..tt]
			xxinv = invsym(X'*X)
			beta = xxinv*(X'*wye)
			s = sqrt(((wye - X*beta)'*(wye - X*beta))/(rows(wye)-cols(X)))
			dful = beta[1]/(s * sqrt(xxinv[1,1]))
			if (dful<.) {
			rss = (wye - X*beta)'*(wye - X*beta)
			sic = log(rss/rows(wye)) + log(rows(wye))*(cols(X)-1)/rows(wye)
			if (sic < minval) {
				minval = sic
				minlag = 0
				}
			}
		for(i=1;i<=maxlag;i++) {
			X = (ly[t..tt],lagd[t..tt,1..i],iota[t..tt,.])
			xxinv = invsym(X'*X)
			beta = xxinv*(X'*wye)
			s = sqrt(((wye - X*beta)'*(wye - X*beta))/(rows(wye)-cols(X)))
			dful = beta[1]/(s * sqrt(xxinv[1,1]))
			if (dful<.) {
// sic = log(e(rss)/e(N)) + (log(e(N))*(e(df_m)))/e(N)
// rss: resid sum sq
// n: rows(wye)
// df_m: cols(X)-1
			rss = (wye - X*beta)'*(wye - X*beta)
			sic = log(rss/rows(wye)) + log(rows(wye))*(cols(X)-1)/rows(wye)
			if (sic < minval) {
				minval = sic
				minlag = i
		//		printf("lag = %4.0g sic = %10.5g minval, minlag = %10.5g %5.0f\n", i, sic, minval, minlag)
			}
			}
//		end loop over lags
			}
//      store optimal lag in results matrix
			optlag[t,tt] = minlag
//	end loop over start/end
			}
			}
// end for krit=3 (sic)
		}
// -----------------------------------------------------------------------------
// krit=4 (gts05), krit=5 (gts10): estimate each and keep track of max lag surpassing threshold

	if (krit==4 | krit==5) {
	crit05 = 0.05
	crit10 = 0.10
// "maxlag, wwid, first, last, LAST"
// maxlag, wwid, first, last, LAST
	for(t=first;t<=last;t++) {
		wo = t + wwid - 1
		for(tt=wo;tt<=LAST;tt++) {

			maxlag4 = 0
			maxlag5 = 0
// loop over lags for each combination of start/end
// add standard DF, 0 lags, unless crit==1
			X = (ly[t..tt],iota[t..tt,.])
			wye = dy[t..tt]
			xxinv = invsym(X'*X)
			beta = xxinv*(X'*wye)
			s = sqrt(((wye - X*beta)'*(wye - X*beta))/(rows(wye)-cols(X)))
			dful = beta[1]/(s * sqrt(xxinv[1,1]))
			if (dful<.) {
				pval = 2 * ttail((rows(wye)-cols(X)),abs(dful))
				if (pval < crit05) maxlag4 = 0
				if (pval < crit10) maxlag5 = 0
				}
		for(i=1;i<=maxlag;i++) {
			X = (ly[t..tt],lagd[t..tt,1..i],iota[t..tt,.])
			xxinv = invsym(X'*X)
			beta = xxinv*(X'*wye)
			s = sqrt(((wye - X*beta)'*(wye - X*beta))/(rows(wye)-cols(X)))
			dful = beta[1]/(s * sqrt(xxinv[1,1]))
			if (dful<.) {
				pval = 2 * ttail((rows(wye)-cols(X)),abs(dful))
				if (pval < crit05) maxlag4 = i
				if (pval < crit10) maxlag5 = i
				}
// 			end loop over lags
			}
// t,tt, maxlag4, maxlag5
//      store optimal lag in results matrix
			if (krit==4) optlag[t,tt] = maxlag4
			if (krit==5) optlag[t,tt] = maxlag5
//		end loop over tt
			}
//	end loop over start/end
		}
// end for krit=4,5
	}

// end function
}

// -----------------------------------------------------------------------------
void  psy(real scalar  enobs,    ///
          real scalar  first,    ///
          real scalar  last,     ///
		  real scalar  wwid,     ///
		  real scalar  LAST,     ///
		  real scalar  kase,     ///
		  string scalar touse,   ///
		  string scalar lyy,     ///
		  string scalar dyy,     ///
		  string scalar lagvv,   ///
		  string scalar newvars)

{
	external real matrix optlag
	external real scalar adfstat, pwystat, psystat
	real rowvector adf, sadf, gsadf
	string scalar etoile
	real matrix sumstats

	st_view(ly=., ., tokens(lyy))
	st_view(dy=., ., tokens(dyy))
	st_view(lagd=., ., tokens(lagvv))
	st_view(rolsupgen=., ., tokens(newvars))

	iota = J(LAST,1,1)
//produce trd over full range
	trd = J(rows(ly),1,.)
	for(t=1;t<=rows(ly);t++){
		trd[t]=t
	}
	rolladf = J(LAST,1,.)
	// case 2: include trd *DISABLE*
//	if(kase==2) iota = iota, trd
	ntim=0
//	sumstats matrix contains the three stats and their CVS
	sumstats = J(3,4,.)
//  results matrix contains subscripts for start, end of estimation window
	results = J(LAST, LAST, .)
//  nobs matrix records number of observations in each estimation
	nobs = J(LAST, LAST, .)
// mark first defined row of results
	fr=0
	for(t=first;t<=last;t++) {
		wo = t + wwid - 1
		for(tt=wo;tt<=LAST;tt++) {
			lord = optlag[t,tt]
// dfuller `yy' in `t'/`tt', lags(`r(lag_`retl')') `trend'
// 		qui regress D.`y'  L1.`y' L(1/`i')D.`y' `tt' if `touse'
// adjust for zero lags as optimal
			if (lord>0){
				X = (ly[t..tt],lagd[t..tt,1..lord],iota[t..tt,.])
			}
			else {
				X = (ly[t..tt],iota[t..tt,.])
			}
			wye = dy[t..tt]
			xxinv = invsym(X'*X)
			beta = xxinv*(X'*wye)
			s = sqrt(((wye - X*beta)'*(wye - X*beta))/(rows(wye)-cols(X)))
			dful = beta[1]/(s * sqrt(xxinv[1,1]))
// grab first test results for each t, label as element t for roladf
			if(tt==wo) rolladf[t] = dful
			results[t,tt] = dful
			nobs[t,tt] = rows(wye)
			if (dful<.) {
				if (fr==0) {
					fr=t
				}
				ntim++
			}
 /*
			t, tt, results[t,tt]
			"--"
			rows(wye)
			wye,X
			beta
			eps
			"dof"
			dof=(rows(wye)-cols(X)-1)
			dof
			s
			dful

	if(t==first & tt==LAST) {
		lord, t, tt, dful
		wye, X
	}
*/
		}
	}
// access adf CVs
	fref = st_global("fradf")
	fk = fopen(fref, "r")
	adfcv = fgetmatrix(fk)
	fclose(fk)
// access sadf CVs
	fref = st_global("frsadf")
	fk = fopen(fref, "r")
	sadfcv = fgetmatrix(fk)
	fclose(fk)
// access gsadf CVs
	fref = st_global("frgsadf")
	fk = fopen(fref, "r")
	gsadfcv = fgetmatrix(fk)
	fclose(fk)
// base on full sample entering touse
	adf = cvret(adfcv, enobs)
	sadf = cvret(sadfcv, enobs)
	gsadf = cvret(gsadfcv, enobs)
	st_matrix("__adfcv",adf)
	st_matrix("__sadfcv",sadf)
	st_matrix("__gsadfcv",gsadf)
	st_numscalar("__tests",ntim)
	st_matrix("__results", results)
	st_matrix("__nobs", nobs)
	adfstat = results[fr,LAST]
	st_numscalar("__adfstat",adfstat)
/*
	etoile="   "
	if (adfstat > adf[2]) etoile="*  "
	if (adfstat > adf[3]) etoile="** "
	if (adfstat > adf[4]) etoile="***"
	st_local("adf_s",etoile)
*/
	sumstats[1,1] = adfstat
	sumstats[1,2..4] = adf[2..4]
//  maximum of rolladf
	rollstat = colmax(rolladf)
	st_numscalar("__rollstat",rollstat)
// identify range of non-missing values in roladf
	e=colminmax(mm_which((rolladf:<.):*trd[1..rows(rolladf)]))
	fr1 = e[1,1]
	fr2 = e[2,1]
// identify range of non-missing values in supadf
	f=colminmax(mm_which((results[fr,.]':<.):*trd[1..LAST]))
	to1 = f[1,1]
	to2 = f[2,1]
//	frow = results[fr,.]'
// load supadf with transposed first row of results
	rolsupgen[1..to2,2] = results[fr,.]'
// adjust first column to align with 2d, 3d
	rolsupgen[to1..to2,1] = rolladf[fr1..fr2,1]
	pwystat = rowmax(results[fr,.])
	st_numscalar("__pwystat",pwystat)
/*
	etoile="   "
	if (pwystat > sadf[2]) etoile="*  "
	if (pwystat > sadf[3]) etoile="** "
	if (pwystat > sadf[4]) etoile="***"
	st_local("sadf_s",etoile)
*/
	sumstats[2,1] = pwystat
	sumstats[2,2..4] = sadf[2..4]
// reverse subscripts
	psyadf = colmax(results)'
// load genadf with psyadf vector
	rolsupgen[1..to2,3] = psyadf
	psystat = colmax(psyadf)
	st_numscalar("__psystat",psystat)
/*
	etoile="   "
	if (psystat > gsadf[2]) etoile="*  "
	if (psystat > gsadf[3]) etoile="** "
	if (psystat > gsadf[4]) etoile="***"
	st_local("gsadf_s",etoile)
*/
	sumstats[3,1] = psystat
	sumstats[3,2..4] = gsadf[2..4]
	st_matrix("radfstats",sumstats)
// pass limits for use in pwy
	st_numscalar("__to1",to1)
	st_numscalar("__to2",to2)
}
// -----------------------------------------------------------------------------
real rowvector cvret(real matrix cvmat, ///
					 real scalar enobs)
{
	real scalar lrow, frac
	real rowvector cvl, cvh, cvnt
	if (enobs <= 600) {
		lrow = enobs - 5
		return(cvmat[lrow,.])
	}
// row 595: 600
// row 609: 2000
	if (enobs >= 2000) {
		return(cvmat[609,.])
	}
// T >600, <2000
	nh = floor(enobs/100) + 589
	frac = mod(enobs,100)/100
	cvl = cvmat[nh,.]
	cvh = cvmat[nh+1,.]
	cvnt = cvl + frac :* (cvh - cvl)
	return(cvnt)
}
// -----------------------------------------------------------------------------
void  pwy(real   scalar  tee,   ///
          string scalar  cvv)
{
	st_view(cvvars=., ., tokens(cvv))

	// access gensupADF CVs from bsadf
	fref = st_global("frball_sadf")
	fk = fopen(fref, "r")
	bscv = fgetmatrix(fk)
	fclose(fk)
	// locate rows with tee in col 5
	chunk = select(bscv, bscv[.,5]:==tee)[.,2..4]
	nr = rows(chunk)
	to1 = st_numscalar("__to1")
	to2 = st_numscalar("__to2")
	need = to2 - to1 + 1
	fr = to1 + (need-nr)
	cvvars[fr..to2,.] = chunk
}
// -----------------------------------------------------------------------------
void bradf(string scalar yy,
		   string scalar yb,
		   string scalar dyy,
           string scalar bb,
		   real scalar maxlag,
		   real scalar enobs,
		   real scalar bootrepl,
		   string scalar resid,
		   string scalar touse,
		   string scalar ldyy,
		   string scalar cmd)
{
		external real scalar adfstat, pwystat, psystat
		real colvector index

		st_view(y=.,.,tokens(yy))
		st_view(yboot=.,.,tokens(yb))
		st_view(ldy=.,.,tokens(ldyy))
		st_view(dy=.,.,tokens(dyy))
		st_view(eps=.,.,tokens(resid),touse)
		b = st_matrix(bb)
		stat = J(bootrepl, 3, .)
/*
		maxlag
		ldy
		b = st_matrix(bb)
		b
		cmd
*/

// lose maxlag+2 to define L(maxlag)D.y
		maxlag2 = maxlag + 2
		nindex = enobs - maxlag2 + 1
		index = J(nindex,1,.)
		e = J(nindex,1,.)
		for(i=1;i<=bootrepl;i++) {
			index = runiformint(nindex, 1, maxlag2, enobs)
			w = rnormal(nindex,1,0,1)
			e = eps[index] :* w
//		maxlag2, enobs
			if (maxlag == 0) {
				xb = e
			}
			else {
				xb = ldy[maxlag2..enobs,.] * b[1..maxlag]' + e
			}
			yboot[.,.] = y
/*
		"yboot"
		yboot[maxlag2..enobs]
		"xb"
		xb
*/
			for(j=maxlag2;j<=enobs;j++) yboot[j] = yboot[j-1] + xb[j-maxlag-1]
//		"yboot"
//		yboot[maxlag2..enobs]
			stata(cmd, 1)
			stat[i,.] = adfstat, pwystat, psystat
		}
//		stat
// mm_quantile from moremata
		p90 = mm_quantile(stat, 1, 0.90)
		p95 = mm_quantile(stat, 1, 0.95)
		p99 = mm_quantile(stat, 1, 0.99)
//		"90, 95, 99 percentiles of bootstrap distribution"
		bsres = (p90 \ p95 \ p99)'
//		bsres
		st_matrix("__bsres",bsres)
		st_local("tee",strofreal(nindex))
}
// -----------------------------------------------------------------------------

end