*! ralscoint 1.0.1 16may2026 Dr Merwan Roudane *! RALS cointegration tests: ECM, ADL, EG, EG2 with RALS augmentation * Reference: Lee, H., Lee, J., & Im, K. (2015) Studies in Nonlinear Dynamics * & Econometrics 19(4): 397-413. * Mirrors the GAUSS code RALS_coint_size_power.g and RALS_coint_crit.g from * the SNDE 2013-0060 supplement (Hyejin Lee 2012 dissertation, Ch. 3). *------------------------------------------------------------------------------ program define ralscoint, rclass version 14.0 capture mata: __rals_loaded() if _rc qui _rals_mata syntax varlist(min=2 ts) [if] [in], [ /// TREND /// METHod(string) /// ecm adl eg eg2 all BETA(real 1) /// prespecified beta (ECM only) BW(integer 0) /// LR-var bandwidth ; 0 = automatic Graph /// noHEADer ] if "`method'"=="" local method "all" local method = lower("`method'") if !inlist("`method'","ecm","adl","eg","eg2","all") { di as error "method() must be ecm, adl, eg, eg2 or all" 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 ralscoint" exit 459 } local timevar `r(timevar)' local y : word 1 of `varlist' local x : list varlist - y local kreg : list sizeof x local model = cond("`trend'"!="", 2, 1) tempname res cv mata: _rals_coint_engine("`y'", "`x'", "`touse'", `model', `bw', `beta', "`method'", "`res'", "`cv'") matrix colnames `res' = ECM_t ECM_rals ECM_rho2 ADL_t ADL_rals ADL_rho2 EG_t EG_rals EG_rho2 EG2_t EG2_rals EG2_rho2 T if "`header'"=="" { _rals_coint_header "`y'" "`x'" `model' `kreg' `bw' } _rals_coint_table `res' `cv' "`method'" `kreg' `model' return scalar T = `res'[1,13] return scalar ECM_t = `res'[1,1] return scalar ECM_rals = `res'[1,2] return scalar ECM_rho2 = `res'[1,3] return scalar ADL_t = `res'[1,4] return scalar ADL_rals = `res'[1,5] return scalar ADL_rho2 = `res'[1,6] return scalar EG_t = `res'[1,7] return scalar EG_rals = `res'[1,8] return scalar EG_rho2 = `res'[1,9] return scalar EG2_t = `res'[1,10] return scalar EG2_rals = `res'[1,11] return scalar EG2_rho2 = `res'[1,12] return matrix cv = `cv' return local yvar = "`y'" return local xvars = "`x'" return local model = cond(`model'==1,"constant","constant+trend") return local method = "`method'" if "`graph'" != "" { _rals_coint_graph "`y'" "`x'" "`timevar'" "`touse'" } end program define _rals_coint_header args y x model kreg bw di as text "" di as text "{c TLC}{hline 78}{c TRC}" di as text "{c |} " as result "RALS cointegration tests (ECM/ADL/EG/EG2){col 79}{c |}" di as text "{c |} Reference : Lee, Lee & Im (2015, SNDE 19(4):397-413){col 79}{c |}" di as text "{c |} Dependent : " as result "`y'{col 79}{c |}" di as text "{c |} Regressors: " as result "`x'{col 79}{c |}" di as text "{c |} k = `kreg' model = " as result cond(`model'==1,"constant","constant+trend") /// as text " LRV bandwidth = " as result cond(`bw'==0,"auto","`bw'") as text "{col 79}{c |}" di as text "{c BLC}{hline 78}{c BRC}" end program define _rals_coint_table args r cv method kreg model di as text "" di as text "{hline 78}" di as text " Test Stage1 t RALS stat rho^2 OLS 5% CV RALS 1%/5%/10% CVs" di as text "{hline 78}" local rows "ECM ADL EG EG2" local i 1 foreach m of local rows { local print = 0 if "`method'"=="all" local print 1 if "`method'"=="`=lower("`m'")'" local print 1 if `print' { local t1 = `r'[1, (`i'-1)*3 + 1] local t2 = `r'[1, (`i'-1)*3 + 2] local r2 = `r'[1, (`i'-1)*3 + 3] local olscv = `cv'[`i', 1] local rcv1 = `cv'[`i', 2] local rcv5 = `cv'[`i', 3] local rcv10 = `cv'[`i', 4] di as text " `m' " as result %8.3f `t1' " " %8.3f `t2' /// as text " " as result %6.3f `r2' /// as text " " as result %8.3f `olscv' /// as text " " as result %7.3f `rcv1' " " %7.3f `rcv5' " " %7.3f `rcv10' } local ++i } di as text "{hline 78}" di as text " Reject H0 of NO cointegration when statistic < CV(5%)." di as text "{hline 78}" end program define _rals_coint_graph args y x timevar touse local xfirst : word 1 of `x' twoway (line `y' `timevar' if `touse', lcolor("32 119 180") lwidth(medthick)) /// (line `xfirst' `timevar' if `touse', lcolor("220 50 50") lpattern(dash)),/// title("Cointegration: `y' and `xfirst'", size(medium)) /// legend(order(1 "`y'" 2 "`xfirst'") row(1) size(small)) /// scheme(s2color) graphregion(color(white)) name(ralscoint_`y', replace) end version 14.0 mata: mata set matastrict off void _rals_coint_engine(string scalar yname, string scalar xnames, string scalar touse, real scalar model, real scalar bw, real scalar beta_pre, string scalar method, 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, p n = rows(y); k = cols(X) p = model - 1 real matrix detm, detm1 if (model==1) { detm = J(n-1,1,1) detm1 = J(n,1,1) } else { detm = J(n-1,1,1), (1::n-1) detm1 = J(n,1,1), (1::n) } // pre-compute helper quantities ------------------------------------------ real colvector dy real matrix dx dy = y[2::n] - y[1::n-1] dx = X[2::n,.] - X[1::n-1,.] // -- 1. ECM (prespecified beta) ----------------------------------------- real colvector zECM zECM = y - X*J(k,1,beta_pre) real colvector zECM1, dECM zECM1 = zECM[1::n-1] dECM = zECM[2::n] - zECM[1::n-1] real matrix ecmX, ecmXX if (model==1) ecmX = zECM1, detm else ecmX = zECM1, detm ecmXX = ecmX, dx // OLS for RALS auxiliary residuals (without dx, see GAUSS code) real matrix XA, XAi real colvector bA, eA, sA XA = ecmX XAi = invsym(quadcross(XA,XA)) bA = XAi*quadcross(XA,dy) eA = dy - XA*bA // OLS for stage-1 ECM stat (with dx) real matrix XB, XBi real colvector bB, eB, sB XB = ecmXX XBi = invsym(quadcross(XB,XB)) bB = XBi*quadcross(XB,dy) eB = dy - XB*bB real scalar sigB sigB = (eB'eB)/(rows(XB)-cols(XB)) sB = sqrt(diagonal(sigB*XBi)) real scalar ECMt ECMt = bB[1,1]/sB[1,1] // RALS augmentation on eA real matrix W_ecm, ECMR, ECMRi real colvector b_ecmR, e_ecmR, s_ecmR real scalar m2v, m3v, sigA m2v = sum(eA:^2)/rows(eA); m3v = sum(eA:^3)/rows(eA) W_ecm = (eA:^2 :- m2v), (eA:^3 :- m3v :- 3*m2v*eA) ECMR = ecmXX, W_ecm ECMRi = invsym(quadcross(ECMR,ECMR)) b_ecmR = ECMRi*quadcross(ECMR,dy) e_ecmR = dy - ECMR*b_ecmR real scalar sigAR sigAR = (e_ecmR'e_ecmR)/(rows(ECMR)-cols(ECMR)) s_ecmR = sqrt(diagonal(sigAR*ECMRi)) real scalar ECMt_rals, rho_ecm ECMt_rals = b_ecmR[1,1]/s_ecmR[1,1] rho_ecm = __rals_rho2_lr(e_ecmR, eA) // -- 2. ADL -------------------------------------------------------------- real matrix ADLx, ADLxi real colvector b_adl, e_adl, s_adl if (model==1) ADLx = y[1::n-1], X[1::n-1,.], dx else ADLx = y[1::n-1], X[1::n-1,.], dx, detm ADLxi = invsym(quadcross(ADLx,ADLx)) b_adl = ADLxi*quadcross(ADLx,dy) e_adl = dy - ADLx*b_adl real scalar sig_adl, ADLt sig_adl = (e_adl'e_adl)/(rows(ADLx)-cols(ADLx)) s_adl = sqrt(diagonal(sig_adl*ADLxi)) ADLt = b_adl[1,1]/s_adl[1,1] // RALS ---- real matrix W_adl, ADLR, ADLRi real colvector b_adlR, e_adlR, s_adlR real scalar m2a, m3a m2a = sum(e_adl:^2)/rows(e_adl); m3a = sum(e_adl:^3)/rows(e_adl) W_adl = (e_adl:^2 :- m2a), (e_adl:^3 :- m3a :- 3*m2a*e_adl) ADLR = ADLx, W_adl ADLRi = invsym(quadcross(ADLR,ADLR)) b_adlR = ADLRi*quadcross(ADLR,dy) e_adlR = dy - ADLR*b_adlR real scalar sig_adlR sig_adlR = (e_adlR'e_adlR)/(rows(ADLR)-cols(ADLR)) s_adlR = sqrt(diagonal(sig_adlR*ADLRi)) real scalar ADLt_rals, rho_adl ADLt_rals = b_adlR[1,1]/s_adlR[1,1] rho_adl = __rals_rho2_lr(e_adlR, e_adl) // -- 3. EG --------------------------------------------------------------- real matrix Xa_eg, Xa_egi real colvector b_eg, z_eg if (model==1) Xa_eg = X else Xa_eg = X, detm1 Xa_egi = invsym(quadcross(Xa_eg,Xa_eg)) b_eg = Xa_egi*quadcross(Xa_eg,y) z_eg = y - Xa_eg*b_eg real colvector d_eg, eg1 d_eg = z_eg[2::n] - z_eg[1::n-1] eg1 = z_eg[1::n-1] real matrix Xeg, Xegi real colvector b_eg2, e_eg2, s_eg2 Xeg = eg1 Xegi = invsym(quadcross(Xeg,Xeg)) b_eg2 = Xegi*quadcross(Xeg,d_eg) e_eg2 = d_eg - Xeg*b_eg2 real scalar sig_eg, EGt sig_eg = (e_eg2'e_eg2)/(rows(Xeg)-cols(Xeg)) s_eg2 = sqrt(diagonal(sig_eg*Xegi)) EGt = b_eg2[1,1]/s_eg2[1,1] // RALS real scalar m2e, m3e m2e = sum(e_eg2:^2)/rows(e_eg2); m3e = sum(e_eg2:^3)/rows(e_eg2) real matrix W_eg, EGR, EGRi W_eg = (e_eg2:^2 :- m2e), (e_eg2:^3 :- m3e :- 3*m2e*e_eg2) if (model==1) { EGR = eg1, W_eg } else { EGR = eg1, W_eg, detm } EGRi = invsym(quadcross(EGR,EGR)) real colvector b_egR, e_egR, s_egR b_egR = EGRi*quadcross(EGR,d_eg) e_egR = d_eg - EGR*b_egR real scalar sig_egR sig_egR = (e_egR'e_egR)/(rows(EGR)-cols(EGR)) s_egR = sqrt(diagonal(sig_egR*EGRi)) real scalar EGt_rals, rho_eg EGt_rals = b_egR[1,1]/s_egR[1,1] rho_eg = __rals_rho2_lr(e_egR, e_eg2) // -- 4. EG2 (Lee 2012) -------------------------------------------------- real matrix EG2x, EG2xi real colvector b_eg22, e_eg22, s_eg22 if (model==1) EG2x = eg1, dx else EG2x = eg1, dx, detm EG2xi = invsym(quadcross(EG2x,EG2x)) b_eg22 = EG2xi*quadcross(EG2x,d_eg) e_eg22 = d_eg - EG2x*b_eg22 real scalar sig_eg22, EG2t sig_eg22 = (e_eg22'e_eg22)/(rows(EG2x)-cols(EG2x)) s_eg22 = sqrt(diagonal(sig_eg22*EG2xi)) EG2t = b_eg22[1,1]/s_eg22[1,1] // RALS on EG2: auxiliary residuals come from regression of d_eg on eg1[+detm] (per GAUSS) real matrix Xaux if (model==1) Xaux = eg1 else Xaux = eg1, detm real matrix Xauxi Xauxi = invsym(quadcross(Xaux,Xaux)) real colvector b_aux, e_aux b_aux = Xauxi*quadcross(Xaux,d_eg) e_aux = d_eg - Xaux*b_aux real scalar m2g, m3g m2g = sum(e_aux:^2)/rows(e_aux); m3g = sum(e_aux:^3)/rows(e_aux) real matrix W_eg2, EG2R, EG2Ri W_eg2 = (e_aux:^2 :- m2g), (e_aux:^3 :- m3g :- 3*m2g*e_aux) EG2R = EG2x, W_eg2 EG2Ri = invsym(quadcross(EG2R,EG2R)) real colvector b_eg2R, e_eg2R, s_eg2R b_eg2R = EG2Ri*quadcross(EG2R,d_eg) e_eg2R = d_eg - EG2R*b_eg2R real scalar sig_eg2R sig_eg2R = (e_eg2R'e_eg2R)/(rows(EG2R)-cols(EG2R)) s_eg2R = sqrt(diagonal(sig_eg2R*EG2Ri)) real scalar EG2t_rals, rho_eg2 EG2t_rals = b_eg2R[1,1]/s_eg2R[1,1] rho_eg2 = __rals_rho2_lr(e_eg2R, e_aux) // ---- Critical values --------------------------------------------------- // Each row of cvr is : [OLS-5% , RALS-1% , RALS-5% , RALS-10%] real matrix cvr cvr = J(4,4,.) cvr[1,1] = __rals_cv_ols(1, model, k) cvr[2,1] = __rals_cv_ols(2, model, k) cvr[3,1] = __rals_cv_ols(3, model, k) cvr[4,1] = __rals_cv_ols(4, model, k) real rowvector c1, c2, c3, c4 c1 = __rals_interp(__rals_cv_ecm(p), rho_ecm) c2 = __rals_interp(__rals_cv_adl(p), rho_adl) c3 = __rals_interp(__rals_cv_eg(p), rho_eg) c4 = __rals_interp(__rals_cv_eg(p), rho_eg2) cvr[1,2..4] = c1 cvr[2,2..4] = c2 cvr[3,2..4] = c3 cvr[4,2..4] = c4 st_matrix(rname, (ECMt, ECMt_rals, rho_ecm, ADLt, ADLt_rals, rho_adl, EGt, EGt_rals, rho_eg, EG2t, EG2t_rals, rho_eg2, n)) st_matrix(cvname, cvr) } end