*! _xtmvu_mata v1.0.6 11oct2026 Yusuf Toyin Yusuf *! Mata engine: panel multivariate ARDL unit root / cointegration test *! (extends Sam, McNown, Goh and Goh 2025). Loaded with -do- by _xtmvu_load. mata: mata set matastrict off // version stamp: the loader reloads the engine when an older copy is still in Mata memory real scalar xmu_loaded() { return(106) } // columns of the F-test (1-based, in the design matrix) real rowvector xmu_fidx(real scalar k, real scalar cs) { real rowvector f f = (2..k+1) if (cs==4) f = f, (k+3) return(f) } real colvector xmu_dy(real colvector y) { real scalar T real colvector dy T = rows(y) dy = J(T,1,.) dy[|2\T|] = y[|2\T|] - y[|1\T-1|] return(dy) } // design matrix: [Ly, Lx(k), const, (trend), Dy lags(p), DX lags(q*k), DX contemporaneous(k)] // dependent rows are s0+1..T (s0 = number of pre-sample rows) real matrix xmu_Z(real colvector y, real matrix X, real scalar p, real scalar q, real scalar cs, real scalar s0) { real scalar T, k, n, j real colvector idx, dy real matrix Z, dX T = rows(y); k = cols(X); n = T - s0 idx = (s0+1::T) dy = xmu_dy(y) dX = J(T,k,.) dX[|2,1\T,k|] = X[|2,1\T,k|] - X[|1,1\T-1,k|] Z = y[idx :- 1], X[idx :- 1, .], J(n,1,1) if (cs==4 | cs==5) Z = Z, idx for (j=1; j<=p; j++) Z = Z, dy[idx :- j] for (j=1; j<=q; j++) Z = Z, dX[idx :- j, .] Z = Z, dX[idx, .] return(Z) } // Breusch-Godfrey LM test (given order) p-value for residuals r of a regression on Z real scalar xmu_bgp(real matrix Z, real colvector r, real scalar order) { real scalar n, o, R2 real matrix W real colvector e, b, res n = rows(r) W = Z for (o=1; o<=order; o++) { e = J(o,1,0) \ r[|1\n-o|] W = W, e } b = invsym(cross(W,W))*cross(W,r) res = r - W*b R2 = 1 - cross(res,res)/cross(r,r) return(chi2tail(order, n*R2)) } // Jarque-Bera normality test p-value real scalar xmu_jbp(real colvector r) { real scalar n, m2, S, K n = rows(r) m2 = mean(r:^2) S = mean(r:^3)/m2^1.5 K = mean(r:^4)/m2^2 return(chi2tail(2, n*(S^2/6 + (K-3)^2/24))) } // centred R-squared of y on a constant and X real scalar xmu_r2c(real colvector y, real matrix X) { real matrix W real colvector b, res, yc W = J(rows(y),1,1), X b = invsym(cross(W,W))*cross(W,y) res = y - W*b yc = y :- mean(y) return(1 - cross(res,res)/cross(yc,yc)) } // Ramsey RESET (squares and cubes of the fitted values) p-value real scalar xmu_resetp(real matrix Z, real colvector dep, real colvector r) { real scalar n, m, rss0, rss1, F real colvector fit, b2, r2 real matrix Z2 n = rows(Z) m = cols(Z) fit = dep - r Z2 = Z, fit:^2, fit:^3 b2 = invsym(cross(Z2,Z2))*cross(Z2,dep) r2 = dep - Z2*b2 rss0 = cross(r,r) rss1 = cross(r2,r2) F = ((rss0-rss1)/2)/(rss1/(n-m-2)) return(Ftail(2, n-m-2, F)) } // residual diagnostics of one unit: BG(1), BG(2), Jarque-Bera, Breusch-Pagan (Koenker), ARCH(1), RESET real rowvector xmu_diagrow(real matrix Z, real colvector dep, real colvector r) { real scalar n, bp, ar real colvector fit, u n = rows(r) fit = dep - r bp = chi2tail(1, n*xmu_r2c(r:^2, fit)) u = r:^2 ar = chi2tail(1, (n-1)*xmu_r2c(u[|2\n|], u[|1\n-1|])) return((xmu_bgp(Z,r,1), xmu_bgp(Z,r,2), xmu_jbp(r), bp, ar, xmu_resetp(Z,dep,r))) } // panel-level residual quantities of one unit: n, skewness, kurtosis, residual variance, df, // variance inflation factors of Ly and Lx(1..k), largest VIF among the stationary regressors, condition number real rowvector xmu_panelrow(real matrix Z, real colvector r, real scalar k, real scalar ntr) { real scalar n, m, m2, S, K, s2, vst, cond real rowvector idx, ssq, ev real colvector vif real matrix C, Zs n = rows(Z) m = cols(Z) m2 = mean(r:^2) S = mean(r:^3)/m2^1.5 K = mean(r:^4)/m2^2 s2 = cross(r,r)/(n-m) idx = (1..k+1), ((k+3)..m) C = correlation(Z[., idx]) vif = diagonal(invsym(C)) vst = max(vif[|k+2+ntr \ rows(vif)|]) ssq = sqrt(colsum(Z:^2)) Zs = Z :/ ssq ev = symeigenvalues(cross(Zs,Zs)) cond = sqrt(max(ev)/min(ev)) return((n, S, K, s2, n-m, vif[|1 \ k+1|]', vst, cond)) } // (t on Ly, F on fidx columns, rss, m) real rowvector xmu_stats(real matrix Z, real colvector dep, real rowvector fidx) { real scalar n, m, rss, s2, t, F real matrix A real colvector b, r, bf n = rows(Z); m = cols(Z) A = invsym(cross(Z,Z)) b = A*cross(Z,dep) r = dep - Z*b rss = cross(r,r) s2 = rss/(n-m) t = b[1]/sqrt(s2*A[1,1]) bf = b[fidx',1] F = (bf' * invsym(A[fidx,fidx]) * bf) / (s2*cols(fidx)) return((t, F, rss, m)) } // lag selection on a common sample (first maxlag+1 rows are pre-sample); crit 1=AIC 2=BIC real rowvector xmu_select(real colvector y, real matrix X, real scalar maxlag, real scalar cs, real scalar crit) { real scalar p, q, s0, n, m, ic, best, pen, bp, bq, rss, T real matrix Z, A real colvector dep, b, r, dy T = rows(y) s0 = maxlag + 1 n = T - s0 pen = (crit==1 ? 2 : ln(n)) dy = xmu_dy(y) dep = dy[|s0+1\T|] best = .; bp = .; bq = . for (p=0; p<=maxlag; p++) { for (q=0; q<=maxlag; q++) { Z = xmu_Z(y,X,p,q,cs,s0) m = cols(Z) if (m > n-8) continue A = invsym(cross(Z,Z)) b = A*cross(Z,dep) r = dep - Z*b rss = cross(r,r) ic = n*ln(rss/n) + pen*m if (best==. | ic < best) { best = ic; bp = p; bq = q } } } return((bp,bq)) } // one-unit bootstrap: joint null (pi=0, delta=0 [, trend=0 in Case 4]), drift-free generating // process, residuals centred/rescaled and resampled with the supplied index matrix ridx (B x n). // resel=1: lags are re-selected inside every draw (all (p,q) up to maxlag, criterion crit). real matrix xmu_bunit(real colvector y, real matrix X, real scalar p, real scalar q, real scalar cs, real scalar s0, real scalar B, real matrix ridx, real scalar resel, real scalar maxlag, real scalar crit) { real scalar T, k, n, m, ntr, b, j, l, t, v, pp_, qq, ic, best, pen, mm, bt, bF real rowvector keep, kidx, posd, posf, fidx, st real colvector dep, br, er, fp, phi, ys, dys, idx, dy real matrix Z, Zr, Zs, Zb, out T = rows(y); k = cols(X); ntr = (cs==4 | cs==5) Z = xmu_Z(y,X,p,q,cs,s0) n = rows(Z); m = cols(Z) dy = xmu_dy(y) dep = dy[|s0+1\T|] idx = (s0+1::T) fidx = xmu_fidx(k,cs) if (cs==4) keep = (k+2), (k+4..m) else keep = (k+2..m) Zr = Z[., keep] br = invsym(cross(Zr,Zr))*cross(Zr,dep) er = dep - Zr*br er = (er :- mean(er)) * sqrt(n/(n-cols(Zr))) kidx = J(1,m,0) kidx[keep] = (1..cols(keep)) posf = (k+3+ntr+p..m) fp = Z[., posf] * br[kidx[posf]',1] if (p>0) { posd = (k+3+ntr..k+2+ntr+p) phi = br[kidx[posd]',1] } else { posd = J(1,0,.) phi = J(0,1,.) } pen = (crit==1 ? 2 : ln(n)) ys = J(T,1,0); dys = J(T,1,0) ys[|1\s0|] = y[|1\s0|] if (s0>=2) dys[|2\s0|] = y[|2\s0|] - y[|1\s0-1|] out = J(B,2,.) for (b=1; b<=B; b++) { for (j=1; j<=n; j++) { t = s0 + j v = fp[j] + er[ridx[b,j]] for (l=1; l<=p; l++) v = v + phi[l]*dys[t-l] dys[t] = v ys[t] = ys[t-1] + v } if (resel==0) { Zs = Z Zs[.,1] = ys[idx :- 1] for (l=1; l<=p; l++) Zs[., posd[l]] = dys[idx :- l] st = xmu_stats(Zs, dys[idx], fidx) out[b,1] = st[1]; out[b,2] = st[2] } else { best = . for (pp_=0; pp_<=maxlag; pp_++) { for (qq=0; qq<=maxlag; qq++) { Zb = xmu_Z(ys,X,pp_,qq,cs,s0) mm = cols(Zb) if (mm > n-8) continue st = xmu_stats(Zb, dys[idx], fidx) ic = n*ln(st[3]/n) + pen*mm if (best==. | ic < best) { best = ic; bt = st[1]; bF = st[2] } } } out[b,1] = bt; out[b,2] = bF } } return(out) } real colvector xmu_cap(real colvector x, real scalar cap) { return(x :* (x :< cap) :+ cap :* (x :>= cap)) } // quantile of a bootstrap distribution (lower=1: lower-tail critical value, 0: upper-tail) real scalar xmu_cv(real colvector x, real scalar a, real scalar lower) { real colvector s real scalar B, i s = sort(x,1) B = rows(s) if (lower) i = max((1, ceil(a*B))) else i = min((B, max((1, ceil((1-a)*B))))) return(s[i]) } real scalar xmu_pval(real scalar obs, real colvector star, real scalar upper) { real scalar B B = rows(star) if (upper) return((1 + sum(star :>= obs))/(B+1)) return((1 + sum(star :<= obs))/(B+1)) } void xmu_panel_run(string scalar vars, string scalar pv, string scalar tv, string scalar touse, real scalar cs, real scalar maxlag, real scalar crit, real scalar fp_, real scalar fq_, real scalar B, real scalar csdmode, real scalar alpha) { external real colvector xmu_tbs, xmu_fbs string rowvector names real matrix data, info, Rc, U, ridx, Cd, tstar, Fstar, ptstar, pFstar, PAN, UNI, X, Z, bs real colvector id, tm, y, tt, dep, dy, kept, ok, ra, rb, r, b, tobs, Fobs, pt, pF, rk, rkF real colvector tbs_, fbs_, ftobs_s, intobs_s, tstat_s real rowvector fidx, pq, st real scalar N0, Nk, i, u, k, ntr, tmin, tmax, Tspan, s0, n, m, p, q, Ti, cstart, a, c, ia, ib, Tij, rho, CD, shared, balanced, nmax, cap, bad real scalar tbar_o, fish_t_o, invn_t_o, Fbar_o, fish_F_o, invn_F_o, nrow, minratio, maxm real colvector tbar_s, fish_t_s, invn_t_s, Fbar_s, fish_F_s, invn_F_s real rowvector cvs real matrix Rho st_numscalar("_xmu_rc", 0) names = tokens(vars) data = st_data(., names, touse) id = st_data(., pv, touse) tm = st_data(., tv, touse) k = cols(data) - 1 ntr = (cs==4 | cs==5) info = panelsetup(id, 1) N0 = rows(info) tmin = min(tm); tmax = max(tm); Tspan = tmax - tmin + 1 fidx = xmu_fidx(k, cs) U = J(N0, 16+k+8, 0) Rc = J(Tspan, N0, .) for (i=1; i<=N0; i++) { y = data[|info[i,1],1 \ info[i,2],1|] X = data[|info[i,1],2 \ info[i,2],k+1|] tt = tm[|info[i,1] \ info[i,2]|] Ti = rows(y) if (Ti > 1) { if (any((tt[|2\Ti|] - tt[|1\Ti-1|]) :!= 1)) { st_numscalar("_xmu_rc", 459) st_numscalar("_xmu_badunit", id[info[i,1]]) return } } if (crit > 0) s0 = maxlag + 1 else s0 = max((fp_, fq_)) + 1 n = Ti - s0 if (crit > 0) { if (n - (2*k+2+ntr) < 8) continue pq = xmu_select(y, X, maxlag, cs, crit) if (pq[1] == .) continue p = pq[1]; q = pq[2] } else { p = fp_; q = fq_ } Z = xmu_Z(y, X, p, q, cs, s0) m = cols(Z) if (n - m < 8) continue dy = xmu_dy(y) dep = dy[|s0+1\Ti|] st = xmu_stats(Z, dep, fidx) b = invsym(cross(Z,Z))*cross(Z,dep) r = dep - Z*b cstart = tt[1] - tmin + 1 Rc[|cstart+s0,i \ cstart+Ti-1,i|] = r U[i,.] = (1, p, q, n, m, st[1], st[2], cstart, s0, Ti, xmu_diagrow(Z, dep, r), xmu_panelrow(Z, r, k, ntr)) } kept = selectindex(U[.,1] :== 1) Nk = rows(kept) if (Nk < 2) { st_numscalar("_xmu_rc", 198) return } // Pesaran (2004) CD test on the ARDL residuals (pairwise, overlap >= 5) CD = 0 for (a=1; a= 5) { ra = Rc[ok,ia]; rb = Rc[ok,ib] Rho = correlation((ra,rb)) rho = Rho[1,2] if (rho != .) CD = CD + sqrt(Tij)*rho } } } CD = sqrt(2/(Nk*(Nk-1))) * CD if (csdmode == 1) shared = 1 else if (csdmode == 0) shared = 0 else shared = (abs(CD) > invnormal(0.975)) balanced = (min(U[kept,8]) == max(U[kept,8])) & (min(U[kept,4]) == max(U[kept,4])) nmax = max(U[kept,4]) if (shared) { if (balanced) Cd = floor(nmax*runiform(B,nmax)) :+ 1 else Cd = floor(Tspan*runiform(B,nmax)) :+ 1 } tstar = J(B, Nk, .) Fstar = J(B, Nk, .) for (u=1; u<=Nk; u++) { i = kept[u] y = data[|info[i,1],1 \ info[i,2],1|] X = data[|info[i,1],2 \ info[i,2],k+1|] p = U[i,2]; q = U[i,3]; n = U[i,4]; s0 = U[i,9]; cstart = U[i,8] if (shared) { if (balanced) ridx = Cd[., (1..n)] else { ridx = Cd[., (1..n)] :- (cstart + s0) :+ 1 bad = 0 ridx = ridx :* ((ridx :>= 1) :& (ridx :<= n)) :+ (floor(n*runiform(B,n)) :+ 1) :* (1 :- ((ridx :>= 1) :& (ridx :<= n))) } } else ridx = floor(n*runiform(B,n)) :+ 1 bs = xmu_bunit(y, X, p, q, cs, s0, B, ridx, (crit>0), maxlag, crit) tstar[.,u] = bs[.,1] Fstar[.,u] = bs[.,2] } tobs = U[kept,6] Fobs = U[kept,7] pt = ((1 :+ colsum(tstar :<= tobs'))/(B+1))' pF = ((1 :+ colsum(Fstar :>= Fobs'))/(B+1))' ptstar = J(B, Nk, .) pFstar = J(B, Nk, .) for (u=1; u<=Nk; u++) { rk = invorder(order(tstar[.,u],1)) ptstar[.,u] = rk/(B+1) rkF = invorder(order(Fstar[.,u],1)) pFstar[.,u] = (B :- rkF :+ 1)/(B+1) } cap = 1 - 0.5/(B+1) tbar_o = mean(tobs) fish_t_o = -2*sum(ln(pt)) invn_t_o = sum(invnormal(xmu_cap(pt,cap)))/sqrt(Nk) Fbar_o = mean(Fobs) fish_F_o = -2*sum(ln(pF)) invn_F_o = sum(invnormal(xmu_cap(pF,cap)))/sqrt(Nk) tbar_s = rowsum(tstar)/Nk fish_t_s = -2*rowsum(ln(ptstar)) invn_t_s = rowsum(invnormal(ptstar))/sqrt(Nk) Fbar_s = rowsum(Fstar)/Nk fish_F_s = -2*rowsum(ln(pFstar)) invn_F_s = rowsum(invnormal(pFstar))/sqrt(Nk) // panel matrix: stat, p-value, cv10, cv05, cv025, cv01, cv(alpha); rows tbar fish_t invn_t Fbar fish_F invn_F PAN = J(6,7,.) PAN[1,1] = tbar_o; PAN[1,2] = xmu_pval(tbar_o, tbar_s, 0) PAN[2,1] = fish_t_o; PAN[2,2] = xmu_pval(fish_t_o, fish_t_s, 1) PAN[3,1] = invn_t_o; PAN[3,2] = xmu_pval(invn_t_o, invn_t_s, 0) PAN[4,1] = Fbar_o; PAN[4,2] = xmu_pval(Fbar_o, Fbar_s, 1) PAN[5,1] = fish_F_o; PAN[5,2] = xmu_pval(fish_F_o, fish_F_s, 1) PAN[6,1] = invn_F_o; PAN[6,2] = xmu_pval(invn_F_o, invn_F_s, 0) cvs = (0.10, 0.05, 0.025, 0.01, alpha) for (c=1; c<=5; c++) { PAN[1,2+c] = xmu_cv(tbar_s, cvs[c], 1) PAN[2,2+c] = xmu_cv(fish_t_s, cvs[c], 0) PAN[3,2+c] = xmu_cv(invn_t_s, cvs[c], 1) PAN[4,2+c] = xmu_cv(Fbar_s, cvs[c], 0) PAN[5,2+c] = xmu_cv(fish_F_s, cvs[c], 0) PAN[6,2+c] = xmu_cv(invn_F_s, cvs[c], 1) } UNI = J(Nk, 10, .) for (u=1; u<=Nk; u++) { i = kept[u] UNI[u,1] = id[info[i,1]] UNI[u,2] = U[i,10] UNI[u,3] = U[i,2] UNI[u,4] = U[i,3] UNI[u,5] = tobs[u] UNI[u,6] = Fobs[u] UNI[u,7] = pt[u] UNI[u,8] = pF[u] UNI[u,9] = 1 + (pt[u] < alpha) + 2*(pF[u] < alpha) UNI[u,10] = U[i,11] } // diagnostics: obs per parameter (selected model) and for the largest model searched minratio = min(U[kept,4] :/ U[kept,5]) maxm = 1 + k + 1 + ntr + maxlag + maxlag*k + k DIAG = U[kept,1], U[kept,11..16] for (u=1; u<=Nk; u++) DIAG[u,1] = id[info[kept[u],1]] st_matrix("_xmu_diag", DIAG) PDIAG = U[kept,1], U[kept,17..16+k+8] for (u=1; u<=Nk; u++) PDIAG[u,1] = id[info[kept[u],1]] st_matrix("_xmu_pdiag", PDIAG) st_numscalar("_xmu_nvif", sum(U[kept,23+k] :> 10)) st_numscalar("_xmu_maxvif", max(U[kept,23+k])) st_matrix("_xmu_unit", UNI) st_matrix("_xmu_panel", PAN) st_numscalar("_xmu_N", Nk) st_numscalar("_xmu_Ndrop", N0 - Nk) st_numscalar("_xmu_CD", CD) st_numscalar("_xmu_CDp", 2*(1 - normal(abs(CD)))) st_numscalar("_xmu_shared", shared) st_numscalar("_xmu_balanced", balanced) st_numscalar("_xmu_Tmin", min(U[kept,10])) st_numscalar("_xmu_Tmax", max(U[kept,10])) st_numscalar("_xmu_Tavg", mean(U[kept,10])) st_numscalar("_xmu_minratio", minratio) st_numscalar("_xmu_searchratio", min(U[kept,4])/maxm) st_numscalar("_xmu_bgn", sum(U[kept,11] :< 0.05)) st_numscalar("_xmu_nlow", sum((U[kept,4]/maxm) :< 3)) st_numscalar("_xmu_maxm", maxm) st_numscalar("_xmu_k", k) xmu_tbs = tbar_s xmu_fbs = Fbar_s } // data for the graph command void xmu_to_data(string scalar umat) { external real colvector xmu_tbs, xmu_fbs real matrix U real scalar B, Nk, nobs real rowvector ix U = st_matrix(umat) B = rows(xmu_tbs) Nk = rows(U) nobs = max((B, Nk)) st_addobs(nobs) ix = st_addvar("double", ("tbs","fbs","uidx","ut","uF","upt","upF")) st_store((1::B), ix[1], xmu_tbs) st_store((1::B), ix[2], xmu_fbs) st_store((1::Nk), ix[3], (1::Nk)) st_store((1::Nk), ix[4], U[.,5]) st_store((1::Nk), ix[5], U[.,6]) st_store((1::Nk), ix[6], U[.,7]) st_store((1::Nk), ix[7], U[.,8]) } end