*! _mc_mata 1.0.2 18may2026 Dr Merwan Roudane (merwanroudane920@gmail.com) *! Mata core for the multicoint package -- pure .mata source file. *! *! Loaded on demand by every multicoint command via *! qui capture mata: __mc_loaded() *! if _rc qui _mc_mata version 14.0 mata mata set matastrict off void __mc_loaded() { /* sentinel */ } // ---------------------------------------------------------------------- // ENGSTED-GONZALO-HALDRUP (1997) critical-value tables // ---------------------------------------------------------------------- real matrix _mc_egh_table(string scalar trend, real scalar m2) { real matrix tab tab = J(5, 20, 0) if (trend == "ct" & m2 == 1) { tab[1,] = (-5.21,-4.72,-4.29,-3.88,-4.66,-4.33,-4.01,-3.67,-4.55,-4.18,-3.90,-3.59,-4.41,-4.04,-3.83,-3.51,-4.33,-4.04,-3.78,-3.49) tab[2,] = (-5.60,-5.10,-4.71,-4.30,-5.11,-4.70,-4.42,-4.08,-4.85,-4.54,-4.26,-3.94,-4.73,-4.43,-4.19,-3.89,-4.73,-4.42,-4.15,-3.87) tab[3,] = (-6.09,-5.57,-5.14,-4.69,-5.47,-5.07,-4.74,-4.38,-5.21,-4.86,-4.58,-4.26,-5.07,-4.79,-4.51,-4.20,-5.00,-4.73,-4.48,-4.18) tab[4,] = (-6.47,-5.95,-5.53,-5.08,-5.89,-5.43,-5.13,-4.76,-5.52,-5.18,-4.91,-4.59,-5.38,-5.05,-4.78,-4.74,-5.34,-5.04,-4.78,-4.50) tab[5,] = (-6.95,-6.37,-5.90,-5.44,-6.35,-5.85,-5.47,-5.10,-5.86,-5.49,-5.20,-4.89,-5.66,-5.35,-5.08,-4.77,-5.63,-5.31,-5.06,-4.76) return(tab) } if (trend == "ct" & m2 == 2) { tab[1,] = (-5.81,-5.25,-4.83,-4.41,-5.14,-4.77,-4.45,-4.10,-4.93,-4.56,-4.31,-3.98,-4.81,-4.49,-4.20,-3.91,-4.75,-4.42,-4.14,-3.84) tab[2,] = (-6.24,-5.68,-5.21,-4.80,-5.62,-5.22,-4.89,-4.51,-5.23,-4.90,-4.62,-4.29,-5.11,-4.77,-4.50,-4.20,-5.05,-4.74,-4.48,-4.18) tab[3,] = (-6.70,-6.17,-5.70,-5.22,-5.98,-5.53,-5.17,-4.79,-5.59,-5.19,-4.93,-4.62,-5.35,-5.07,-4.80,-4.51,-5.34,-5.02,-4.75,-4.46) tab[4,] = (-7.19,-6.63,-6.08,-5.89,-6.23,-5.81,-5.48,-5.12,-5.97,-5.58,-5.25,-4.92,-5.69,-5.37,-5.07,-4.80,-5.67,-5.33,-5.06,-4.76) tab[5,] = (-7.61,-6.93,-6.43,-5.91,-6.64,-6.18,-5.82,-5.41,-6.09,-5.76,-5.50,-5.16,-5.95,-5.61,-5.34,-5.04,-5.92,-5.56,-5.29,-5.02) return(tab) } if (trend == "ctt" & m2 == 1) { tab[1,] = (-5.77,-5.28,-4.86,-4.43,-5.20,-4.81,-4.47,-4.12,-4.94,-4.60,-4.32,-3.98,-4.77,-4.47,-4.21,-3.92,-4.73,-4.43,-4.17,-3.88) tab[2,] = (-6.21,-5.69,-5.27,-4.83,-5.56,-5.16,-4.83,-4.47,-5.29,-4.93,-4.64,-4.32,-5.11,-4.79,-4.52,-4.20,-5.05,-4.75,-4.49,-4.20) tab[3,] = (-6.66,-6.10,-5.65,-5.20,-5.92,-5.50,-5.17,-4.82,-5.57,-5.23,-4.95,-4.63,-5.42,-5.08,-4.82,-4.52,-5.36,-5.04,-4.77,-4.48) tab[4,] = (-7.12,-6.51,-6.05,-5.55,-6.27,-5.85,-5.50,-5.12,-5.90,-5.54,-5.25,-4.91,-5.71,-5.38,-5.11,-4.81,-5.60,-5.30,-5.04,-4.76) tab[5,] = (-7.61,-6.93,-6.43,-5.91,-6.56,-6.15,-5.79,-5.41,-6.18,-5.81,-5.52,-5.19,-5.96,-5.64,-5.36,-5.05,-5.87,-5.57,-5.30,-5.26) return(tab) } tab[1,] = (-6.44,-5.85,-5.42,-4.96,-5.61,-5.21,-4.88,-4.52,-5.33,-4.97,-4.67,-4.34,-5.13,-4.79,-4.52,-4.23,-5.07,-4.76,-4.50,-4.21) tab[2,] = (-6.85,-6.30,-5.82,-5.33,-5.99,-5.58,-5.22,-4.86,-5.63,-5.27,-4.98,-4.65,-5.43,-5.09,-4.84,-4.54,-5.35,-5.05,-4.78,-4.49) tab[3,] = (-7.32,-6.68,-6.21,-5.69,-6.35,-5.90,-5.54,-5.16,-5.90,-5.54,-5.25,-4.92,-5.69,-5.37,-5.10,-4.80,-5.61,-5.29,-5.04,-4.76) tab[4,] = (-7.68,-7.06,-6.55,-6.03,-6.63,-6.23,-5.86,-5.46,-6.19,-5.85,-5.55,-5.22,-5.96,-5.64,-5.37,-5.07,-5.85,-5.55,-5.30,-5.02) tab[5,] = (-8.18,-7.47,-6.93,-6.38,-7.00,-6.55,-6.16,-5.76,-6.47,-6.10,-5.80,-5.47,-6.21,-5.87,-5.60,-5.31,-6.12,-5.80,-5.54,-5.26) return(tab) } void _mc_egh_cv(string scalar trend, real scalar m1, real scalar m2, real scalar T, string scalar cv01, string scalar cv025, string scalar cv05, string scalar cv10) { real matrix tab real scalar irow, T0, w, Tlo, Thi, v01, v025, v05, v10, mm1, mm2 real rowvector Tvec, qrow mm2 = m2 mm1 = m1 if (mm2 < 1) mm2 = 1 if (mm2 > 2) mm2 = 2 if (mm1 < 0) mm1 = 0 if (mm1 > 4) mm1 = 4 tab = _mc_egh_table(trend, mm2) irow = mm1 + 1 Tvec = (25, 50, 100, 250, 500) qrow = tab[irow,] if (T <= 25) T0 = 1 else if (T <= 50) T0 = 2 else if (T <= 100) T0 = 3 else if (T <= 250) T0 = 4 else T0 = 5 if (T0 == 1) { v01 = qrow[1] v025 = qrow[2] v05 = qrow[3] v10 = qrow[4] } else { Tlo = Tvec[T0-1] Thi = Tvec[T0] w = (T - Tlo)/(Thi - Tlo) if (w < 0) w = 0 if (w > 1) w = 1 v01 = (1-w)*qrow[(T0-2)*4 + 1] + w*qrow[(T0-1)*4 + 1] v025 = (1-w)*qrow[(T0-2)*4 + 2] + w*qrow[(T0-1)*4 + 2] v05 = (1-w)*qrow[(T0-2)*4 + 3] + w*qrow[(T0-1)*4 + 3] v10 = (1-w)*qrow[(T0-2)*4 + 4] + w*qrow[(T0-1)*4 + 4] } st_numscalar(cv01, v01) st_numscalar(cv025, v025) st_numscalar(cv05, v05) st_numscalar(cv10, v10) } // ---------------------------------------------------------------------- // GRANGER-LEE rough p-value (MacKinnon Engle-Granger c.v. surrogate) // ---------------------------------------------------------------------- void _mc_gl_pval(real scalar t, string scalar trend, real scalar k, real scalar T, string scalar pvname, string scalar cv05name) { real matrix asymp05 real scalar trow, kc, cv5, z, pv asymp05 = J(4, 7, 0) asymp05[1,] = (-1.95, -3.34, -3.74, -4.10, -4.42, -4.71, -4.98) asymp05[2,] = (-2.86, -3.34, -3.74, -4.10, -4.42, -4.71, -4.98) asymp05[3,] = (-3.41, -3.78, -4.12, -4.43, -4.72, -4.99, -5.25) asymp05[4,] = (-3.96, -4.32, -4.62, -4.89, -5.14, -5.38, -5.61) if (trend == "none") trow = 1 else if (trend == "c") trow = 2 else if (trend == "ct") trow = 3 else trow = 4 kc = k + 1 if (kc > 7) kc = 7 if (kc < 1) kc = 1 cv5 = asymp05[trow, kc] z = (t - cv5) / 0.8 pv = normal(z) if (pv < 0) pv = 0 if (pv > 1) pv = 1 st_numscalar(pvname, pv) st_numscalar(cv05name, cv5) } // ---------------------------------------------------------------------- // TAOLS estimator -- basis-transform OLS // ---------------------------------------------------------------------- void _mc_taols_compute(string scalar yvar, string scalar Xcums, string scalar Xflows, string scalar Dxs, string scalar trend, real scalar K, string scalar touse, string scalar bname, string scalar Vname, string scalar rssname, string scalar r2name, string scalar Nname, string scalar resvar) { real matrix Y, Xc, Xf, Dx real matrix V_Y, V_Xc, V_Xf, V_Dx, V_l, V_tr real matrix Phi, Xdesign, XX, bvec, resid, Vmat real matrix bout, Vout, bxc, bxf, bdet real colvector uhat, touseflag, full real rowvector keep real scalar T, i, t, nXc, nXf, nDx, nd, sigma2, RSS, TSS, R2, ybar, c, ti, ki Y = st_data(., yvar, touse) Xc = st_data(., tokens(Xcums), touse) Xf = st_data(., tokens(Xflows), touse) Dx = st_data(., tokens(Dxs), touse) T = rows(Y) Phi = J(T, K, 0) for (i=1; i<=K; i++) { for (t=1; t<=T; t++) { Phi[t,i] = sqrt(2)*sin( (i - 0.5)*pi()*(t/T) ) } } V_Y = (Phi' * Y) / sqrt(T) V_Xc = (Phi' * Xc) / sqrt(T) V_Xf = (Phi' * Xf) / sqrt(T) V_Dx = (Phi' * Dx) / sqrt(T) V_l = colsum(Phi)' / sqrt(T) V_tr = J(K, 0, 0) if (trend == "ct" | trend == "ctt") V_tr = (V_tr, (Phi'*(1::T))/sqrt(T)/T) if (trend == "ctt") V_tr = (V_tr, (Phi'*((1::T):^2))/sqrt(T)/T^2) nXc = cols(V_Xc) nXf = cols(V_Xf) nDx = cols(V_Dx) Xdesign = V_Xc, V_Xf, V_Dx if (trend != "none") Xdesign = Xdesign, V_l if (trend == "ct" | trend == "ctt") Xdesign = Xdesign, V_tr[,1] if (trend == "ctt") Xdesign = Xdesign, V_tr[,2] XX = quadcross(Xdesign, Xdesign) bvec = qrsolve(XX, quadcross(Xdesign, V_Y)) resid = V_Y - Xdesign * bvec sigma2 = (resid' * resid) / (K - cols(Xdesign)) if (sigma2 <= 0 | sigma2 == .) sigma2 = (resid' * resid)/K Vmat = sigma2 * invsym(XX) nd = 0 if (trend != "none") nd = nd + 1 if (trend == "ct" | trend == "ctt") nd = nd + 1 if (trend == "ctt") nd = nd + 1 keep = 1..(nXc + nXf) if (nd > 0) keep = (keep, (nXc+nXf+nDx+1)..(nXc+nXf+nDx+nd)) bout = bvec[keep,]' Vout = Vmat[keep,keep] bxc = bvec[1..nXc,] bxf = bvec[(nXc+1)..(nXc+nXf),] uhat = Y - Xc*bxc - Xf*bxf if (nd > 0) { bdet = bvec[(nXc+nXf+nDx+1)..(nXc+nXf+nDx+nd),] c = 1 if (trend != "none") { uhat = uhat :- bdet[c,1] c = c + 1 } if (trend == "ct" | trend == "ctt") { uhat = uhat - bdet[c,1] * (1::T) / T c = c + 1 } if (trend == "ctt") { uhat = uhat - bdet[c,1] * ((1::T):^2) / T^2 } } RSS = uhat' * uhat ybar = mean(Y) TSS = (Y :- ybar)' * (Y :- ybar) R2 = 1 - RSS/TSS st_matrix(bname, bout) st_matrix(Vname, Vout) st_numscalar(rssname, RSS) st_numscalar(r2name, R2) st_numscalar(Nname, T) if (resvar != "") { touseflag = st_data(., touse) full = J(st_nobs(), 1, .) ki = 1 for (ti=1; ti<=st_nobs(); ti++) { if (touseflag[ti] == 1) { full[ti] = uhat[ki] ki = ki + 1 } } (void) st_addvar("double", resvar) st_store(., resvar, full) } } // ---------------------------------------------------------------------- // TAOLS adaptive F-test (Sun et al 2026) // ---------------------------------------------------------------------- real colvector __mc_diff_demean(real colvector x) { real colvector dx real scalar n n = rows(x) if (n < 2) return(J(1,1,0)) dx = x[2..n] - x[1..(n-1)] return(dx :- mean(dx)) } void _mc_taols_test_compute(string scalar Ycumvar, string scalar yflowvar, string scalar Xcums, string scalar Xflows, string scalar Dxs, string scalar trend, real scalar K, real scalar nx, string scalar touse, string scalar Fm_n, string scalar Fmp_n, string scalar Fc_n, string scalar Fcp_n, string scalar Fa_n, string scalar Fap_n, string scalar w_n) { real matrix Y, yf, Xc, Xf, Dx real matrix Phi, Vyc, Vy, Vxc, Vxf, Vdx, Vl, Vtr real matrix Xm, b_m, e_m, XXm, sigma_m, V_m real matrix Xc2, b_c, e_c, XXc, sigma_c, V_c real matrix gam_m, S_m, Wm, gam_c, S_c, Wc real scalar T, i, t, dfm, dfc, Omega_zz, Sigma_zz, ratio real scalar kappa, weight, Wa, dfden_m, dfden_c, p_m, p_c, p_a Y = st_data(., Ycumvar, touse) yf = st_data(., yflowvar, touse) Xc = st_data(., tokens(Xcums), touse) Xf = st_data(., tokens(Xflows), touse) Dx = st_data(., tokens(Dxs), touse) T = rows(Y) Phi = J(T, K, 0) for (i=1; i<=K; i++) { for (t=1; t<=T; t++) { Phi[t,i] = sqrt(2)*sin( (i - 0.5)*pi()*(t/T) ) } } Vyc = (Phi'*Y) / sqrt(T) Vy = (Phi'*yf) / sqrt(T) Vxc = (Phi'*Xc) / sqrt(T) Vxf = (Phi'*Xf) / sqrt(T) Vdx = (Phi'*Dx) / sqrt(T) Vl = colsum(Phi)'/sqrt(T) Vtr = J(K, 0, 0) if (trend == "ct" | trend == "ctt") Vtr = (Vtr, (Phi'*(1::T))/sqrt(T)/T) if (trend == "ctt") Vtr = (Vtr, (Phi'*((1::T):^2))/sqrt(T)/T^2) Xm = Vxc, Vxf, Vdx if (trend != "none") Xm = Xm, Vl if (trend == "ct" | trend == "ctt") Xm = Xm, Vtr[,1] if (trend == "ctt") Xm = Xm, Vtr[,2] XXm = quadcross(Xm, Xm) b_m = qrsolve(XXm, quadcross(Xm, Vyc)) e_m = Vyc - Xm*b_m dfm = K - cols(Xm) if (dfm <= 0) dfm = 1 sigma_m = (e_m'*e_m)/dfm V_m = sigma_m * invsym(XXm) Xc2 = Vxf, Vdx if (trend != "none") Xc2 = Xc2, Vl if (trend == "ct" | trend == "ctt") Xc2 = Xc2, Vtr[,1] if (trend == "ctt") Xc2 = Xc2, Vtr[,2] XXc = quadcross(Xc2, Xc2) b_c = qrsolve(XXc, quadcross(Xc2, Vy)) e_c = Vy - Xc2*b_c dfc = K - cols(Xc2) if (dfc <= 0) dfc = 1 sigma_c = (e_c'*e_c)/dfc V_c = sigma_c * invsym(XXc) gam_m = b_m[(nx+1)..(2*nx),1] S_m = V_m[(nx+1)..(2*nx),(nx+1)..(2*nx)] Wm = (gam_m'*invsym(S_m)*gam_m)/nx gam_c = b_c[(nx+1)..(2*nx),1] S_c = V_c[(nx+1)..(2*nx),(nx+1)..(2*nx)] Wc = (gam_c'*invsym(S_c)*gam_c)/nx Omega_zz = sum(e_m:^2)/dfm Sigma_zz = sum(__mc_diff_demean(yf):^2)/T if (Sigma_zz <= 0) Sigma_zz = 1 ratio = Omega_zz/Sigma_zz kappa = 0.5 weight = exp(-(T^kappa) * ratio) if (weight < 0) weight = 0 if (weight > 1) weight = 1 Wa = weight*Wm + (1-weight)*Wc dfden_m = K - 3*nx - 1 dfden_c = K - 2*nx if (dfden_m < 1) dfden_m = 1 if (dfden_c < 1) dfden_c = 1 p_m = 1 - F(nx, dfden_m, Wm) p_c = 1 - F(nx, dfden_c, Wc) p_a = 1 - F(nx, dfden_m, Wa) if (p_m == .) p_m = 1 if (p_c == .) p_c = 1 if (p_a == .) p_a = 1 st_numscalar(Fm_n, Wm) st_numscalar(Fmp_n, p_m) st_numscalar(Fc_n, Wc) st_numscalar(Fcp_n, p_c) st_numscalar(Fa_n, Wa) st_numscalar(Fap_n, p_a) st_numscalar(w_n, weight) } end