*! adfmaxur v1.0 JOtero 09dec2016 *! v1.1 CFBaum 24mar2017 Mata logic added *! v1.2 CFBaum 25mar2017 reverse logic corrected, onepanel enabled *! v1.3 CFBaum 26jan2018 guard against rounding values in alphas *! v1.3 CFBaum 02apr2018 fix to guard against very low test statistics for which p-value approximation may not be good *! v1.3 CFBaum 24sep2022 remove mata clear capture program drop adfmaxur program adfmaxur, rclass version 13 syntax varname(ts) [if] [in] [, TREND MAXLag(integer -1) noPRINT] loc qq qui 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)'a/adfmaxur.mtx" tempvar y yr trd trdr 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 p1t1 = (`maxlag'/`t') local p2t2 = (`maxlag'/`t')^2 local p3t3 = (`maxlag'/`t')^3 local p4t4 = (`maxlag'/`t')^4 local firstobs = `maxlag'+1 local case = cond("`trend'" == "", 1, 2) // cfb add if qualifier quietly gen double `y' = `varlist' if `touse' // Generate the reverse realisation of the time series before markout for lags qui gen double `yr' = . mata: st_view(y=.,.,"`y'","`touse'"); st_view(yr=.,.,"`yr'","`touse'"); yr[.,.]=y[rows(y)..1,.] qui gen byte `sspl' = 1 if `touse' markout `touse' L(1/`firstobs').`sspl' tempvar lag pval adff adfr adfmax 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 `adff' = . quietly gen `adfr' = . quietly gen `adfmax' = . 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 adfmax regression with no lags // Store adf forward, reverse and max t-statistics // Case 1: Model includes constant if `case'==1 { `qq' regress D.`y' L1.`y' if `touse' 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 `adff' = _b[L1.`y']/_se[L1.`y'] in 1 `qq' regress D.`yr' L1.`yr' if `touse' quietly replace `adfr' = _b[L1.`yr']/_se[L1.`yr'] in 1 quietly replace `adfmax' = max(`adff',`adfr') in 1 } // Case 2: Model includes constant and trend else if `case'==2 { `qq' regress D.`y' L1.`y' `trd' if `touse' 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 `adff' = _b[L1.`y']/_se[L1.`y'] in 1 `qq' regress D.`yr' L1.`yr' `trd' if `touse' quietly replace `adfr' = _b[L1.`yr']/_se[L1.`yr'] in 1 quietly replace `adfmax' = max(`adff',`adfr') in 1 } // Run the adfmax regression augmented with lags (model includes constant) if `maxlag'>0 { forvalues i=1/`maxlag' { local ii = `i' + 1 // Case 1: Model includes constant if `case'==1 { `qq' regress D.`y' L1.`y' L(1/`i')D.`y' if `touse' 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 `adff' = _b[L1.`y']/_se[L1.`y'] in `ii' `qq' regress D.`yr' L1.`yr' L(1/`i')D.`yr' if `touse' quietly replace `adfr' = _b[L1.`yr']/_se[L1.`yr'] in `ii' quietly replace `adfmax' = max(`adff',`adfr') in `ii' } // Case 2: Model includes constant and trend else if `case'==2 { `qq' regress D.`y' L1.`y' L(1/`i')D.`y' `trd' if `touse' 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 `adff' = _b[L1.`y']/_se[L1.`y'] in `ii' `qq' regress D.`yr' L1.`yr' L(1/`i')D.`yr' `trd' if `touse' quietly replace `adfr' = _b[L1.`yr']/_se[L1.`yr'] in `ii' quietly replace `adfmax' = max(`adff',`adfr') 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 adfmaxm1 = `adfmax'[`store_fix'+1] local adfmaxm2 = `adfmax'[`store_aic'+1] local adfmaxm3 = `adfmax'[`store_sic'+1] local adfmaxm4 = `adfmax'[`store_gts05'+1] local adfmaxm5 = `adfmax'[`store_gts10'+1] // Compute the 1, 5 and 10% critical values along with the associated p-value of the adfmax 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 "adfmaxur.mtx" loc qhl 1 `t1' `t2' `t3' `t4' `p1t1' `p2t2' `p3t3' `p4t4' loc adfmaxv `adfmaxm1' `adfmaxm2' `adfmaxm3' `adfmaxm4' `adfmaxm5' mata: adfmaxur1("`qhl'", "`adfmaxv'") 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] = `adfmax'[`store_`w''+1] } if "`print'" != "noprint" { display as result _n "Leybourne (1995) 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 "Model includes constant" loc treat "constant" } else if `case'==2 { display "Model includes constant and trend" loc treat "constant and trend" } display as text "{hline 74}" display "Criteria{col 12}Lags{col 20}ADFmax 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 adfmax_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 = `lp' - `fp' + 1 return local varname = "`varlist'" restore end // DISABLE mata: mata clear version 13 mata void adfmaxur1(string scalar consts, string scalar adfmaxv) { fref = st_global("fileref") fk = fopen(fref, "r") adfmaxur = fgetmatrix(fk) // vl = fget(fk) // vars = tokens(vl)' fclose(fk) alphas = adfmaxur[.,1] fix_inormal = invnormal(alphas) kon = strtoreal(tokens(consts))' adfmax = strtoreal(tokens(adfmaxv)) 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 = adfmaxur[.,(fc..lc)] fix_qhat = c1_fix_cb * kon fix_qhat2 = fix_qhat:^2 fix_dist = abs(fix_qhat :- adfmax[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 } // fix to guard against very low test statistics for which p-value approximation may not be good pv_fix = 0 if (adfmax[c1] > fix_qhat[1]) { 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 = adfmax[c1] , adfmax[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