*! ralsfadl 1.0.1 16may2026 Dr Merwan Roudane *! RALS-Fourier ADL cointegration test * Reference: Yilanci, V., Ulucak, R., Zhang, Y., & Andreoni, V. (2022). * Sustainable Development 31(2): 812-824. * (extends Banerjee, Arcabic & Lee 2017 ADL Fourier model with * the RALS w-augmentation of Im & Schmidt 2008.) * Mirrors the Eviews programs ralsfadl1.prg .. ralsfadl4.prg in the * RALS-FADL.zip supplement and supports 1-4 regressors. *------------------------------------------------------------------------------ program define ralsfadl, rclass version 14.0 capture mata: __rals_loaded() if _rc qui _rals_mata syntax varlist(min=2 ts) [if] [in], [ /// TREND /// MAXLags(integer 3) /// FMAX(integer 5) /// FREQuency(real -1) /// Graph /// noHEADer ] if `fmax' < 1 | `fmax' > 5 { di as error "fmax() must be between 1 and 5" exit 198 } marksample touse qui count if `touse' if r(N) < 25 { di as error "sample too small (N=" r(N) ")" exit 2001 } capture tsset if _rc { di as error "data must be tsset before running ralsfadl" exit 459 } local timevar `r(timevar)' local y : word 1 of `varlist' local x : list varlist - y local kreg : list sizeof x if `kreg' < 1 | `kreg' > 4 { di as error "ralsfadl supports 1-4 regressors (you supplied `kreg')" exit 198 } local model = cond("`trend'"!="", 2, 1) tempname res cv mata: _rals_fadl_engine("`y'", "`x'", "`touse'", `model', `maxlags', `fmax', `frequency', "`res'", "`cv'") local tFADL = `res'[1,1] local tRALS = `res'[1,2] local rho2 = `res'[1,3] local k_ = `res'[1,4] local T = `res'[1,5] local lagy = `res'[1,6] local minAIC = `res'[1,7] if "`header'"=="" { _rals_fadl_header "`y'" "`x'" `model' `kreg' `T' } _rals_fadl_table `tFADL' `tRALS' `rho2' `cv' `k_' `lagy' `minAIC' `kreg' `model' return scalar tauFADL = `tFADL' return scalar tauRALS = `tRALS' return scalar rho2 = `rho2' return scalar kfreq = `k_' return scalar T = `T' return scalar lag = `lagy' return scalar AIC = `minAIC' return scalar cv01 = `cv'[1,1] return scalar cv05 = `cv'[1,2] return scalar cv10 = `cv'[1,3] return local test = "RALS-FADL" if "`graph'" != "" { _rals_fadl_graph "`y'" "`x'" "`timevar'" "`touse'" `k_' `T' } end program define _rals_fadl_header args y x model kreg T di as text "" di as text "{c TLC}{hline 78}{c TRC}" di as text "{c |} " as result "RALS-Fourier ADL cointegration test{col 79}{c |}" di as text "{c |} Reference : Yilanci, Ulucak, Zhang & Andreoni (2022, Sustainable Dev.){col 79}{c |}" di as text "{c |} Dependent : " as result "`y'" as text "{col 79}{c |}" di as text "{c |} Regressors: " as result "`x'" as text "{col 79}{c |}" di as text "{c |} k = `kreg' T = `T' model = " as result cond(`model'==1,"constant","constant+trend") as text "{col 79}{c |}" di as text "{c BLC}{hline 78}{c BRC}" end program define _rals_fadl_table args tFADL tRALS rho2 cv kfreq lag aic kreg model di as text "" di as text "{hline 78}" di as text " Test Statistic rho^2 1% 5% 10%" di as text "{hline 78}" di as text " Fourier ADL (stage 1) " as result %10.4f `tFADL' /// as text " . . . ." di as text " RALS-FADL (stage 2) " as result %10.4f `tRALS' /// as text " " as result %8.4f `rho2' /// as text " " as result %8.3f `cv'[1,1] " " %8.3f `cv'[1,2] " " %8.3f `cv'[1,3] di as text "{hline 78}" di as text " Optimal Fourier frequency k = " as result `kfreq' /// as text " Δy lag = " as result `lag' /// as text " minAIC = " as result %8.4f `aic' local concl = cond(`tRALS' < `cv'[1,2], "{bf:reject} H0: no cointegration at 5%", /// "fail to reject H0: no cointegration at 5%") di as result " Conclusion: `concl'." di as text "{hline 78}" di as text "" end program define _rals_fadl_graph args y x timevar touse kfreq T local xfirst : word 1 of `x' tempvar fc fs local pi = _pi qui gen `fc' = sin(2*`pi'*_n*`kfreq'/`T')+cos(2*`pi'*_n*`kfreq'/`T') if `touse' twoway (line `y' `timevar' if `touse', lcolor("32 119 180") lwidth(medthick) yaxis(1)) /// (line `xfirst' `timevar' if `touse', lcolor("220 50 50") lpattern(dash) yaxis(1)) /// (line `fc' `timevar' if `touse', lcolor("80 80 80") lpattern(solid) yaxis(2)), /// title("RALS-Fourier ADL: `y' on `xfirst'", size(medium)) /// subtitle("Optimal Fourier frequency k=`kfreq'", size(small)) /// legend(order(1 "`y'" 2 "`xfirst'" 3 "Fourier component") row(1) size(small)) /// ytitle("levels") xtitle("`timevar'") /// scheme(s2color) graphregion(color(white)) name(ralsfadl_`y', replace) end version 14.0 mata: mata set matastrict off void _rals_fadl_engine(string scalar yname, string scalar xnames, string scalar touse, real scalar model, real scalar pmax, real scalar fmax, real scalar fix_k, string scalar rname, string scalar cvname) { real colvector y real matrix X st_view(y, ., yname, touse) st_view(X, ., xnames, touse) real scalar n, k, kk, kbest, popt, j_ n = rows(y) k = cols(X) real colvector dy real matrix dX dy = y[2::n] - y[1::n-1] dX = X[2::n,.] - X[1::n-1,.] real scalar kfrom, kto if (fix_k>0) { kfrom = fix_k kto = fix_k } else { kfrom = 1 kto = fmax } real scalar bestAIC, AICv, p, jx bestAIC = 1e+20; kbest = 1; popt = 1 real colvector best_dep real matrix best_X real colvector best_b, best_e real scalar best_n // Grid over Fourier frequency and lag of Δy (Δx uses same lag for parity // with the Eviews loops in ralsfadl{1..4}.prg). for (kk=kfrom; kk<=kto; kk++) { real colvector ss, cc ss = sin(2*pi()*(2::n)*kk/n) cc = cos(2*pi()*(2::n)*kk/n) for (p=1; p<=pmax; p++) { // build design ---------------------------------------------------- // y on detm + y(-1) + X(-1) + Δx(-1..-p) + Δy(-1..-p) + s + c real matrix Lm_dy Lm_dy = J(rows(dy), p, 0) for (j_=1; j_<=p; j_++) Lm_dy[j_+1::rows(dy), j_] = dy[1::rows(dy)-j_] real matrix Lm_dx Lm_dx = J(rows(dy), p*k, 0) for (jx=1; jx<=k; jx++) { for (j_=1; j_<=p; j_++) { Lm_dx[j_+1::rows(dy), (jx-1)*p + j_] = dX[1::rows(dy)-j_, jx] } } real matrix detm if (model==1) detm = J(rows(dy),1,1) else detm = J(rows(dy),1,1), (2::n) real matrix Xall Xall = detm, y[1::n-1], X[1::n-1,.], Lm_dx, Lm_dy, ss, cc // trim leading p rows real colvector depf depf = dy[p+1::rows(dy)] real matrix Xf Xf = Xall[p+1::rows(Xall), .] // OLS ------------------------------------------------------------ real matrix XX, XXi real colvector bv, ev XX = quadcross(Xf,Xf) XXi = invsym(XX) bv = XXi*quadcross(Xf,depf) ev = depf - Xf*bv real scalar nObs, kobs, AICc nObs = rows(Xf); kobs = cols(Xf) AICc = ln(ev'ev/nObs) + 2*(kobs+2)/nObs if (AICc < bestAIC) { bestAIC = AICc; kbest = kk; popt = p best_dep = depf; best_X = Xf; best_b = bv; best_e = ev; best_n = nObs } } } // stage-1 statistic: t on X(-1)[,1] -> column position is detcols + 1 (y(-1)) + 1 (first X(-1)) real scalar detcols detcols = (model==1 ? 1 : 2) real scalar pos pos = detcols + 1 + 1 // y(-1) then first X(-1) --> Stata uses the same convention // We compute t on the X(-1)[,1] regressor (the first regressor's level lag) // as the EViews code does: @tstats(2) where position 2 is x(-1) under // a constant-only model. In general (constant + trend), x(-1) is at // position 3 (counting from 1 in the design matrix). // pos already encodes this. real matrix XXi_ XXi_ = invsym(quadcross(best_X,best_X)) real colvector se real scalar sig2 sig2 = (best_e'best_e)/(rows(best_X)-cols(best_X)) se = sqrt(diagonal(sig2*XXi_)) real scalar tFADL tFADL = best_b[pos,1] / se[pos,1] // RALS augmentation real scalar m2v, m3v, sig2A m2v = sum(best_e:^2)/rows(best_e) m3v = sum(best_e:^3)/rows(best_e) real matrix W, XR, XRi real colvector bR, eR, sR W = (best_e:^2 :- m2v), (best_e:^3 :- m3v :- 3*m2v*best_e) XR = best_X, W XRi = invsym(quadcross(XR,XR)) bR = XRi*quadcross(XR,best_dep) eR = best_dep - XR*bR sig2A = (eR'eR)/(rows(XR)-cols(XR)) sR = sqrt(diagonal(sig2A*XRi)) real scalar tRALS, rho2 tRALS = bR[pos,1]/sR[pos,1] rho2 = sig2A/sig2 if (rho2<0) rho2=0 if (rho2>1) rho2=1 // CV: use RALS-Fourier ADF table parameterised by # of regressors real rowvector cvr cvr = __rals_fadf_cv(n-popt, k+1, model, rho2) st_matrix(rname, (tFADL, tRALS, rho2, kbest, n, popt, bestAIC)) st_matrix(cvname, cvr) } end