*! xtnonlincoint 1.0.0 21jun2026 *! Nonlinear panel cointegration tests with structural breaks and CSD *! Author: Merwan Roudane (merwanroudane920@gmail.com) *! github.com/merwanroudane *! *! Subcommands *! ecm : Nonlinear error-correction-based panel cointegration test *! Omay, Emirmahmutoglu & Denaux (2017, Economics Letters 157, 1-4) *! *! fffff : Fractional Frequency Flexible Fourier Form panel cointegration test *! Olayeni, Tiwari & Wohar (2021, Applied Economics Letters 28, 482-486) *! *! all : run both tests on the same model program define xtnonlincoint, rclass version 14.0 gettoken sub 0 : 0, parse(" ,") if ("`sub'"=="") { di as error "specify a subcommand: {bf:ecm}, {bf:fffff} or {bf:all}" di as error "see {bf:help xtnonlincoint}" exit 198 } if ("`sub'"=="ecm") { xtnonlincoint_ecm `0' } else if ("`sub'"=="fffff") { xtnonlincoint_fffff `0' } else if ("`sub'"=="all") { xtnonlincoint_all `0' } else { di as error "unknown subcommand {bf:`sub'}" di as error "valid subcommands are {bf:ecm}, {bf:fffff}, {bf:all}" exit 198 } return add end *----------------------------------------------------------------------------- * ECM : nonlinear error-correction based panel cointegration test (Omay 2017) *----------------------------------------------------------------------------- program define xtnonlincoint_ecm, rclass version 14.0 syntax varlist(min=2 numeric) [if] [in] [, /// Lags(integer 1) /// VARlags(integer 1) /// Breps(integer 299) /// Seed(integer 12345) /// Level(cilevel) /// TRend /// GRaph /// noPRINT ] gettoken dv xv : varlist qui xtset local ivar "`r(panelvar)'" local tvar "`r(timevar)'" if ("`ivar'"=="" | "`tvar'"=="") { di as error "data are not {bf:xtset}; run {bf:xtset panelvar timevar} first" exit 459 } marksample touse markout `touse' `dv' `xv' `ivar' `tvar' local trd = ("`trend'"!="") if (`lags'<0) local lags 0 if (`varlags'<1) local varlags 1 if (`breps'<19) local breps 19 mata: _xtnlc_ecm_driver("`dv'","`xv'","`ivar'","`tvar'","`touse'", /// `lags', `varlags', `breps', `seed', `trd') if (scalar(__bal)==0) { di as error "{bf:xtnonlincoint ecm} requires a balanced panel" di as error "(every panel must have the same number of time periods)" exit 459 } matrix colnames __ind = panel tau pvalue cv10 cv5 cv1 if ("`print'"=="") { _xtnlc_ecm_table, dv(`dv') xv(`xv') lags(`lags') varlags(`varlags') /// breps(`breps') level(`level') trend(`trd') } if ("`graph'"!="") { _xtnlc_ecm_graph, level(`level') } * ---- returns ---- return scalar N = scalar(__N) return scalar T = scalar(__T) return scalar lags = `lags' return scalar varlags = `varlags' return scalar breps = `breps' return scalar stat = scalar(__gstat) return scalar p = scalar(__pval) return scalar cv10 = scalar(__gcv10) return scalar cv5 = scalar(__gcv5) return scalar cv1 = scalar(__gcv1) return local depvar "`dv'" return local indepvars "`xv'" return local test "Nonlinear ECM-based panel cointegration (Omay et al. 2017)" return local cmd "xtnonlincoint ecm" return matrix indstat = __ind return matrix bootdist = __bdist capture scalar drop __bal __N __T __gstat __pval __gcv10 __gcv5 __gcv1 capture matrix drop __ind __bdist end *----------------------------------------------------------------------------- * FFFFF : fractional frequency flexible Fourier form test (Olayeni 2021) *----------------------------------------------------------------------------- program define xtnonlincoint_fffff, rclass version 14.0 syntax varlist(min=2 numeric) [if] [in] [, /// MAXLags(integer 1) /// KStep(real 0.1) /// Breps(integer 299) /// BLock(integer 0) /// Seed(integer 12345) /// Level(cilevel) /// TRend /// SPSM /// GRaph /// noPRINT ] gettoken dv xv : varlist qui xtset local ivar "`r(panelvar)'" local tvar "`r(timevar)'" if ("`ivar'"=="" | "`tvar'"=="") { di as error "data are not {bf:xtset}; run {bf:xtset panelvar timevar} first" exit 459 } marksample touse markout `touse' `dv' `xv' `ivar' `tvar' local trd = ("`trend'"!="") if (`maxlags'<0) local maxlags 0 if (`kstep'<=0) local kstep 0.1 if (`breps'<19) local breps 19 mata: _xtnlc_fffff_driver("`dv'","`xv'","`ivar'","`tvar'","`touse'", /// `maxlags', `kstep', `breps', `block', `seed', `trd') if (scalar(__bal)==0) { di as error "{bf:xtnonlincoint fffff} requires a balanced panel" di as error "(every panel must have the same number of time periods)" exit 459 } matrix colnames __ind = panel kss khat phat pvalue cv10 cv5 cv1 matrix colnames __spsm = step gmstat pvalue minkss khat panel if ("`print'"=="") { _xtnlc_fffff_table, dv(`dv') xv(`xv') maxlags(`maxlags') /// breps(`breps') level(`level') trend(`trd') spsm(`=("`spsm'"!="")') } if ("`graph'"!="") { _xtnlc_fffff_graph, level(`level') } * ---- returns ---- return scalar N = scalar(__N) return scalar T = scalar(__T) return scalar maxlags = `maxlags' return scalar breps = `breps' return scalar stat = scalar(__gstat) return scalar p = scalar(__pval) return scalar cv10 = scalar(__gcv10) return scalar cv5 = scalar(__gcv5) return scalar cv1 = scalar(__gcv1) return scalar nstat = scalar(__nstat) return local depvar "`dv'" return local indepvars "`xv'" return local test "FFFFF panel cointegration (Olayeni et al. 2021)" return local cmd "xtnonlincoint fffff" return matrix indstat = __ind return matrix spsm = __spsm return matrix bootdist = __bdist capture scalar drop __bal __N __T __gstat __pval __gcv10 __gcv5 __gcv1 __nstat capture matrix drop __ind __spsm __bdist end *----------------------------------------------------------------------------- * ALL : run both tests *----------------------------------------------------------------------------- program define xtnonlincoint_all, rclass version 14.0 syntax varlist(min=2 numeric) [if] [in] [, /// Lags(integer 1) /// VARlags(integer 1) /// MAXLags(integer 1) /// KStep(real 0.1) /// Breps(integer 299) /// BLock(integer 0) /// Seed(integer 12345) /// Level(cilevel) /// TRend /// SPSM /// GRaph ] di as text _n "{hline 78}" di as text "{bf:[1/2]} Nonlinear ECM-based test" di as text "{hline 78}" xtnonlincoint_ecm `varlist' `if' `in', lags(`lags') varlags(`varlags') /// breps(`breps') seed(`seed') level(`level') `trend' `graph' local e_stat = r(stat) local e_p = r(p) di as text _n "{hline 78}" di as text "{bf:[2/2]} FFFFF test" di as text "{hline 78}" xtnonlincoint_fffff `varlist' `if' `in', maxlags(`maxlags') kstep(`kstep') /// breps(`breps') block(`block') seed(`seed') level(`level') /// `trend' `spsm' `graph' return scalar ecm_stat = `e_stat' return scalar ecm_p = `e_p' return scalar fffff_stat = r(stat) return scalar fffff_p = r(p) return local cmd "xtnonlincoint all" end *----------------------------------------------------------------------------- * ECM results table *----------------------------------------------------------------------------- program define _xtnlc_ecm_table syntax , dv(string) xv(string) lags(int) varlags(int) breps(int) /// level(int) trend(int) local nr = rowsof(__ind) di as text _n "{hline 74}" di as text "Nonlinear ECM-based panel cointegration test" di as text "Omay, Emirmahmutoglu & Denaux (2017)" di as text "{hline 74}" di as text "Dependent variable : " as result "`dv'" di as text "Regressors : " as result "`xv'" local dt "demeaned" if (`trend') local dt "detrended" di as text "Deterministics : " as result "`dt'" /// as text " ECM lags: " as result "`lags'" /// as text " VAR lags: " as result "`varlags'" di as text "Bootstrap reps : " as result "`breps'" /// as text " (sieve, cross-section dependent)" di as text "Panels (N) : " as result %5.0f scalar(__N) /// as text " Periods (T): " as result %5.0f scalar(__T) di as text "{hline 74}" di as text "H0: no cointegration in any panel" di as text "H1: nonlinear (LSTAR) error correction in some panel" di as text "{hline 74}" di as text %12s "Panel" " " %10s "tau(c)" " " %8s "p-value" /// " " %8s "cv10%" " " %8s "cv5%" " " %8s "cv1%" di as text "{hline 74}" forvalues i = 1/`nr' { local pid = __ind[`i',1] local ta = __ind[`i',2] local pv = __ind[`i',3] local c10 = __ind[`i',4] local c5 = __ind[`i',5] local c1 = __ind[`i',6] local st "" if (`pv'<0.10) local st "*" if (`pv'<0.05) local st "**" if (`pv'<0.01) local st "***" di as text %12.0g `pid' " " as result %10.4f `ta' "`st'" /// _col(28) as result %8.3f `pv' " " as text %8.3f `c10' /// " " %8.3f `c5' " " %8.3f `c1' } di as text "{hline 74}" local gs = scalar(__gstat) local gp = scalar(__pval) local gst "" if (`gp'<0.10) local gst "*" if (`gp'<0.05) local gst "**" if (`gp'<0.01) local gst "***" di as result %12s "Group mean" " " as result %10.4f `gs' "`gst'" /// _col(28) as result %8.3f `gp' " " as text %8.3f scalar(__gcv10) /// " " %8.3f scalar(__gcv5) " " %8.3f scalar(__gcv1) di as text "{hline 74}" di as text "Stars: * p<.10 ** p<.05 *** p<.01 (bootstrap, upper tail)." di as text "Reject H0 (=> cointegration) when tau(c) exceeds the critical value." end *----------------------------------------------------------------------------- * FFFFF results table *----------------------------------------------------------------------------- program define _xtnlc_fffff_table syntax , dv(string) xv(string) maxlags(int) breps(int) level(int) /// trend(int) spsm(int) local nr = rowsof(__ind) di as text _n "{hline 78}" di as text "Fractional Frequency Flexible Fourier Form (FFFFF) cointegration test" di as text "Olayeni, Tiwari & Wohar (2021)" di as text "{hline 78}" di as text "Dependent variable : " as result "`dv'" di as text "Regressors : " as result "`xv'" local dt "level (constant)" if (`trend') local dt "trend" di as text "Deterministics : " as result "`dt'" /// as text " Max lags: " as result "`maxlags'" di as text "Bootstrap reps : " as result "`breps'" /// as text " (stationary, cross-section dependent)" di as text "Panels (N) : " as result %5.0f scalar(__N) /// as text " Periods (T): " as result %5.0f scalar(__T) di as text "{hline 78}" di as text "H0: no cointegration in any panel (residual has a unit root)" di as text "H1: nonlinear (KSS/ESTAR) cointegration with smooth Fourier breaks" di as text "{hline 78}" di as text %12s "Panel" " " %9s "KSS" " " %6s "k" " " %5s "lag" /// " " %8s "p-value" " " %8s "cv10%" " " %8s "cv5%" di as text "{hline 78}" forvalues i = 1/`nr' { local pid = __ind[`i',1] local ks = __ind[`i',2] local kh = __ind[`i',3] local ph = __ind[`i',4] local pv = __ind[`i',5] local c10 = __ind[`i',6] local c5 = __ind[`i',7] local st "" if (`pv'<0.10) local st "*" if (`pv'<0.05) local st "**" if (`pv'<0.01) local st "***" di as text %12.0g `pid' " " as result %9.3f `ks' "`st'" /// _col(26) as result %6.1f `kh' " " %5.0f `ph' /// " " as result %8.3f `pv' " " as text %8.3f `c10' " " %8.3f `c5' } di as text "{hline 78}" local gs = scalar(__gstat) local gp = scalar(__pval) local gst "" if (`gp'<0.10) local gst "*" if (`gp'<0.05) local gst "**" if (`gp'<0.01) local gst "***" di as result %12s "Group mean" " " as result %9.3f `gs' "`gst'" /// _col(40) as result %8.3f `gp' " " as text %8.3f scalar(__gcv10) /// " " %8.3f scalar(__gcv5) di as text "{hline 78}" di as text "Stars: * p<.10 ** p<.05 *** p<.01 (bootstrap, lower tail)." di as text "Reject H0 (=> cointegration) when KSS is below the critical value." if (`spsm') { local ns = rowsof(__spsm) di as text _n "{hline 64}" di as text "Sequential Panel Selection Method (SPSM)" di as text "{hline 64}" di as text %6s "Step" " " %10s "group KSS" " " %8s "p-value" /// " " %9s "min KSS" " " %6s "k" " " %10s "panel" di as text "{hline 64}" forvalues i = 1/`ns' { local sp = __spsm[`i',1] local gm = __spsm[`i',2] local pv = __spsm[`i',3] local mk = __spsm[`i',4] local kh = __spsm[`i',5] local pid = __spsm[`i',6] local st "" if (`pv'<0.10) local st "*" if (`pv'<0.05) local st "**" if (`pv'<0.01) local st "***" di as result %6.0f `sp' " " %10.3f `gm' " " %8.3f `pv' "`st'" /// _col(40) as result %9.3f `mk' " " %6.1f `kh' " " %10.0g `pid' } di as text "{hline 64}" di as text "Series removed (top rows, while significant) drive the" di as text "panel cointegration; later non-significant series do not." } end *----------------------------------------------------------------------------- * ECM graph *----------------------------------------------------------------------------- program define _xtnlc_ecm_graph syntax , level(int) preserve clear qui svmat double __ind, names(col) qui gen long _ix = _n qui svmat double __bdist, name(bdist) local gs = scalar(__gstat) local gc = scalar(__gcv5) capture { twoway (bar tau _ix, barwidth(0.6) color(navy%70)) /// (connected cv5 _ix, lpattern(dash) lcolor(cranberry) /// msymbol(none)), /// legend(order(1 "tau(c)" 2 "5% critical value") rows(1) size(small)) /// ytitle("MWALD statistic") xtitle("Panel index") /// title("Individual cointegration statistics", size(medium)) /// name(xtnlc_ecm_ind, replace) nodraw } capture { twoway (histogram bdist1, frequency color(navy%50)) , /// xline(`gs', lcolor(cranberry) lpattern(solid) lwidth(medthick)) /// xline(`gc', lcolor(forest_green) lpattern(dash)) /// xtitle("Bootstrap group-mean tau(c)") ytitle("Frequency") /// title("Bootstrap null distribution", size(medium)) /// note("Solid red = observed; dashed green = 5% critical value") /// name(xtnlc_ecm_boot, replace) nodraw } capture graph combine xtnlc_ecm_ind xtnlc_ecm_boot, /// title("Nonlinear ECM-based panel cointegration test", size(medium)) /// name(xtnlc_ecm, replace) if (_rc) { di as text "(graph could not be produced in this environment)" } restore end *----------------------------------------------------------------------------- * FFFFF graph *----------------------------------------------------------------------------- program define _xtnlc_fffff_graph syntax , level(int) preserve clear qui svmat double __ind, names(col) qui gen long _ix = _n qui svmat double __spsm, name(sp) local gs = scalar(__gstat) capture { twoway (bar kss _ix, barwidth(0.6) color(maroon%70)) /// (connected cv5 _ix, lpattern(dash) lcolor(navy) msymbol(none)), /// legend(order(1 "KSS" 2 "5% critical value") rows(1) size(small)) /// ytitle("KSS statistic") xtitle("Panel index") /// title("Individual cointegration statistics", size(medium)) /// name(xtnlc_fffff_ind, replace) nodraw } capture { twoway (connected sp2 sp1, lcolor(navy) mcolor(navy)) , /// ytitle("Group-mean KSS") xtitle("SPSM step") /// title("Sequential panel selection", size(medium)) /// note("Series peeled off one per step (most stationary first)") /// name(xtnlc_fffff_spsm, replace) nodraw } capture graph combine xtnlc_fffff_ind xtnlc_fffff_spsm, /// title("FFFFF panel cointegration test", size(medium)) /// name(xtnlc_fffff, replace) if (_rc) { di as text "(graph could not be produced in this environment)" } restore end *============================================================================= * Mata implementation *============================================================================= version 14.0 mata: struct _xtnlc_pv { real matrix B real matrix W real matrix Z real colvector dnu real scalar khat real scalar phat } // ---- quantile (type-7 linear interpolation) ---- real scalar _xtnlc_quantile(real colvector x, real scalar pr) { real colvector s real scalar n, h, fl n = rows(x) if (n==0) { return(.) } if (n==1) { return(x[1]) } s = sort(x,1) h = pr*(n-1)+1 fl = floor(h) if (fl<1) { fl = 1 } if (fl>=n) { return(s[n]) } return(s[fl] + (h-fl)*(s[fl+1]-s[fl])) } // ---- lag row builder: [M[s-1,], M[s-2,], ..., M[s-q,]] ---- real rowvector _xtnlc_lagrow(real matrix M, real scalar s, real scalar q) { real rowvector out real scalar j out = M[s-1,] for (j=2; j<=q; j=j+1) { out = (out, M[s-j,]) } return(out) } // ---- panel boundaries: returns N x 2 [start,end], requires sorted id ---- real matrix _xtnlc_bounds(real colvector ID) { real colvector ids real matrix bnd real scalar N, i, a, n, cnt ids = uniqrows(ID) N = rows(ids) n = rows(ID) bnd = J(N, 2, 0) a = 1 for (i=1; i<=N; i=i+1) { bnd[i,1] = a cnt = 0 while (a<=n) { if (ID[a]!=ids[i]) { break } cnt = cnt + 1 a = a + 1 } bnd[i,2] = bnd[i,1] + cnt - 1 } return(bnd) } // ======================= ECM (Omay et al. 2017) ======================= // MWALD tau(c) for one panel; returns . if not computable real scalar _xtnlc_ecm_panel(real colvector y, real matrix X, real scalar p, /// real scalar trend) { real scalar T, k, t0, m, i, tt, c, j, base, ncol real colvector u, dep, b, e, bb real matrix D, R, XtXi, V real scalar s2, r1, r2, v11, v22, v21, t1, t2, denom T = rows(y) k = cols(X) if (trend) { D = (J(T,1,1), range(1,T,1), X) } else { D = (J(T,1,1), X) } bb = invsym(quadcross(D,D))*quadcross(D,y) u = y - D*bb t0 = p + 2 if (t0>T) { return(.) } m = T - t0 + 1 ncol = 3 + k + p*(1+k) if (m<=ncol) { return(.) } R = J(m, ncol, 0) dep = J(m, 1, 0) for (i=1; i<=m; i=i+1) { tt = t0 + i - 1 dep[i] = y[tt] - y[tt-1] R[i,1] = 1 R[i,2] = u[tt-1] R[i,3] = u[tt-1]*u[tt-1] for (c=1; c<=k; c=c+1) { R[i,3+c] = X[tt,c] - X[tt-1,c] } base = 3 + k for (j=1; j<=p; j=j+1) { R[i,base+1] = y[tt-j] - y[tt-j-1] for (c=1; c<=k; c=c+1) { R[i,base+1+c] = X[tt-j,c] - X[tt-j-1,c] } base = base + 1 + k } } XtXi = invsym(quadcross(R,R)) b = XtXi*quadcross(R,dep) e = dep - R*b s2 = quadcross(e,e)/(m-ncol) V = s2*XtXi r1 = b[2] r2 = b[3] v11 = V[2,2] v22 = V[3,3] v21 = V[3,2] if (v22<=0 | v11<=0) { return(.) } denom = v11 - v21*v21/v22 if (denom<=0) { return(.) } t1 = (r1 - r2*v21/v22) t1 = t1*t1/denom t2 = 0 if (r2<0) { t2 = r2*r2/v22 } return(t1+t2) } // estimate VAR(q) on dz; fill Bmat (q*kk x kk) and centered residuals W (L x kk) void _xtnlc_var(real matrix dz, real scalar q, real matrix Bmat, real matrix W) { real scalar R, kk, s, L real matrix Xl, Yl R = rows(dz) kk = cols(dz) L = R - q if (L<1) { Bmat = J(q*kk, kk, 0) W = J(1, kk, 0) return } Xl = J(L, q*kk, 0) Yl = J(L, kk, 0) for (s=q+1; s<=R; s=s+1) { Xl[s-q,] = _xtnlc_lagrow(dz, s, q) Yl[s-q,] = dz[s,] } Bmat = invsym(quadcross(Xl,Xl))*quadcross(Xl,Yl) W = Yl - Xl*Bmat W = W :- mean(W) } // reconstruct z* (T x kk) from VAR coefficients and bootstrap innovations real matrix _xtnlc_recon(real matrix Bmat, real matrix Wstar, real scalar q, /// real scalar T) { real scalar R, kk, s real matrix dzs, zs real rowvector lr R = T - 1 kk = cols(Wstar) dzs = J(R, kk, 0) for (s=q+1; s<=R; s=s+1) { lr = _xtnlc_lagrow(dzs, s, q) dzs[s,] = lr*Bmat + Wstar[s-q,] } zs = J(T, kk, 0) for (s=2; s<=T; s=s+1) { zs[s,] = zs[s-1,] + dzs[s-1,] } return(zs) } void _xtnlc_ecm_driver(string scalar dv, string scalar xv, string scalar idv, /// string scalar tvv, string scalar tousev, real scalar p, real scalar q, /// real scalar B, real scalar seed, real scalar trend) { real colvector Y, ID, TT, ord, ids, tau, bdist, bdv, ip, c10, c5, c1 real matrix X, bnd, BTAU, OUT real scalar N, i, b, T0, s, e, L, gstat real colvector yp, ys, idx, col, colf, rowv, rowf real matrix Xp, Xs, zp, dz, zs, Wst struct _xtnlc_pv colvector PV real matrix Bm, Wm Y = st_data(., dv, tousev) X = st_data(., tokens(xv), tousev) ID = st_data(., idv, tousev) TT = st_data(., tvv, tousev) ord = order((ID,TT), (1,2)) Y = Y[ord] X = X[ord,] ID = ID[ord] TT = TT[ord] bnd = _xtnlc_bounds(ID) ids = uniqrows(ID) N = rows(bnd) // balanced check T0 = bnd[1,2] - bnd[1,1] + 1 for (i=1; i<=N; i=i+1) { if ((bnd[i,2]-bnd[i,1]+1)!=T0) { st_numscalar("__bal", 0) return } } st_numscalar("__bal", 1) // observed statistics + per-panel VAR for bootstrap tau = J(N,1,.) PV = _xtnlc_pv(N) for (i=1; i<=N; i=i+1) { s = bnd[i,1] e = bnd[i,2] yp = Y[|s \ e|] Xp = X[|s,. \ e,.|] tau[i] = _xtnlc_ecm_panel(yp, Xp, p, trend) zp = (yp, Xp) dz = zp[|2,. \ rows(zp),.|] - zp[|1,. \ rows(zp)-1,.|] Bm = J(0,0,0) Wm = J(0,0,0) _xtnlc_var(dz, q, Bm, Wm) PV[i].B = Bm PV[i].W = Wm } bdv = select(tau, tau:<.) gstat = mean(bdv) // bootstrap L = rows(PV[1].W) BTAU = J(B, N, .) rseed(seed) for (b=1; b<=B; b=b+1) { idx = ceil(runiform(L,1):*L) idx = idx + (idx:<1) idx = idx - (idx:>L):*(idx:-L) for (i=1; i<=N; i=i+1) { Wst = (PV[i].W)[idx,] zs = _xtnlc_recon(PV[i].B, Wst, q, T0) ys = zs[,1] Xs = zs[|1,2 \ rows(zs),cols(zs)|] BTAU[b,i] = _xtnlc_ecm_panel(ys, Xs, p, trend) } } // bootstrap group distribution bdist = J(B,1,.) for (b=1; b<=B; b=b+1) { rowv = BTAU[b,]' rowf = select(rowv, rowv:<.) if (rows(rowf)>0) { bdist[b] = mean(rowf) } } bdv = select(bdist, bdist:<.) // individual p-values and critical values (upper tail) ip = J(N,1,.) c10 = J(N,1,.) c5 = J(N,1,.) c1 = J(N,1,.) for (i=1; i<=N; i=i+1) { col = BTAU[,i] colf = select(col, col:<.) if (rows(colf)>0 & tau[i]<.) { ip[i] = mean(colf:>=tau[i]) c10[i] = _xtnlc_quantile(colf, 0.90) c5[i] = _xtnlc_quantile(colf, 0.95) c1[i] = _xtnlc_quantile(colf, 0.99) } } OUT = (ids, tau, ip, c10, c5, c1) st_matrix("__ind", OUT) st_matrix("__bdist", bdv) st_numscalar("__N", N) st_numscalar("__T", T0) st_numscalar("__gstat", gstat) st_numscalar("__pval", mean(bdv:>=gstat)) st_numscalar("__gcv10", _xtnlc_quantile(bdv,0.90)) st_numscalar("__gcv5", _xtnlc_quantile(bdv,0.95)) st_numscalar("__gcv1", _xtnlc_quantile(bdv,0.99)) } // ======================= FFFFF (Olayeni et al. 2021) ======================= // stage-1 cointegrating-regression residual real colvector _xtnlc_resid(real colvector y, real matrix Z, real scalar trend) { real scalar T real matrix D real colvector b T = rows(y) if (trend) { D = (J(T,1,1), range(1,T,1), Z) } else { D = (J(T,1,1), Z) } b = invsym(quadcross(D,D))*quadcross(D,y) return(y - D*b) } // nonlinear ADF + Fourier fit; returns (t_gamma, ssr, m, ncol) real rowvector _xtnlc_kss_fit(real colvector nu, real scalar k, real scalar p, /// real scalar trend) { real scalar T, t0, m, i, tt, j, base, ncol, s2, gamma, segamma real colvector dep, b, e real matrix R, XtXi, V T = rows(nu) t0 = p + 2 if (t0>T) { return((.,.,.,.)) } m = T - t0 + 1 ncol = 4 + p if (trend) { ncol = ncol + 1 } if (m<=ncol) { return((.,.,.,.)) } R = J(m, ncol, 0) dep = J(m, 1, 0) for (i=1; i<=m; i=i+1) { tt = t0 + i - 1 dep[i] = nu[tt] - nu[tt-1] R[i,1] = 1 R[i,2] = nu[tt-1]*nu[tt-1]*nu[tt-1] for (j=1; j<=p; j=j+1) { R[i,2+j] = nu[tt-j] - nu[tt-j-1] } base = 2 + p R[i,base+1] = sin(2*pi()*k*tt/T) R[i,base+2] = cos(2*pi()*k*tt/T) if (trend) { R[i,base+3] = tt } } XtXi = invsym(quadcross(R,R)) b = XtXi*quadcross(R,dep) e = dep - R*b s2 = quadcross(e,e)/(m-ncol) V = s2*XtXi if (V[2,2]<=0) { return((.,.,.,.)) } gamma = b[2] segamma = sqrt(V[2,2]) return((gamma/segamma, quadcross(e,e), m, ncol)) } // grid-search over (k,p); returns (t_gamma, khat, phat, aic) real rowvector _xtnlc_kss_search(real colvector nu, real scalar pmax, /// real scalar trend, real colvector kgrid) { real scalar ng, ik, p, aic, bestaic, bestt, bestk, bestp real rowvector fit ng = rows(kgrid) bestaic = . bestt = . bestk = . bestp = . for (ik=1; ik<=ng; ik=ik+1) { for (p=0; p<=pmax; p=p+1) { fit = _xtnlc_kss_fit(nu, kgrid[ik], p, trend) if (fit[2]<. & fit[3]>fit[4] & fit[2]>0) { aic = fit[3]*log(fit[2]/fit[3]) + 2*fit[4] if (bestaic==. | aicL) { cur = 1 } } idx[t] = cur } return(idx) } void _xtnlc_fffff_driver(string scalar dv, string scalar xv, string scalar idv, /// string scalar tvv, string scalar tousev, real scalar pmax, real scalar kstep, /// real scalar B, real scalar blockopt, real scalar seed, real scalar trend) { real colvector Y, ID, TT, ord, ids, KSS, Khat, Phat, kgrid real colvector bm, bmf, bdv, ip, c10, c5, c1, dnustar, nustar, idx real colvector idxset, curK, curKf, active, rr, rrf real matrix X, bnd, BKSS, OUT, SP real scalar N, i, b, T0, s, e, L, gstat, meanbl, t, nstep, mn, mi, jj, ii real scalar om, pv real colvector yp, nu real matrix Zp real rowvector sr struct _xtnlc_pv colvector PV Y = st_data(., dv, tousev) X = st_data(., tokens(xv), tousev) ID = st_data(., idv, tousev) TT = st_data(., tvv, tousev) ord = order((ID,TT), (1,2)) Y = Y[ord] X = X[ord,] ID = ID[ord] TT = TT[ord] bnd = _xtnlc_bounds(ID) ids = uniqrows(ID) N = rows(bnd) T0 = bnd[1,2] - bnd[1,1] + 1 for (i=1; i<=N; i=i+1) { if ((bnd[i,2]-bnd[i,1]+1)!=T0) { st_numscalar("__bal", 0) return } } st_numscalar("__bal", 1) kgrid = range(0.1, 2, kstep) KSS = J(N,1,.) Khat = J(N,1,.) Phat = J(N,1,.) PV = _xtnlc_pv(N) for (i=1; i<=N; i=i+1) { s = bnd[i,1] e = bnd[i,2] yp = Y[|s \ e|] Zp = X[|s,. \ e,.|] nu = _xtnlc_resid(yp, Zp, trend) sr = _xtnlc_kss_search(nu, pmax, trend, kgrid) KSS[i] = sr[1] Khat[i] = sr[2] Phat[i] = sr[3] PV[i].dnu = nu[|2 \ rows(nu)|] - nu[|1 \ rows(nu)-1|] PV[i].Z = Zp PV[i].khat = sr[2] PV[i].phat = sr[3] } curKf = select(KSS, KSS:<.) gstat = mean(curKf) // bootstrap: stationary, common index across panels (CSD) L = T0 - 1 meanbl = blockopt if (meanbl<2) { meanbl = ceil(sqrt(T0)) } if (meanbl<2) { meanbl = 2 } BKSS = J(B, N, .) rseed(seed) for (b=1; b<=B; b=b+1) { idx = _xtnlc_sb_index(L, meanbl) for (i=1; i<=N; i=i+1) { dnustar = (PV[i].dnu)[idx] nustar = J(T0,1,0) for (t=2; t<=T0; t=t+1) { nustar[t] = nustar[t-1] + dnustar[t-1] } // re-impose the first-stage projection on the original Z so the // bootstrap residual carries the same spurious mean reversion as // the observed residual (Olayeni et al. 2021, steps 5-7), then // re-run the SAME (k,p) selection so the frequency search is part // of the bootstrapped statistic (correct size; otherwise the // selection advantage of the observed stat inflates rejections) nustar = _xtnlc_resid(nustar, PV[i].Z, trend) sr = _xtnlc_kss_search(nustar, pmax, trend, kgrid) BKSS[b,i] = sr[1] } } // group p-value (lower tail) + individual bm = J(B,1,.) for (b=1; b<=B; b=b+1) { rr = BKSS[b,]' rrf = select(rr, rr:<.) if (rows(rrf)>0) { bm[b] = mean(rrf) } } bdv = select(bm, bm:<.) ip = J(N,1,.) c10 = J(N,1,.) c5 = J(N,1,.) c1 = J(N,1,.) for (i=1; i<=N; i=i+1) { rr = BKSS[,i] rrf = select(rr, rr:<.) if (rows(rrf)>0 & KSS[i]<.) { ip[i] = mean(rrf:<=KSS[i]) c10[i] = _xtnlc_quantile(rrf, 0.10) c5[i] = _xtnlc_quantile(rrf, 0.05) c1[i] = _xtnlc_quantile(rrf, 0.01) } } // SPSM active = J(N,1,1) SP = J(N, 6, .) nstep = 0 while (sum(active)>0) { idxset = select((1::N), active:==1) curK = KSS[idxset] curKf = select(curK, curK:<.) if (rows(curKf)==0) { break } om = mean(curKf) bm = J(B,1,.) for (b=1; b<=B; b=b+1) { rr = (BKSS[b,])' rr = rr[idxset] rrf = select(rr, rr:<.) if (rows(rrf)>0) { bm[b] = mean(rrf) } } bmf = select(bm, bm:<.) pv = mean(bmf:<=om) mn = . mi = 0 for (jj=1; jj<=rows(idxset); jj=jj+1) { ii = idxset[jj] if (KSS[ii]<.) { if (mi==0) { mn = KSS[ii] mi = ii } else { if (KSS[ii]0) { SP = SP[|1,1 \ nstep,6|] } OUT = (ids, KSS, Khat, Phat, ip, c10, c5, c1) st_matrix("__ind", OUT) st_matrix("__spsm", SP) st_matrix("__bdist", bdv) st_numscalar("__N", N) st_numscalar("__T", T0) st_numscalar("__nstat", rows(select(KSS, KSS:<.))) st_numscalar("__gstat", gstat) st_numscalar("__pval", mean(bdv:<=gstat)) st_numscalar("__gcv10", _xtnlc_quantile(bdv,0.10)) st_numscalar("__gcv5", _xtnlc_quantile(bdv,0.05)) st_numscalar("__gcv1", _xtnlc_quantile(bdv,0.01)) } end