*! _xtmcg_mata 1.0.0  22may2026  Dr Merwan Roudane (merwanroudane920@gmail.com)
*! Mata core for the xtmulticointgrat package.
*! Implements: per-i first-level OLS, cumulation, ADF lag selection (t-sig/AIC/BIC),
*! Pedroni between-dimension Z_rho and Z_t statistics, PANIC factor extraction
*! (PC on first differences), Bai-Ng (2004) IC for factor selection, MQ_c test
*! for the number of stochastic trends in the factors, and pooled idiosyncratic ADF.
*! Loaded on demand via:
*!     qui capture mata: __xtmcg_loaded()
*!     if _rc qui _xtmcg_mata

version 14.0
mata
mata set matastrict off

void __xtmcg_loaded() { /* sentinel */ }


// ===========================================================================
// I.  Berenguer-Rico & Carrion-i-Silvestre (2006) MOMENTS  (Tables 1, 2, 3)
// ===========================================================================
// Tables 1 & 2 give finite-sample moments of  N^{-1/2} Z_rho_NT and Z_t_NT
// for four deterministic specs x five m1 (number of I(1) regressors, 0..4) x
// two m2 (number of I(2) regressors, 1..2) x four sample sizes T (50,100,250,1000).
//
// Returns a (5 x 16) matrix laid out as
//      [row m1=0..4 ;  col block (T=50,100,250,1000) x (Theta, Psi)]
// for the chosen (det, stat, m2).

real matrix _xtmcg_mom_block(string scalar det, string scalar stat, real scalar m2)
{
    real matrix M
    M = J(5,16,.)

    // Each row holds 4 (Theta,Psi) pairs = 8 columns; cols 9..16 unused.
    // Layout: [T50 Theta, T50 Psi, T100 Theta, T100 Psi, T250 Theta, T250 Psi,
    //          T1000 Theta, T1000 Psi]

    // ---------- NON-DETERMINISTIC ----------
    if (det=="none" & stat=="rho" & m2==1) {
        M[1,1..8] = (-5.654,26.739, -5.373,25.178, -5.213,24.894, -5.093,24.789)
        M[2,1..8] = (-10.037,45.394, -9.806,43.489, -9.620,42.471, -9.456,41.596)
        M[3,1..8] = (-14.189,65.113, -13.946,59.888, -13.849,59.566, -13.624,57.703)
        M[4,1..8] = (-18.157,85.567, -17.968,76.561, -17.899,74.861, -17.684,73.107)
        M[5,1..8] = (-21.994,108.195, -21.899,93.358, -21.898,90.320, -21.706,88.715)
        return(M)
    }
    if (det=="none" & stat=="rho" & m2==2) {
        M[1,1..8] = (-10.290,49.626, -9.616,43.187, -9.185,41.569, -8.823,40.649)
        M[2,1..8] = (-14.626,71.889, -13.989,61.523, -13.603,58.983, -13.214,56.987)
        M[3,1..8] = (-18.761,94.887, -18.142,79.496, -17.796,75.546, -17.414,72.969)
        M[4,1..8] = (-22.745,120.241, -22.166,97.678, -21.847,91.441, -21.473,88.537)
        M[5,1..8] = (-26.596,147.699, -26.094,116.281, -25.887,107.664, -25.507,103.787)
        return(M)
    }
    if (det=="none" & stat=="t" & m2==1) {
        M[1,1..8] = (-1.492,0.925, -1.377,0.875, -1.302,0.870, -1.250,0.908)
        M[2,1..8] = (-2.168,0.860, -2.059,0.778, -1.984,0.727, -1.935,0.706)
        M[3,1..8] = (-2.673,0.844, -2.554,0.735, -2.479,0.676, -2.421,0.630)
        M[4,1..8] = (-3.089,0.840, -2.963,0.730, -2.879,0.650, -2.816,0.594)
        M[5,1..8] = (-3.452,0.849, -3.319,0.725, -3.227,0.638, -3.160,0.577)
        return(M)
    }
    if (det=="none" & stat=="t" & m2==2) {
        M[1,1..8] = (-2.248,0.860, -2.070,0.750, -1.942,0.720, -1.846,0.743)
        M[2,1..8] = (-2.748,0.866, -2.579,0.740, -2.462,0.675, -2.376,0.649)
        M[3,1..8] = (-3.164,0.860, -2.995,0.729, -2.875,0.656, -2.790,0.606)
        M[4,1..8] = (-3.524,0.859, -3.353,0.731, -3.227,0.645, -3.140,0.584)
        M[5,1..8] = (-3.841,0.858, -3.671,0.726, -3.547,0.640, -3.454,0.571)
        return(M)
    }

    // ---------- CONSTANT ----------
    if (det=="c" & stat=="rho" & m2==1) {
        M[1,1..8] = (-10.381,44.055, -9.819,38.759, -9.453,36.351, -9.276,35.796)
        M[2,1..8] = (-14.259,64.275, -13.707,56.068, -13.329,52.227, -13.126,51.462)
        M[3,1..8] = (-18.144,85.923, -17.582,72.963, -17.251,68.480, -17.029,67.094)
        M[4,1..8] = (-22.060,110.768, -21.496,90.987, -21.185,84.285, -20.953,82.385)
        M[5,1..8] = (-25.878,136.910, -25.359,109.777, -25.125,100.824, -24.824,97.521)
        return(M)
    }
    if (det=="c" & stat=="rho" & m2==2) {
        M[1,1..8] = (-15.617,76.183, -14.387,59.183, -13.613,53.168, -13.240,51.578)
        M[2,1..8] = (-19.499,100.986, -18.286,77.465, -17.550,69.541, -17.101,66.796)
        M[3,1..8] = (-23.363,126.880, -22.163,95.559, -21.462,86.355, -20.987,82.306)
        M[4,1..8] = (-27.251,155.956, -26.081,116.434, -25.393,102.557, -24.895,97.005)
        M[5,1..8] = (-31.048,188.191, -29.933,137.205, -29.326,119.582, -28.764,111.988)
        return(M)
    }
    if (det=="c" & stat=="t" & m2==1) {
        M[1,1..8] = (-2.363,0.804, -2.218,0.723, -2.124,0.675, -2.071,0.663)
        M[2,1..8] = (-2.780,0.816, -2.641,0.722, -2.544,0.650, -2.489,0.617)
        M[3,1..8] = (-3.155,0.825, -3.015,0.720, -2.915,0.637, -2.855,0.595)
        M[4,1..8] = (-3.498,0.829, -3.356,0.722, -3.250,0.633, -3.184,0.577)
        M[5,1..8] = (-3.806,0.833, -3.664,0.724, -3.556,0.634, -3.480,0.568)
        return(M)
    }
    if (det=="c" & stat=="t" & m2==2) {
        M[1,1..8] = (-2.955,0.821, -2.742,0.702, -2.592,0.631, -2.502,0.608)
        M[2,1..8] = (-3.298,0.824, -3.100,0.713, -2.955,0.629, -2.863,0.587)
        M[3,1..8] = (-3.615,0.825, -3.426,0.713, -3.282,0.631, -3.187,0.578)
        M[4,1..8] = (-3.908,0.823, -3.728,0.715, -3.584,0.631, -3.486,0.565)
        M[5,1..8] = (-4.175,0.824, -4.007,0.718, -3.864,0.634, -3.759,0.559)
        return(M)
    }

    // ---------- LINEAR TREND ----------
    if (det=="ct" & stat=="rho" & m2==1) {
        M[1,1..8] = (-15.978,75.296, -14.712,58.648, -14.118,53.217, -13.706,51.180)
        M[2,1..8] = (-19.677,98.989, -18.407,76.561, -17.815,69.352, -17.423,67.149)
        M[3,1..8] = (-23.383,124.728, -22.166,94.558, -21.630,86.156, -21.222,82.757)
        M[4,1..8] = (-27.090,151.184, -25.947,113.290, -25.483,101.903, -25.062,97.387)
        M[5,1..8] = (-30.873,181.851, -29.763,133.997, -29.345,118.118, -28.907,112.825)
        return(M)
    }
    if (det=="ct" & stat=="rho" & m2==2) {
        M[1,1..8] = (-21.468,120.172, -19.259,83.279, -18.155,71.121, -17.472,66.520)
        M[2,1..8] = (-25.244,148.714, -23.055,103.193, -22.003,87.481, -21.281,82.437)
        M[3,1..8] = (-29.028,180.650, -26.858,122.963, -25.881,104.295, -25.115,97.409)
        M[4,1..8] = (-32.800,214.248, -30.689,144.101, -29.756,120.904, -28.988,111.843)
        M[5,1..8] = (-36.592,250.352, -34.532,167.071, -33.635,138.068, -32.839,127.237)
        return(M)
    }
    if (det=="ct" & stat=="t" & m2==1) {
        M[1,1..8] = (-2.972,0.791, -2.759,0.674, -2.634,0.609, -2.548,0.573)
        M[2,1..8] = (-3.300,0.804, -3.099,0.694, -2.974,0.619, -2.889,0.573)
        M[3,1..8] = (-3.604,0.814, -3.416,0.701, -3.291,0.627, -3.204,0.570)
        M[4,1..8] = (-3.889,0.816, -3.711,0.708, -3.588,0.628, -3.497,0.562)
        M[5,1..8] = (-4.158,0.823, -3.989,0.715, -3.864,0.629, -3.767,0.559)
        return(M)
    }
    if (det=="ct" & stat=="t" & m2==2) {
        M[1,1..8] = (-3.471,0.827, -3.195,0.700, -3.019,0.616, -2.895,0.568)
        M[2,1..8] = (-3.754,0.821, -3.501,0.708, -3.332,0.621, -3.210,0.567)
        M[3,1..8] = (-4.019,0.813, -3.786,0.708, -3.626,0.626, -3.501,0.560)
        M[4,1..8] = (-4.269,0.809, -4.055,0.709, -3.899,0.628, -3.774,0.553)
        M[5,1..8] = (-4.507,0.811, -4.309,0.711, -4.157,0.631, -4.028,0.552)
        return(M)
    }

    // ---------- QUADRATIC TREND ----------
    if (det=="ctt" & stat=="rho" & m2==1) {
        M[1,1..8] = (-21.887,117.198, -19.737,83.195, -18.556,69.827, -17.895,65.924)
        M[2,1..8] = (-25.485,145.813, -23.395,102.971, -22.239,86.462, -21.590,81.553)
        M[3,1..8] = (-29.169,177.125, -27.123,123.170, -25.992,103.017, -25.336,96.851)
        M[4,1..8] = (-32.826,210.241, -30.888,144.739, -29.814,119.555, -29.145,112.164)
        M[5,1..8] = (-36.543,248.073, -34.677,168.437, -33.646,136.887, -32.971,126.942)
        return(M)
    }
    if (det=="ctt" & stat=="rho" & m2==2) {
        M[1,1..8] = (-27.853,182.413, -24.477,115.187, -22.642,88.917, -21.549,80.452)
        M[2,1..8] = (-31.606,217.124, -28.277,137.932, -26.448,106.112, -25.343,96.868)
        M[3,1..8] = (-35.452,256.907, -32.103,160.992, -30.062,121.549, -28.911,109.766)
        M[4,1..8] = (-39.185,296.390, -35.962,185.851, -34.141,140.566, -32.996,126.828)
        M[5,1..8] = (-42.988,341.124, -39.786,211.252, -38.004,158.116, -36.848,141.614)
        return(M)
    }
    if (det=="ctt" & stat=="t" & m2==1) {
        M[1,1..8] = (-3.511,0.776, -3.238,0.669, -3.052,0.585, -2.934,0.545)
        M[2,1..8] = (-3.774,0.785, -3.529,0.686, -3.351,0.600, -3.236,0.549)
        M[3,1..8] = (-4.029,0.788, -3.806,0.692, -3.633,0.611, -3.518,0.549)
        M[4,1..8] = (-4.271,0.797, -4.068,0.699, -3.903,0.615, -3.785,0.549)
        M[5,1..8] = (-4.501,0.804, -4.317,0.702, -4.156,0.619, -4.037,0.547)
        return(M)
    }
    if (det=="ctt" & stat=="t" & m2==2) {
        M[1,1..8] = (-3.951,0.816, -3.630,0.709, -3.396,0.604, -3.233,0.550)
        M[2,1..8] = (-4.186,0.806, -3.898,0.709, -3.678,0.613, -3.519,0.550)
        M[3,1..8] = (-4.412,0.796, -4.153,0.703, -3.942,0.618, -3.786,0.548)
        M[4,1..8] = (-4.625,0.792, -4.396,0.704, -4.195,0.620, -4.038,0.546)
        M[5,1..8] = (-4.828,0.793, -4.622,0.700, -4.434,0.622, -4.278,0.545)
        return(M)
    }

    return(M)
}


// Linearly interpolate (Theta, Psi) for a given T from Table 1 / Table 2.
real rowvector _xtmcg_mom_lookup(string scalar det, string scalar stat,
                                  real scalar m1, real scalar m2, real scalar T)
{
    real matrix M
    real rowvector row, Tgrid, out
    real scalar mm1, mm2, T0, w, lo, hi

    mm1 = m1
    mm2 = m2
    if (mm1 < 0) mm1 = 0
    if (mm1 > 4) mm1 = 4
    if (mm2 < 1) mm2 = 1
    if (mm2 > 2) mm2 = 2

    M = _xtmcg_mom_block(det, stat, mm2)
    row = M[mm1+1, 1..8]
    Tgrid = (50, 100, 250, 1000)

    out = J(1, 2, 0)
    if (T <= Tgrid[1]) {
        out[1] = row[1]
        out[2] = row[2]
        return(out)
    }
    if (T >= Tgrid[4]) {
        out[1] = row[7]
        out[2] = row[8]
        return(out)
    }
    if (T <= Tgrid[2]) {
        lo = 1
        hi = 2
    } else if (T <= Tgrid[3]) {
        lo = 2
        hi = 3
    } else {
        lo = 3
        hi = 4
    }
    w = (T - Tgrid[lo])/(Tgrid[hi] - Tgrid[lo])
    out[1] = (1-w)*row[2*lo-1] + w*row[2*hi-1]
    out[2] = (1-w)*row[2*lo]   + w*row[2*hi]
    return(out)
}


// Table 3: idiosyncratic ADF moments (constant vs. trend specs)
//   "c"  -> ADF_e_c   "t" -> ADF_e_tau
real rowvector _xtmcg_mom_idio(string scalar spec, real scalar T)
{
    real matrix M
    real rowvector row, Tgrid, out
    real scalar lo, hi, w

    M = (-0.401, 1.167, -0.410, 1.054, -0.420, 0.996, -0.421, 0.970 \
         -1.563, 0.415, -1.554, 0.378, -1.540, 0.357, -1.529, 0.339)
    if (spec == "c") row = M[1, .]
    else             row = M[2, .]
    Tgrid = (50, 100, 250, 1000)
    out = J(1,2,0)
    if (T <= Tgrid[1]) {
        out = (row[1], row[2])
        return(out)
    }
    if (T >= Tgrid[4]) {
        out = (row[7], row[8])
        return(out)
    }
    if (T <= Tgrid[2]) {
        lo = 1
        hi = 2
    } else if (T <= Tgrid[3]) {
        lo = 2
        hi = 3
    } else {
        lo = 3
        hi = 4
    }
    w = (T - Tgrid[lo])/(Tgrid[hi] - Tgrid[lo])
    out[1] = (1-w)*row[2*lo-1] + w*row[2*hi-1]
    out[2] = (1-w)*row[2*lo]   + w*row[2*hi]
    return(out)
}


// ===========================================================================
// II.  Helper utilities
// ===========================================================================

// OLS coefficients (no intercept added internally)
real colvector _xtmcg_ols(real matrix Y, real matrix X)
{
    return(qrsolve(quadcross(X,X), quadcross(X,Y)))
}

// Detrend a column vector x by regressing on (constant) or (constant, trend)
real colvector _xtmcg_detrend(real colvector x, string scalar det)
{
    real scalar T
    real matrix W
    real colvector b, r
    T = rows(x)
    if (det == "none") return(x)
    if (det == "c") W = J(T,1,1)
    else if (det == "ct") W = (J(T,1,1), (1::T))
    else                  W = (J(T,1,1), (1::T), (1::T):^2)
    b = qrsolve(quadcross(W,W), quadcross(W,x))
    r = x - W*b
    return(r)
}

// Build deterministic regressor matrix for given length T
real matrix _xtmcg_det(real scalar T, string scalar det)
{
    if (det == "none") return(J(T,0,0))
    if (det == "c")    return(J(T,1,1))
    if (det == "ct")   return((J(T,1,1), (1::T)))
    return((J(T,1,1), (1::T), (1::T):^2))
}

// Cumulative sum of a column
real colvector _xtmcg_cumsum(real colvector x)
{
    real colvector y
    real scalar i, n
    n = rows(x)
    y = J(n,1,0)
    y[1] = x[1]
    for (i=2; i<=n; i++) y[i] = y[i-1] + x[i]
    return(y)
}


// ===========================================================================
// III.  ADF regression with lag selection (t-sig, AIC, BIC, fixed)
// ===========================================================================
// Returns (1 x 4) rowvector: (rho, t_rho, sum_phi, p_used)
real rowvector _xtmcg_adf(real colvector y, real scalar pmax, string scalar lagsel,
                          string scalar det)
{
    real scalar T, p, j, k, dfres, sigma2, rho, t_rho, sum_phi, best_p, best_ic
    real scalar tlast, ic, RSS
    real colvector dy, lev, b, e
    real matrix X, Z, det_mat
    real rowvector out

    T = rows(y)
    dy = y[2..T] :- y[1..(T-1)]    // length T-1
    lev = y[1..(T-1)]              // y_{t-1}

    if (lagsel == "fixed") {
        best_p = pmax
    } else if (lagsel == "tsig") {
        // Ng & Perron (1995) sequential t-rule: start at pmax, drop one
        // lag at a time while |t| on highest lag < 1.96.
        best_p = pmax
        for (p=pmax; p>=1; p--) {
            X = lev[(p+1)..(T-1)]
            for (j=1; j<=p; j++) {
                X = (X, dy[(p+1-j)..(T-1-j)])
            }
            det_mat = _xtmcg_det(rows(X), det)
            if (cols(det_mat) > 0) X = (X, det_mat)
            Z = dy[(p+1)..(T-1)]
            b = qrsolve(quadcross(X,X), quadcross(X,Z))
            e = Z - X*b
            dfres = rows(X) - cols(X)
            if (dfres <= 0) continue
            sigma2 = (e'e)/dfres
            if (sigma2 <= 0) continue
            tlast = b[p+1] / sqrt(sigma2 * invsym(quadcross(X,X))[p+1,p+1])
            if (abs(tlast) >= 1.96) {
                best_p = p
                break
            }
            if (p == 1) best_p = 0
        }
    } else {
        // AIC / BIC / HQIC: try p=0..pmax, pick min IC
        best_p = 0
        best_ic = .
        for (p=0; p<=pmax; p++) {
            X = lev[(pmax+1)..(T-1)]
            for (j=1; j<=p; j++) {
                X = (X, dy[(pmax+1-j)..(T-1-j)])
            }
            det_mat = _xtmcg_det(rows(X), det)
            if (cols(det_mat) > 0) X = (X, det_mat)
            Z = dy[(pmax+1)..(T-1)]
            b = qrsolve(quadcross(X,X), quadcross(X,Z))
            e = Z - X*b
            k = cols(X)
            if (rows(X) <= k) continue
            RSS = e'e
            if (RSS <= 0) continue
            sigma2 = RSS/rows(X)
            if (lagsel == "aic") {
                ic = ln(sigma2) + 2*k/rows(X)
            } else if (lagsel == "bic") {
                ic = ln(sigma2) + ln(rows(X))*k/rows(X)
            } else if (lagsel == "hqic") {
                ic = ln(sigma2) + 2*k*ln(ln(rows(X)))/rows(X)
            } else {
                ic = .
            }
            if (ic < best_ic | best_ic == .) {
                best_ic = ic
                best_p = p
            }
        }
    }

    // Final ADF regression with best_p lags
    p = best_p
    if (p > T-3) p = T-3
    if (p < 0) p = 0
    X = lev[(p+1)..(T-1)]
    for (j=1; j<=p; j++) {
        X = (X, dy[(p+1-j)..(T-1-j)])
    }
    det_mat = _xtmcg_det(rows(X), det)
    if (cols(det_mat) > 0) X = (X, det_mat)
    Z = dy[(p+1)..(T-1)]
    b = qrsolve(quadcross(X,X), quadcross(X,Z))
    e = Z - X*b
    dfres = rows(X) - cols(X)
    sigma2 = (e'e)/dfres
    rho = b[1]
    t_rho = rho / sqrt(sigma2 * invsym(quadcross(X,X))[1,1])
    sum_phi = 0
    if (p > 0) {
        for (j=1; j<=p; j++) sum_phi = sum_phi + b[1+j]
    }
    out = (rho, t_rho, sum_phi, p)
    return(out)
}


// ===========================================================================
// IV.  Per-i first-level cointegration: u = y - det α - X β,  S = Σu
// ===========================================================================
// y, X: (T x 1) and (T x k) matrices for one individual
// det:  "none","c","ct","ctt"
// Returns S_t (T x 1) cumulated residual.
real colvector _xtmcg_firstlevel(real colvector y, real matrix X, string scalar det)
{
    real matrix W, R
    real colvector b, u
    real scalar T, dx, dd
    T = rows(y)
    dx = cols(X)
    R = _xtmcg_det(T, det)
    dd = cols(R)
    W = (R, X)
    if (cols(W) == 0) {
        u = y
    } else {
        b = qrsolve(quadcross(W,W), quadcross(W,y))
        u = y - W*b
    }
    return(_xtmcg_cumsum(u))
}


// ===========================================================================
// V.  PANIC factor extraction on first differences of u-hat
// ===========================================================================
// U: T x N matrix of stage-2 residuals (one column per individual)
// det: deterministic spec for the second-stage regression
// rmax: maximum number of factors
// icname: "ic1","ic2","ic3","bic3"  -> Bai-Ng (2002) panel IC for factor number
//
// Returns Fhat, Lambdahat, ehat (idiosyncratic, recovered by cumsum), r_selected.
// Output via st_matrix for use back in Stata.

void _xtmcg_panic(string scalar Uname, string scalar det, real scalar rmax,
                   string scalar icname,
                   string scalar Fname, string scalar Lname, string scalar Ename,
                   string scalar rsel_name)
{
    real matrix U, dU, F, Lam, e, Fhat, Lhat, ehat, ehat_lev
    real matrix V, evecs
    real colvector evals
    real scalar T, N, Tm, k, r, r_star, ic_min, ic_val, sigma2, NT, c, i
    real rowvector sums

    U = st_matrix(Uname)
    T = rows(U)
    N = cols(U)
    if (rmax > min((T,N))-1) rmax = min((T,N))-1
    if (rmax < 1) rmax = 1

    // First-difference each column
    dU = U[2..T,.] - U[1..(T-1),.]
    Tm = rows(dU)

    // Demean columns (PANIC routine works on demeaned ΔU)
    for (i=1; i<=N; i++) dU[.,i] = dU[.,i] :- mean(dU[.,i])

    // PCA: eigen-decomposition of (1/(T*N)) ΔU * ΔU'
    // Use SVD of dU for numerical stability
    // dU = U * S * V', columns of U are principal components of rows
    V = dU * dU'        // Tm x Tm
    V = V / (Tm * N)
    symeigensystem(V, evecs, evals)   // returns DESCENDING (largest first)
    // Build factor estimates Fhat (Tm x k) and loadings Lhat (N x k) for each k
    // Select k via panel IC of Bai-Ng (2002).
    NT = Tm*N
    c  = min((Tm,N))

    r_star = 1
    ic_min = .
    for (k=1; k<=rmax; k++) {
        F = sqrt(Tm) * evecs[., 1..k]               // Tm x k
        Lam = (F' * dU)' / Tm                        // N x k
        e = dU - F * Lam'                            // Tm x N
        sigma2 = sum(e:^2) / NT
        if (icname == "ic1") ic_val = ln(sigma2) + k*((Tm+N)/NT)*ln(NT/(Tm+N))
        else if (icname == "ic2") ic_val = ln(sigma2) + k*((Tm+N)/NT)*ln(c)
        else if (icname == "ic3") ic_val = ln(sigma2) + k*(ln(c)/c)
        else if (icname == "bic3") ic_val = sigma2 + k*sigma2*((Tm+N-k)*ln(NT))/NT
        else ic_val = .
        if (ic_val < ic_min | ic_min == .) {
            ic_min = ic_val
            r_star = k
        }
    }

    // Final extraction at r_star
    r = r_star
    Fhat = sqrt(Tm) * evecs[., 1..r]
    Lhat = (Fhat' * dU)' / Tm
    ehat = dU - Fhat * Lhat'

    // Recover LEVELS of idiosyncratic + factor: cumulate by column
    ehat_lev = J(T, N, 0)
    for (i=1; i<=N; i++) ehat_lev[2..T, i] = _xtmcg_cumsum(ehat[., i])

    // Cumulate factors too
    real matrix Fhat_lev
    Fhat_lev = J(T, r, 0)
    for (i=1; i<=r; i++) Fhat_lev[2..T, i] = _xtmcg_cumsum(Fhat[., i])

    st_matrix(Fname, Fhat_lev)
    st_matrix(Lname, Lhat)
    st_matrix(Ename, ehat_lev)
    st_numscalar(rsel_name, r_star)
}


// ===========================================================================
// VI.  Pedroni between-dimension panel statistics for multicointegration
//      (cross-section independent case, Theorem 1).
// ===========================================================================
// Ymat: T x N panel of dep var (one column per i), Xmat: T x (N*m1) flow regressors,
// Xcum: T x (N*m2) cumulated I(2) regressors. m1, m2 fixed across i.
// det:  "none","c","ct","ctt"
// pmax, lagsel: ADF lag controls.
//
// Returns per-i ADF table and pooled stats via st_matrix / st_numscalar.

void _xtmcg_pedroni_indep(string scalar Yname, string scalar Xname, string scalar Xcumname,
                           real scalar m1, real scalar m2, string scalar det,
                           real scalar pmax, string scalar lagsel,
                           string scalar adfname, string scalar poolname)
{
    real matrix Y, X, Xc, R, W, ADFtab
    real colvector y_i, u_i, b
    real matrix X_i, Xc_i, Yfull
    real scalar T, N, i, p_used, rho_i, t_i, sumphi_i
    real scalar Z_rho_raw, Z_t_raw, Nused
    real rowvector adfres, pool
    real matrix Worig

    Y  = st_matrix(Yname)
    X  = st_matrix(Xname)
    Xc = st_matrix(Xcumname)
    T  = rows(Y)
    N  = cols(Y)

    ADFtab = J(N, 4, .)
    Z_rho_raw = 0
    Z_t_raw   = 0
    Nused     = 0

    for (i=1; i<=N; i++) {
        y_i = Y[.,i]
        if (m1 == 0) X_i = J(T,0,0)
        else         X_i = X[., (1+(i-1)*m1)..(i*m1)]
        if (m2 == 0) Xc_i = J(T,0,0)
        else         Xc_i = Xc[., (1+(i-1)*m2)..(i*m2)]

        // Build full regressor matrix: det, Xc_i (I(2)), X_i (I(1))
        R = _xtmcg_det(T, det)
        W = (R, Xc_i, X_i)
        if (cols(W) == 0) {
            u_i = y_i
        } else {
            b   = qrsolve(quadcross(W,W), quadcross(W,y_i))
            u_i = y_i - W*b
        }

        // ADF on u_i with NO deterministics (already partialled out)
        adfres = _xtmcg_adf(u_i, pmax, lagsel, "none")
        rho_i    = adfres[1]
        t_i      = adfres[2]
        sumphi_i = adfres[3]
        p_used   = adfres[4]
        ADFtab[i, 1] = rho_i
        ADFtab[i, 2] = t_i
        ADFtab[i, 3] = sumphi_i
        ADFtab[i, 4] = p_used

        if (rho_i != . & t_i != . & (1 - sumphi_i) != 0) {
            Z_rho_raw = Z_rho_raw + (T * rho_i)/(1 - sumphi_i)
            Z_t_raw   = Z_t_raw + t_i
            Nused = Nused + 1
        }
    }

    // Pool: N^{-1/2} Z = (1/sqrt(N)) * Z_raw
    real scalar Z_rho_pool, Z_t_pool, Theta1, Psi1, Theta2, Psi2, std_rho, std_t
    real rowvector mom_r, mom_t

    Z_rho_pool = Z_rho_raw / sqrt(Nused)
    Z_t_pool   = Z_t_raw   / sqrt(Nused)

    mom_r = _xtmcg_mom_lookup(det, "rho", m1, m2, T)
    mom_t = _xtmcg_mom_lookup(det, "t",   m1, m2, T)
    Theta1 = mom_r[1]; Psi1 = mom_r[2]
    Theta2 = mom_t[1]; Psi2 = mom_t[2]

    // Standardized statistics ~ N(0,1) under H0
    std_rho = (Z_rho_pool - Theta1*sqrt(Nused)) / sqrt(Psi1)
    std_t   = (Z_t_pool   - Theta2*sqrt(Nused)) / sqrt(Psi2)

    pool = (Z_rho_pool, Z_t_pool, std_rho, std_t, Theta1, Psi1, Theta2, Psi2, Nused, T)

    st_matrix(adfname, ADFtab)
    st_matrix(poolname, pool)
}


// ===========================================================================
// VII.  Common-factor branch: PANIC + pooled idiosyncratic ADF + MQ_c on F
// ===========================================================================
// For each i:  1) y_i on (det, x_i) -> ϑ̂_i, then S_i = cumsum(ϑ̂_i)
//              2) y_i on (det, S_i)  -> u_i  (Granger-Lee two-step, eq 11)
// Then assemble U (T x N), pass to PANIC factor extraction.
// Apply pooled ADF to idiosyncratic (LEVELS) using Table 3 moments.
// For factors: ADF on each (if r=1) or MQ_c test (if r>1).

void _xtmcg_twostep_resid(string scalar Yname, string scalar Xname, real scalar m1,
                           string scalar det,
                           string scalar Sname, string scalar Uname)
{
    real matrix Y, X, S, U
    real colvector y_i, S_i, b, u_i
    real matrix R, X_i, W
    real scalar T, N, i, dx

    Y = st_matrix(Yname)
    X = st_matrix(Xname)
    T = rows(Y)
    N = cols(Y)

    S = J(T, N, 0)
    U = J(T, N, 0)

    for (i=1; i<=N; i++) {
        y_i = Y[.,i]
        if (m1 == 0) X_i = J(T,0,0)
        else         X_i = X[., (1+(i-1)*m1)..(i*m1)]
        // Stage 1: y on (det, x) -> ϑ̂, S_t = cumsum
        S_i = _xtmcg_firstlevel(y_i, X_i, det)
        S[., i] = S_i
        // Stage 2: y on (det, S) -> u_t
        R = _xtmcg_det(T, det)
        W = (R, S_i)
        if (cols(W) == 0) {
            u_i = y_i
        } else {
            b   = qrsolve(quadcross(W,W), quadcross(W,y_i))
            u_i = y_i - W*b
        }
        U[., i] = u_i
    }

    st_matrix(Sname, S)
    st_matrix(Uname, U)
}


// Pooled idiosyncratic ADF: read ehat (T x N), run ADF on each col, pool.
void _xtmcg_pool_idio_adf(string scalar Ename, string scalar det,
                           real scalar pmax, string scalar lagsel,
                           string scalar adfname, string scalar poolname)
{
    real matrix E, ADFtab
    real colvector e_i
    real scalar T, N, i, Nused, sum_t, Theta_e, Psi_e, T_eff
    real rowvector adfres, pool, mom
    string scalar spec

    E = st_matrix(Ename)
    T = rows(E)
    N = cols(E)
    ADFtab = J(N, 4, .)
    sum_t = 0
    Nused = 0
    for (i=1; i<=N; i++) {
        e_i = E[.,i]
        adfres = _xtmcg_adf(e_i, pmax, lagsel, det)
        ADFtab[i,] = adfres
        if (adfres[2] != .) {
            sum_t = sum_t + adfres[2]
            Nused = Nused + 1
        }
    }
    T_eff = T
    spec = (det == "c" ? "c" : "t")
    mom = _xtmcg_mom_idio(spec, T_eff)
    Theta_e = mom[1]
    Psi_e   = mom[2]

    real scalar Zbar, Zstd
    Zbar = sum_t / sqrt(Nused)
    Zstd = (sum_t/Nused - Theta_e) * sqrt(Nused/Psi_e)

    pool = (Zbar, Zstd, Theta_e, Psi_e, Nused, T_eff)
    st_matrix(adfname, ADFtab)
    st_matrix(poolname, pool)
}


// MQ_c parametric (Bai-Ng 2004) - determine the number of stochastic trends
// among r estimated factors. We implement a simplified version: for each
// candidate number of common stochastic trends q (0..r), run an ADF on
// successively orthogonalised factors and apply the joint procedure.
//
// Returned: vector of t-statistics for each cumulative factor (1..r) after
// detrending; the user then reads the column showing how many > 5% c.v.

void _xtmcg_mq_test(string scalar Fname, string scalar det, real scalar pmax,
                     string scalar lagsel, string scalar outname)
{
    real matrix F, OUT
    real colvector f_j
    real scalar T, r, j
    real rowvector adfres

    F = st_matrix(Fname)
    T = rows(F)
    r = cols(F)
    OUT = J(r, 4, .)
    for (j=1; j<=r; j++) {
        f_j = F[., j]
        adfres = _xtmcg_adf(f_j, pmax, lagsel, det)
        OUT[j,] = adfres
    }
    st_matrix(outname, OUT)
}


// ===========================================================================
// VIII.  Build T x N matrices from a (long) Stata panel
// ===========================================================================
void _xtmcg_build_matrices(string scalar yvar, string scalar xvars,
                           string scalar ivar, string scalar tvar,
                           string scalar touse,
                           string scalar Yname, string scalar Xname,
                           string scalar Xcumname)
{
    real matrix data, Y, X, Xc
    real colvector ids, vx_col
    real scalar N, T, i, j, m1
    string rowvector xs
    real matrix sub

    xs = tokens(xvars)
    m1 = length(xs)

    data = st_data(., (ivar, tvar, yvar, xs), touse)
    ids  = uniqrows(data[.,1])
    N    = length(ids)
    T    = sum(data[.,1] :== ids[1])

    Y  = J(T, N, .)
    X  = J(T, N*m1, .)
    Xc = J(T, N*m1, .)

    for (i=1; i<=N; i++) {
        sub = select(data, data[.,1] :== ids[i])
        sub = sort(sub, 2)
        Y[., i] = sub[., 3]
        for (j=1; j<=m1; j++) {
            vx_col = sub[., 3+j]
            X[., (i-1)*m1 + j]  = vx_col
            Xc[., (i-1)*m1 + j] = _xtmcg_cumsum(vx_col)
        }
    }
    st_matrix(Yname,    Y)
    st_matrix(Xname,    X)
    st_matrix(Xcumname, Xc)
}


// ===========================================================================
// IX.  Persist factor / residual matrices back into the Stata dataset
// ===========================================================================
void _xtmcg_persist(string scalar Sname, string scalar Uname,
                    string scalar Ename, string scalar Fname,
                    string scalar ivar, string scalar tvar, string scalar touse)
{
    real matrix S, U, E, F
    real scalar T, N, r, i, k, t, nobs
    real colvector ids, idfull, tfull, tousefull
    real colvector slot_S, slot_U, slot_E
    real matrix slot_F
    real colvector mask, tvals_i, idx_full
    real matrix sortmat

    S = st_matrix(Sname)
    U = st_matrix(Uname)
    E = st_matrix(Ename)
    F = st_matrix(Fname)
    T = rows(S)
    N = cols(S)
    r = cols(F)

    nobs      = st_nobs()
    idfull    = st_data(., ivar)
    tfull     = st_data(., tvar)
    tousefull = st_data(., touse)
    ids       = uniqrows(select(idfull, tousefull))

    slot_S = J(nobs, 1, .)
    slot_U = J(nobs, 1, .)
    slot_E = J(nobs, 1, .)
    slot_F = J(nobs, r, .)

    for (i=1; i<=N; i++) {
        mask     = (idfull :== ids[i]) :& (tousefull :== 1)
        idx_full = select(range(1, nobs, 1), mask)
        tvals_i  = select(tfull, mask)
        sortmat  = sort((tvals_i, idx_full), 1)
        for (t=1; t<=rows(sortmat); t++) {
            slot_S[sortmat[t,2], 1] = S[t, i]
            slot_U[sortmat[t,2], 1] = U[t, i]
            slot_E[sortmat[t,2], 1] = E[t, i]
            for (k=1; k<=r; k++) slot_F[sortmat[t,2], k] = F[t, k]
        }
    }

    (void) st_addvar("double", "_xtmcg_S_i")
    (void) st_addvar("double", "_xtmcg_u_i")
    (void) st_addvar("double", "_xtmcg_e_i")
    st_store(., "_xtmcg_S_i", slot_S)
    st_store(., "_xtmcg_u_i", slot_U)
    st_store(., "_xtmcg_e_i", slot_E)
    for (k=1; k<=r; k++) {
        (void) st_addvar("double", "_xtmcg_F" + strofreal(k))
        st_store(., "_xtmcg_F" + strofreal(k), slot_F[., k])
    }
}

end