*! ralsfadf 1.0.1 16may2026 Dr Merwan Roudane *! RALS-Fourier ADF unit-root test (Yilanci, Aydin & Aydin 2019) * MPRA Paper No. 96797. Critical values from Tables 1a / 1b of the paper, * Hansen (1995) interpolation across rho^2 and bilinear in T. *------------------------------------------------------------------------------ program define ralsfadf, rclass version 14.0 capture mata: __rals_loaded() if _rc qui _rals_mata syntax varname(ts) [if] [in], [ /// TREND /// MAXLags(integer 8) /// IC(string) /// 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 } if `frequency' != -1 & (`frequency'<1 | `frequency'>5) { di as error "frequency() must be between 1 and 5" exit 198 } marksample touse qui count if `touse' if r(N) < 25 { di as error "sample too small for RALS-FADF (N=" r(N) ")" exit 2001 } capture tsset if _rc { di as error "data must be {bf:tsset} before running ralsfadf" exit 459 } local timevar `r(timevar)' if "`ic'"=="" local ic "tstat" local ic = lower("`ic'") if !inlist("`ic'","aic","bic","tstat") { di as error "ic() must be aic, bic, or tstat" exit 198 } local model = cond("`trend'"!="", 2, 1) tempname res cv mata: _rals_fadf_engine("`varlist'", "`touse'", `model', `maxlags', "`ic'", `fmax', `frequency', "`res'", "`cv'") local tauFADF = `res'[1,1] local tauRALS = `res'[1,2] local rho2 = `res'[1,3] local T = `res'[1,4] local lag = `res'[1,5] local kfreq = `res'[1,6] local ssr = `res'[1,7] if "`header'"=="" { _rals_print_header "RALS-FADF" "Yilanci, Aydin & Aydin (2019, MPRA 96797)" "`varlist'" `T' `lag' "`ic'" `model' } _rals_print_fadf_table `tauFADF' `tauRALS' `rho2' `cv' `kfreq' `ssr' return scalar tauFADF = `tauFADF' return scalar tauRALS = `tauRALS' return scalar rho2 = `rho2' return scalar T = `T' return scalar lag = `lag' return scalar kfreq = `kfreq' return scalar ssr = `ssr' return scalar cv01_FADF = `cv'[1,1] return scalar cv05_FADF = `cv'[1,2] return scalar cv10_FADF = `cv'[1,3] return scalar cv01 = `cv'[2,1] return scalar cv05 = `cv'[2,2] return scalar cv10 = `cv'[2,3] return local test = "RALS-FADF" return local cmdline = "ralsfadf `varlist', `trend' maxlags(`maxlags') ic(`ic') fmax(`fmax')" if "`graph'" != "" { _rals_graph_fadf "`varlist'" "`timevar'" "`touse'" `tauRALS' `cv'[2,2] `kfreq' `T' } end program define _rals_print_fadf_table args tauFADF tauRALS rho2 cv kfreq ssr _rals_print_two_row "Fourier ADF (stage 1)" "RALS-FADF (stage 2)" /// `tauFADF' `tauRALS' `rho2' /// `cv'[1,1] `cv'[1,2] `cv'[1,3] `cv'[2,1] `cv'[2,2] `cv'[2,3] di as text " Optimal Fourier frequency k = " as result `kfreq' /// as text " SSR minimum = " as result %10.4f `ssr' local concl_f = cond(`tauFADF' < `cv'[1,2], "{bf:reject} H0 (FADF, 5%)", /// "fail to reject H0 (FADF, 5%)") local concl_r = cond(`tauRALS' < `cv'[2,2], "{bf:reject} H0 (RALS-FADF, 5%)", /// "fail to reject H0 (RALS-FADF, 5%)") di as text " Decision (5%)" di as text " {c -}{c -} FADF : " as result "`concl_f'" di as text " {c -}{c -} RALS-FADF : " as result "`concl_r'" di as text "{hline 80}" di as text " Stage-1 CVs : Enders & Lee (2012) Fourier-ADF table." di as text " Stage-2 CVs : Yilanci, Aydin & Aydin (2019) RALS-FADF table." di as text "" end * _rals_graph_fadf is defined in its own .ado (shared with ralsfkss). version 14.0 mata: mata set matastrict off void _rals_fadf_engine(string scalar yname, string scalar touse, real scalar model, real scalar pmax, string scalar icmode, real scalar fmax, real scalar fix_k, string scalar rname, string scalar cvname) { real colvector y, dy, ly, t, s, c real scalar T, p, sig2, sig2A, tauFADF, tauRALS, m2, m3, rho2 real scalar kbest, ssrmin, popt, bestic, j, k real matrix lmat st_view(y, ., yname, touse) T = rows(y) dy = y[2::T] - y[1::T-1] ly = y[1::T-1] t = (2::T) real scalar Tobs Tobs = T // pre-build lag matrix of dy lmat = J(rows(dy), pmax, 0) for (j=1; j<=pmax; j++) lmat[j+1::rows(dy), j] = dy[1::rows(dy)-j] // Step A: grid search for optimal Fourier frequency using OLS SSR -------- ssrmin = 1e+20; kbest = 1 real scalar kfrom, kto if (fix_k>0) { kfrom = fix_k kto = fix_k } else { kfrom = 1 kto = fmax } real colvector pi_ pi_ = J(1,1,pi()) // unused; kept for parity real scalar p_force p_force = pmax for (k=kfrom; k<=kto; k++) { s = sin(2*pi()*(2::T)*k/Tobs) c = cos(2*pi()*(2::T)*k/Tobs) real matrix Xfull if (model==1) Xfull = ly, J(rows(ly),1,1), s, c, lmat else Xfull = ly, J(rows(ly),1,1), (2::T), s, c, lmat real scalar ndet ndet = (model==1 ? 4 : 5) // y(-1) + const + s + c [+ trend] real colvector depf depf = dy real matrix Xf Xf = Xfull if (pmax>0) { depf = depf[pmax+1::rows(dy)] Xf = Xf[pmax+1::rows(dy), .] } // pick lag via icmode ---------------------------------------------- real scalar popt_k real colvector icvec, tvec icvec = J(pmax+1,1,.); tvec = J(pmax+1,1,.) for (p=0; p<=pmax; p++) { real matrix Xp,XXi real colvector bp,ep,sep Xp = Xf[., 1::(ndet+p)] XXi = invsym(quadcross(Xp,Xp)) bp = XXi*quadcross(Xp,depf) ep = depf - Xp*bp if (icmode=="bic") { icvec[p+1] = ln(ep'ep/rows(Xp)) + cols(Xp)*ln(rows(Xp))/rows(Xp) } else { icvec[p+1] = ln(ep'ep/rows(Xp)) + 2*cols(Xp)/rows(Xp) } if (p>0) { sep = sqrt(diagonal((ep'ep/(rows(Xp)-cols(Xp)))*XXi)) tvec[p+1] = abs(bp[cols(Xp),1]/sep[cols(Xp),1]) } } if (icmode=="tstat") { popt_k = 0 for (p=pmax; p>=1; p--) { if (tvec[p+1] >= 1.645) { popt_k = p break } } } else { popt_k = 0 bestic = icvec[1] for (p=1; p<=pmax; p++) { if (icvec[p+1] < bestic) { bestic = icvec[p+1] popt_k = p } } } real matrix Xs Xs = Xf[., 1::(ndet+popt_k)] real matrix XXi_s real colvector bs_, es_ XXi_s = invsym(quadcross(Xs,Xs)) bs_ = XXi_s*quadcross(Xs,depf) es_ = depf - Xs*bs_ real scalar ssrv ssrv = es_'es_ if (ssrv < ssrmin) { ssrmin = ssrv kbest = k popt = popt_k } } // Step B: at kbest fit the final FADF and RALS extension ------------------ s = sin(2*pi()*(2::T)*kbest/Tobs) c = cos(2*pi()*(2::T)*kbest/Tobs) real matrix Xfull2, Xf2 if (model==1) Xfull2 = ly, J(rows(ly),1,1), s, c, lmat else Xfull2 = ly, J(rows(ly),1,1), (2::T), s, c, lmat real scalar ndet2 ndet2 = (model==1 ? 4 : 5) real colvector depf2 depf2 = dy; Xf2 = Xfull2 if (pmax>0) { depf2 = depf2[pmax+1::rows(dy)] Xf2 = Xf2[pmax+1::rows(dy), .] } real matrix Xs2, XXi2 real colvector bs2, es2, sebs2 Xs2 = Xf2[., 1::(ndet2+popt)] XXi2 = invsym(quadcross(Xs2,Xs2)) bs2 = XXi2*quadcross(Xs2,depf2) es2 = depf2 - Xs2*bs2 sig2 = (es2'es2)/(rows(Xs2)-cols(Xs2)) sebs2 = sqrt(diagonal(sig2*XXi2)) tauFADF = bs2[1,1]/sebs2[1,1] // RALS augmentation m2 = sum(es2:^2)/rows(es2) m3 = sum(es2:^3)/rows(es2) real matrix W, Xa, XaI real colvector ba, ea, sea W = (es2:^2 :- m2), (es2:^3 :- m3 :- 3*m2*es2) Xa = Xs2, W XaI = invsym(quadcross(Xa,Xa)) ba = XaI*quadcross(Xa,depf2) ea = depf2 - Xa*ba sig2A = (ea'ea)/(rows(Xa)-cols(Xa)) sea = sqrt(diagonal(sig2A*XaI)) tauRALS = ba[1,1]/sea[1,1] rho2 = sig2A/sig2 if (rho2<0) rho2=0 if (rho2>1) rho2=1 // CV lookup (n approx = T-popt, k=1, since 1 regressor i.e. y itself) ---- real rowvector cvr, cvFADF cvr = __rals_fadf_cv(T-popt, 1, model, rho2) cvFADF = __rals_cv_fadf(T-popt, kbest, model) real matrix CV2 CV2 = cvFADF \ cvr st_matrix(rname, (tauFADF, tauRALS, rho2, T, popt, kbest, ssrmin)) st_matrix(cvname, CV2) } end