*! cupfm.ado - Panel Cointegration with Common Factors
*! Version 1.0.4 - 2026-04-20 (Fix: mata engine embedded at top level per Baum SSC guidelines)
*! Author: Dr. Merwan Roudane (merwanroudane920@gmail.com)
*!
*! Implements all estimators from:
*!   Bai, J., Kao, C. & Ng, S. (2009). Panel cointegration with global
*!     stochastic trends. Journal of Econometrics, 149(1), 82-99.
*!   Bai, J. & Kao, C. (2005). On the estimation and inference of a panel
*!     cointegration model with cross-sectional dependence.
*!     CPR Working Paper No. 75, Syracuse University. SSRN-1815227.
*!
*! Estimators (exactly as in GAUSS bkn source code):
*!   LSDV     : Within/fixed-effects (baseline, biased under CS dependence)
*!   Bai FM   : One-shot 2-step FM  [Bai & Kao 2005, Eq. 7-8]
*!   CupFM    : Continuously-Updated FM  [BKN 2009, Theorem 3, Eq. 16]
*!   CupFM-bar: CupFM with Z-bar instrument  [BKN 2009, alternative]
*!   CupBC    : Continuously-Updated BC  [BKN 2009, Theorem 2]
*!
*! Data: Stata long format (xtset panelvar timevar), balanced panel required
*! Long->Wide conversion handled internally in Mata (matches GAUSS bkn)
*!
*! Mata engine is embedded below at the top level: compiled once when
*! cupfm.ado is first sourced, persists until discard.

// ── MATA ENGINE (compiled on first load of cupfm.ado) ──────────────────────
capture mata: mata drop cupfm_*()
mata:

// ===================================================================
// 1. UTILITY: Stata long format → GAUSS wide "grouped-by-unit"
// ===================================================================
function cupfm_long2wide(V_long, N, T, m) {
    V_wide = J(T, N*m, .)
    for (i=1; i<=N; i++) {
        idx = (i-1)*T+1 :: i*T
        V_wide[., (i-1)*m+1 :: i*m] = V_long[idx, .]
    }
    return(V_wide)
}

// ===================================================================
// 2. UTILITY: Wide → Long
// ===================================================================
function cupfm_wide2long(V_wide, N, T, m) {
    V_long = J(N*T, m, .)
    for (i=1; i<=N; i++) {
        idx = (i-1)*T+1 :: i*T
        V_long[idx, .] = V_wide[., (i-1)*m+1 :: i*m]
    }
    return(V_long)
}

// ===================================================================
// 3. UTILITY: First differences
// ===================================================================
function cupfm_fdif(X) {
    T = rows(X)
    return(X[2::T, .] - X[1::T-1, .])
}

// ===================================================================
// 4. UTILITY: Within-demean each unit
// ===================================================================
function cupfm_demean(X_long, N, T) {
    Xdm = X_long :+ 0
    for (i=1; i<=N; i++) {
        idx = (i-1)*T+1 :: i*T
        Xdm[idx, .] = X_long[idx, .] :- mean(X_long[idx, .])
    }
    return(Xdm)
}

// ===================================================================
// 5. FACTOR EXTRACTION - pass-by-reference
// ===================================================================
void cupfm_factor(U_wide, r, F_out, L_out) {
    T = rows(U_wide)
    N = cols(U_wide)
    if (T >= N) {
        _P = J(T, T, .)
        _D = J(min((T,N)), 1, .)
        _Q = J(N, N, .)
        svd(U_wide, _P, _D, _Q)
        F_out = _P[., 1::r] :* T
        L_out = U_wide' * F_out :/ T^2
    }
    else {
        _P = J(N, N, .)
        _D = J(N, 1, .)
        _Q = J(N, N, .)
        svd(U_wide'*U_wide, _P, _D, _Q)
        L_out = _P[., 1::r] :* sqrt(N)
        F_out = U_wide * L_out :/ N
    }
}

// ===================================================================
// 6. ROTATION WEIGHTS a_ik  (BKN 2009 Section 2)
// ===================================================================
function cupfm_delta(Lambda) {
    N    = rows(Lambda)
    SigL = Lambda'*Lambda :/ N
    return(Lambda * invsym(SigL) * Lambda')
}

// ===================================================================
// 7. BAI-NG (2002) IC CRITERION FOR NUMBER OF FACTORS
// ===================================================================
function cupfm_ic_bn(U_wide, rmax) {
    T    = rows(U_wide)
    N    = cols(U_wide)
    NT   = N*T
    best_r  = 1
    best_ic = 1e15
    F_r = J(T, 1, .)
    L_r = J(N, 1, .)
    for (ri=1; ri<=rmax; ri++) {
        cupfm_factor(U_wide, ri, F_r, L_r)
        resid = U_wide - F_r * L_r'
        V_r   = sum(resid:^2) / NT
        if (V_r <= 0) continue
        ic_r = log(V_r) + ri * (N+T)/NT * log(NT/(N+T))
        if (ic_r < best_ic) {
            best_ic = ic_r
            best_r  = ri
        }
    }
    return(best_r)
}

// ===================================================================
// 8. BARTLETT KERNEL LONG-RUN COVARIANCE
// ===================================================================
function cupfm_kernel_unit(U, bw, Delta_plus) {
    T1    = rows(U)
    m     = cols(U)
    Sigma = U'*U :/ T1
    Gamma = J(m, m, 0)
    for (j=1; j<=bw; j++) {
        if (T1-j < 1) break
        wj      = 1 - j/(bw+1)
        Gamma_j = U[j+1::T1, .]' * U[1::T1-j, .] :/ T1
        Gamma   = Gamma + wj * Gamma_j
    }
    Delta_plus = Sigma + Gamma
    return(Sigma + Gamma + Gamma')
}

// ===================================================================
// 9. CUP PLAIN LS
// ===================================================================
function cupfm_plsbeta(X_long, y_long, F, N, T, k) {
    invFtF = invsym(F'*F)
    XX = J(k, k, 0)
    Xy = J(k, 1, 0)
    for (i=1; i<=N; i++) {
        idx = (i-1)*T+1 :: i*T
        xi  = X_long[idx, .]
        yi  = y_long[idx, .]
        FtX = F'*xi
        FtY = F'*yi
        XX  = XX + xi'*xi - FtX'*invFtF*FtX
        Xy  = Xy + xi'*yi - FtX'*invFtF*FtY
    }
    return(invsym(XX) * Xy)
}

// ===================================================================
// 10. FULL FM BETA - pass-by-reference for all outputs
// ===================================================================
void cupfm_fmbeta(y_wide_T1, X_long_T1, X_wide_T1, aik,
                  F_T1, dF, u1_wide_T1,
                  du2_T1_long, bw, N, T1, k, r,
                  beta_fm1, beta_fm2, tstat_fm1, tstat_fm2, Omega_bar, cond_var) {

    m = 1 + k + r

    du2_T1_wide = cupfm_long2wide(du2_T1_long, N, T1, k)
    du2N_wide   = J(T1, N*k, 0)
    for (ip=1; ip<=N; ip++) {
        sumd = J(T1, k, 0)
        for (kp=1; kp<=N; kp++) {
            sumd = sumd + du2_T1_wide[., (kp-1)*k+1::kp*k] :* aik[ip, kp]
        }
        du2N_wide[., (ip-1)*k+1::ip*k] = du2_T1_wide[., (ip-1)*k+1::ip*k] - sumd:/N
    }

    xbar_wide = J(T1, N*k, 0)
    for (ip=1; ip<=N; ip++) {
        sumx = J(T1, k, 0)
        for (kp=1; kp<=N; kp++) {
            sumx = sumx + X_wide_T1[., (kp-1)*k+1::kp*k] :* aik[ip, kp]
        }
        xbar_wide[., (ip-1)*k+1::ip*k] = X_wide_T1[., (ip-1)*k+1::ip*k] - sumx:/N
    }

    Omega_sum = J(m, m, 0)
    Dplus_sum = J(m, m, 0)
    Dp_k = J(m, m, 0)
    for (ku=1; ku<=N; ku++) {
        u1k   = u1_wide_T1[., ku]
        du2Nk = du2N_wide[., (ku-1)*k+1::ku*k]
        W_k   = u1k, du2Nk, dF
        Om_k  = cupfm_kernel_unit(W_k, bw, Dp_k)
        Omega_sum = Omega_sum + Om_k
        Dplus_sum = Dplus_sum + Dp_k
    }
    Omega_bar = Omega_sum :/ N
    Dplus_bar = Dplus_sum :/ N

    Omega_ub = Omega_bar[1, 2::m]
    Omega_b  = Omega_bar[2::m, 2::m]
    invOmb   = invsym(Omega_b)
    gab      = Omega_ub * invOmb

    y_plus = J(T1, N, .)
    for (ip=1; ip<=N; ip++) {
        dN_ip = du2N_wide[., (ip-1)*k+1::ip*k]
        for (ti=1; ti<=T1; ti++) {
            temp5          = (dN_ip[ti, .], dF[ti, .])'
            y_plus[ti, ip] = y_wide_T1[ti, ip] - gab*temp5
        }
    }

    abu  = Dplus_bar[2::m, 1]
    ab_  = Dplus_bar[2::m, 2::m]
    obu  = Omega_bar[2::m, 1]
    obup = abu - ab_*invOmb*obu

    invFF    = invsym(F_T1'*F_T1)
    db_all   = J(N, r*k, 0)
    for (ip=1; ip<=N; ip++) {
        x_ip = X_wide_T1[., (ip-1)*k+1::ip*k]
        db_i = invFF * (F_T1'*x_ip)
        db_all[ip, .] = rowshape(db_i, 1)
    }
    dbsum_all = J(N, r*k, 0)
    for (ip=1; ip<=N; ip++) {
        for (kp=1; kp<=N; kp++) {
            dbsum_all[ip, .] = dbsum_all[ip, .] + db_all[kp, .] :* aik[ip, kp]
        }
    }
    dbsum1_all = db_all - dbsum_all:/N
    db_bar     = rowshape(colsum(dbsum1_all), r)

    db2 = J(r, N*k, 0)
    for (ip=1; ip<=N; ip++) {
        db2[., (ip-1)*k+1::ip*k] = rowshape(dbsum1_all[ip, .], r)
    }

    obup_x = obup[1::k]
    obup_f = obup[k+1::k+r]
    temp3  = T1 :* (N*obup_x - db_bar'*obup_f)

    invFF1 = invsym(F_T1'*F_T1)
    XX = J(k, k, 0)
    Xy = J(k, 1, 0)
    for (ip=1; ip<=N; ip++) {
        xi  = X_long_T1[(ip-1)*T1+1::ip*T1, .]
        yi  = y_plus[., ip]
        FtX = F_T1'*xi
        FtY = F_T1'*yi
        XX  = XX + xi'*xi - FtX'*invFF1*FtX
        Xy  = Xy + xi'*yi - FtX'*invFF1*FtY
    }
    Xy     = Xy - temp3
    invXX  = invsym(XX)
    beta_fm1  = invXX * Xy

    Omega_uu = Omega_bar[1, 1]
    cond_var = Omega_uu - Omega_ub * invOmb * Omega_ub'
    tstat_fm1 = beta_fm1 :/ sqrt(diagonal(invXX) :* cond_var)

    z1     = xbar_wide - F_T1*db2
    z_long = cupfm_wide2long(z1, N, T1, k)

    XXz = J(k, k, 0)
    Xyz = J(k, 1, 0)
    for (ip=1; ip<=N; ip++) {
        zi  = z_long[(ip-1)*T1+1::ip*T1, .]
        yi  = y_plus[., ip]
        XXz = XXz + zi'*zi
        Xyz = Xyz + zi'*yi
    }
    Xyz      = Xyz - temp3
    invXXz   = invsym(XXz)
    beta_fm2 = invXXz * Xyz
    tstat_fm2 = beta_fm2 :/ sqrt(diagonal(invXXz) :* cond_var)
}

// ===================================================================
// 11. MAIN ENGINE - called from cupfm.ado
// ===================================================================
void cupfm_main(string scalar depvar,
                string scalar indepvars,
                string scalar touse,
                real scalar N,
                real scalar T,
                real scalar r_in,
                real scalar bw,
                string scalar ktype,
                real scalar maxiter,
                real scalar do_autoR,
                real scalar autormax,
                real scalar verbose) {

    real matrix y_long, X_long
    st_view(y_long, ., depvar, touse)
    st_view(X_long, ., tokens(indepvars), touse)
    k  = cols(X_long)
    NT = N * T
    r  = r_in

    Xdm_long  = cupfm_demean(X_long, N, T)
    ydm_long  = cupfm_demean(y_long, N, T)
    XX_lsdv   = Xdm_long'*Xdm_long
    beta_lsdv = invsym(XX_lsdv) * (Xdm_long'*ydm_long)

    uhat_lsdv   = ydm_long - Xdm_long*beta_lsdv
    sigma2_lsdv = (uhat_lsdv'*uhat_lsdv)[1,1] / (NT - N - k)
    if (sigma2_lsdv <= 0) sigma2_lsdv = 1e-12
    var_lsdv    = diagonal(sigma2_lsdv :* invsym(XX_lsdv))
    se_lsdv     = sqrt(var_lsdv)
    t_lsdv      = beta_lsdv :/ se_lsdv

    uhat_long = y_long - X_long*beta_lsdv
    U_wide    = cupfm_long2wide(uhat_long, N, T, 1)

    real matrix F_hat, L_hat
    F_hat = J(T, 1, .)
    L_hat = J(N, 1, .)

    if (do_autoR) {
        rmax_use = floor(min((autormax, min((N, T)) / 2)))
        if (rmax_use < 1) rmax_use = 1
        r = cupfm_ic_bn(U_wide, rmax_use)
    }

    cupfm_factor(U_wide, r, F_hat, L_hat)

    T1          = T - 1
    du2_T1_long = J(N*T1, k, .)
    for (i=1; i<=N; i++) {
        idx_i = (i-1)*T+1 :: i*T
        du2_T1_long[(i-1)*T1+1 :: i*T1, .] = cupfm_fdif(X_long[idx_i, .])
    }

    dF = cupfm_fdif(F_hat)

    FL_wide = J(T, N, 0)
    for (i=1; i<=N; i++) {
        FL_wide[., i] = F_hat * L_hat[i, .]'
    }
    FL_long = cupfm_wide2long(FL_wide, N, T, 1)
    u1_wide = cupfm_long2wide(y_long - X_long*beta_lsdv - FL_long, N, T, 1)

    aik = cupfm_delta(L_hat)

    F_T1       = F_hat[1::T1, .]
    u1_wide_T1 = u1_wide[1::T1, .]
    y_wide     = cupfm_long2wide(y_long, N, T, 1)
    y_wide_T1  = y_wide[1::T1, .]

    X_long_T1 = J(N*T1, k, .)
    for (i=1; i<=N; i++) {
        idx_i = (i-1)*T+1 :: i*T
        X_long_T1[(i-1)*T1+1 :: i*T1, .] = X_long[idx_i, .][1::T1, .]
    }
    X_wide_T1 = cupfm_long2wide(X_long_T1, N, T1, k)

    real matrix beta_bfm, beta_bfm2, ts_bfm, ts_bfm2, Om_bfm
    real scalar cv_bfm
    beta_bfm  = J(k, 1, .)
    beta_bfm2 = J(k, 1, .)
    ts_bfm    = J(k, 1, .)
    ts_bfm2   = J(k, 1, .)
    Om_bfm    = J(1+k+r, 1+k+r, .)
    cv_bfm    = .

    cupfm_fmbeta(y_wide_T1, X_long_T1, X_wide_T1, aik,
                 F_T1, dF, u1_wide_T1,
                 du2_T1_long, bw, N, T1, k, r,
                 beta_bfm, beta_bfm2, ts_bfm, ts_bfm2, Om_bfm, cv_bfm)

    real matrix beta_cur, beta_cur2, ts_cur, ts_cur2, Om_cur
    real scalar cv_cur
    beta_cur  = beta_bfm
    beta_cur2 = beta_bfm2
    ts_cur    = ts_bfm
    ts_cur2   = ts_bfm2
    Om_cur    = Om_bfm
    cv_cur    = cv_bfm

    real matrix _b1, _b2, _t1, _t2, _Om
    real scalar _cv
    _b1 = J(k, 1, .)
    _b2 = J(k, 1, .)
    _t1 = J(k, 1, .)
    _t2 = J(k, 1, .)
    _Om = J(1+k+r, 1+k+r, .)
    _cv = .

    itr = 1
    do {
        uhat_long = y_long - X_long*beta_cur
        U_wide    = cupfm_long2wide(uhat_long, N, T, 1)
        cupfm_factor(U_wide, r, F_hat, L_hat)
        aik = cupfm_delta(L_hat)

        F_T1 = F_hat[1::T1, .]
        dF   = cupfm_fdif(F_hat)

        FL_wide = J(T, N, 0)
        for (i=1; i<=N; i++) {
            FL_wide[., i] = F_hat * L_hat[i, .]'
        }
        FL_long    = cupfm_wide2long(FL_wide, N, T, 1)
        u1_wide    = cupfm_long2wide(y_long - X_long*beta_cur - FL_long, N, T, 1)
        u1_wide_T1 = u1_wide[1::T1, .]

        cupfm_fmbeta(y_wide_T1, X_long_T1, X_wide_T1, aik,
                     F_T1, dF, u1_wide_T1,
                     du2_T1_long, bw, N, T1, k, r,
                     _b1, _b2, _t1, _t2, _Om, _cv)
        beta_cur  = _b1
        beta_cur2 = _b2
        ts_cur    = _t1
        ts_cur2   = _t2
        Om_cur    = _Om
        cv_cur    = _cv

        if (verbose) printf("CupFM iter: %g\n", itr)
        itr++
    } while (itr <= maxiter)

    beta_cupfm  = beta_cur
    beta_cupfm2 = beta_cur2
    ts_cupfm    = ts_cur
    ts_cupfm2   = ts_cur2
    Om_cupfm    = Om_cur
    cv_cupfm    = cv_cur
    niter_cupfm = itr - 1

    uhat_long = y_long - X_long*beta_lsdv
    U_wide    = cupfm_long2wide(uhat_long, N, T, 1)

    real matrix F_bc, L_bc
    F_bc = J(T, r, .)
    L_bc = J(N, r, .)
    cupfm_factor(U_wide, r, F_bc, L_bc)

    itr_bc = 1
    real matrix beta_bc
    beta_bc = J(k, 1, .)
    do {
        beta_bc   = cupfm_plsbeta(X_long, y_long, F_bc, N, T, k)
        uhat_long = y_long - X_long*beta_bc
        U_wide    = cupfm_long2wide(uhat_long, N, T, 1)
        cupfm_factor(U_wide, r, F_bc, L_bc)
        if (verbose) printf("CupBC iter: %g\n", itr_bc)
        itr_bc++
    } while (itr_bc <= maxiter)

    aik_bc  = cupfm_delta(L_bc)
    dF_bc   = cupfm_fdif(F_bc)
    F_T1_bc = F_bc[1::T1, .]

    FL_wide_bc = J(T, N, 0)
    for (i=1; i<=N; i++) {
        FL_wide_bc[., i] = F_bc * L_bc[i, .]'
    }
    FL_long_bc    = cupfm_wide2long(FL_wide_bc, N, T, 1)
    u1_wide_bc    = cupfm_long2wide(y_long - X_long*beta_bc - FL_long_bc, N, T, 1)
    u1_wide_T1_bc = u1_wide_bc[1::T1, .]
    X_wide_T1_bc  = cupfm_long2wide(X_long_T1, N, T1, k)

    real matrix _b1bc, _b2bc, _t1bc, _t2bc, _Ombc
    real scalar _cvbc
    _b1bc = J(k, 1, .)
    _b2bc = J(k, 1, .)
    _t1bc = J(k, 1, .)
    _t2bc = J(k, 1, .)
    _Ombc = J(1+k+r, 1+k+r, .)
    _cvbc = .

    cupfm_fmbeta(y_wide_T1, X_long_T1, X_wide_T1_bc, aik_bc,
                 F_T1_bc, dF_bc, u1_wide_T1_bc,
                 du2_T1_long, bw, N, T1, k, r,
                 _b1bc, _b2bc, _t1bc, _t2bc, _Ombc, _cvbc)

    beta_cupbc = _b1bc
    ts_cupbc   = _t1bc
    cv_cupbc   = _cvbc
    Om_cupbc   = _Ombc

    st_matrix("_cupfm_b_lsdv",   beta_lsdv')
    st_matrix("_cupfm_b_baifm",  beta_bfm')
    st_matrix("_cupfm_b_cupfm",  beta_cupfm')
    st_matrix("_cupfm_b_cupfm2", beta_cupfm2')
    st_matrix("_cupfm_b_cupbc",  beta_cupbc')
    st_matrix("_cupfm_t_lsdv",   t_lsdv')
    st_matrix("_cupfm_t_baifm",  ts_bfm')
    st_matrix("_cupfm_t_cupfm",  ts_cupfm')
    st_matrix("_cupfm_t_cupfm2", ts_cupfm2')
    st_matrix("_cupfm_t_cupbc",  ts_cupbc')
    st_matrix("_cupfm_omega",    Om_cupfm)
    st_matrix("_cupfm_omega_bc", Om_cupbc)
    st_matrix("_cupfm_f",        F_hat)
    st_matrix("_cupfm_lambda",   L_hat)
    st_matrix("_cupfm_aik",      aik)
    st_numscalar("_cupfm_r",         r)
    st_numscalar("_cupfm_bw",        bw)
    st_numscalar("_cupfm_niter",     niter_cupfm)
    st_numscalar("_cupfm_cvar",      cv_cupfm)
    st_numscalar("_cupfm_cvar_bc",   cv_cupbc)
    st_numscalar("_cupfm_N",         N)
    st_numscalar("_cupfm_T",         T)
    st_numscalar("_cupfm_Nobs",      N*T)
    st_numscalar("_cupfm_converged", (niter_cupfm < maxiter ? 1 : 0))
}

end
// ── END MATA ENGINE ────────────────────────────────────────────────────────

capture program drop cupfm
program define cupfm, eclass
    version 14

    // --- Syntax ---------------------------------------------------------------
    // STATA 17 BATCH FIX: integer() and real() options cause r(197) "invalid
    // syntax" in Stata 17 batch mode. ALL numeric options are declared as
    // string() and converted manually below. This does NOT affect any calling
    // syntax - users still write nfactors(2), bandwidth(5) etc. as normal.
    // Also: LEVEL() is reserved by Stata eclass - do NOT declare it.
    // Also: fv removed from varlist - cupfm uses numeric panel data only.
    syntax varlist(min=2 numeric ts) [if] [in], ///
        [NFactors(string)              /// integer: 0 = auto-select via Bai-Ng (2002) IC
         BANDwidth(string)             /// integer: Bartlett kernel bandwidth (default=5)
         KERnel(string)                /// string:  bartlett|parzen (default: bartlett)
         MAXiter(string)               /// integer: max CupFM/CupBC iterations (default=20)
         TOLerance(string)             /// real:    convergence tolerance (default=0.0001)
         AUTORMAX(string)              /// integer: max r in auto-selection (default=8)
         NOICsummary                  /// suppress factor diagnostics table
         PLot                         /// produce all visualization plots
         PLotcoef                     /// plot: coefficient comparison only
         PLotfactors                  /// plot: estimated factors only
         PLotloadings                 /// plot: factor loadings only
         EXPort(string)               /// excel|latex|csv|all
         SAVing(string)               /// filename (no extension) for export/plots
         NOIsily                      /// verbose Mata output
         ]

    // --- Convert string options -> numeric & apply defaults -------------------
    if "`nfactors'"  == "" local nfactors  0
    else                    local nfactors  = real("`nfactors'")
    if "`bandwidth'" == "" local bandwidth 5
    else                    local bandwidth = real("`bandwidth'")
    if "`maxiter'"   == "" local maxiter   20
    else                    local maxiter   = real("`maxiter'")
    if "`tolerance'" == "" local tolerance 0.0001
    else                    local tolerance = real("`tolerance'")
    if "`autormax'"  == "" local autormax  8
    else                    local autormax  = real("`autormax'")
    local level 95   // hardcoded; LEVEL is reserved by Stata eclass

    // --- Validate numeric conversions -----------------------------------------
    foreach opt in nfactors bandwidth maxiter autormax {
        if missing(``opt'') {
            di as error "Option `opt'() must be a valid integer."
            exit 198
        }
        local `opt' = int(``opt'')
    }
    if missing(`tolerance') {
        di as error "Option tolerance() must be a valid real number."
        exit 198
    }

    // --- Sample ---------------------------------------------------------------
    marksample touse
    markout `touse' `varlist'

    // --- Parse varlist --------------------------------------------------------
    local depvar    : word 1 of `varlist'
    local indepvars : list varlist - depvar
    local nv        : word count `indepvars'

    if `nv' < 1 {
        di as error "cupfm requires at least one independent variable."
        exit 102
    }

    // --- Validate -------------------------------------------------------------
    _cupfm_validate `varlist' if `touse', ///
        nfactors(`nfactors') bandwidth(`bandwidth') ///
        maxiter(`maxiter') tolerance(`tolerance')

    local N_g      = r(N_g)
    local T        = r(T)
    local N_obs    = r(N_obs)
    local rmax     = r(rmax)
    local panelvar = "`r(panelvar)'"
    local timevar  = "`r(timevar)'"

    // --- Factor number logic --------------------------------------------------
    local do_autoR = 0
    local r_use    = `nfactors'
    if `nfactors' == 0 {
        local do_autoR = 1
        local r_use    = 1   // Mata will update this
    }

    // --- Kernel ---------------------------------------------------------------
    if "`kernel'" == "" local kernel "bartlett"
    if !inlist("`kernel'", "bartlett", "parzen") {
        di as error "kernel() must be 'bartlett' or 'parzen'."
        exit 198
    }

    // --- Sort data ------------------------------------------------------------
    sort `panelvar' `timevar'

    // --- DISPLAY: Running message ---------------------------------------------
    di ""
    di as text "  {bf:cupfm} - Panel Cointegration with Common Factors"
    di as text "  Running `=`N_g'' units x `=`T'' periods = `=`N_obs'' observations"
    if `do_autoR' {
        di as text "  Factor number: auto-select via Bai-Ng (2002) IC (rmax=`autormax')"
    }
    else {
        di as text "  Factor number: specified r = `r_use'"
    }
    di as text "  Bandwidth: `bandwidth' (`kernel' kernel)"
    di as text "  Max iterations: `maxiter'"
    di as text "  Estimating..."
    di ""

    // Re-sort (dataset sort order must be panelvar timevar for Mata engine)
    sort `panelvar' `timevar'

    // --- CALL MATA ENGINE - results stored in _cupfm_* global matrices --------
    mata: cupfm_main("`depvar'", "`indepvars'", "`touse'", ///
                     `N_g', `T', `r_use', `bandwidth', ///
                     "`kernel'", `maxiter', `do_autoR', `autormax', ///
                     `=cond("`noisily'"!="",1,0)')

    // --- Transfer _cupfm_ Stata scalars to local macros ----------------------
    local r_final     = _cupfm_r
    local niter_cupfm = _cupfm_niter
    local cvar_cupfm  = _cupfm_cvar
    local cvar_cupbc  = _cupfm_cvar_bc
    local converged   = _cupfm_converged

    // --- DISPLAY TABLES -------------------------------------------------------
    _cupfm_display, ///
        depvar(`depvar') indepvars(`indepvars') ///
        panelvar(`panelvar') timevar(`timevar') ///
        ng(`N_g') tperiods(`T') nobs(`N_obs') ///
        rfactors(`r_final') bwidth(`bandwidth') mxiter(`maxiter') ///
        niter(`niter_cupfm') ///
        `noicsummary'

    // --- PLOTS ----------------------------------------------------------------
    // IMPORTANT: must run BEFORE ereturn matrix (which DROPS _cupfm_b_* globals)
    if "`plot'" != "" | "`plotcoef'" != "" | "`plotfactors'" != "" | "`plotloadings'" != "" {
        local graph_opt ""
        if "`plotcoef'" != ""     local graph_opt "coef"
        if "`plotfactors'" != ""  local graph_opt "factors"
        if "`plotloadings'" != "" local graph_opt "loadings"
        if "`plot'" != ""         local graph_opt "all"

    // Pre-build option strings - avoids nested-quote r(132) in cond() expressions
        local saving_opt ""
        if "`saving'" != "" local saving_opt "saving(`saving')"
        local graph_arg "graph(all)"
        if "`graph_opt'" != "" local graph_arg "graph(`graph_opt')"

        _cupfm_plot, ///
            depvar(`depvar') indepvars(`indepvars') ///
            timevar(`timevar') ng(`N_g') tobs(`T') rfact(`r_final') ///
            `saving_opt' `graph_arg'
    }

    // --- EXPORT ---------------------------------------------------------------
    // IMPORTANT: must run BEFORE ereturn matrix (which DROPS _cupfm_b_* globals)
    if "`export'" != "" {
        if "`saving'" == "" {
            local saving "cupfm_results"
        }
        _cupfm_export, ///
            depvar(`depvar') indepvars(`indepvars') ///
            ng(`N_g') tobs(`T') rfact(`r_final') bwuse(`bandwidth') ///
            niter(`niter_cupfm') ///
            format(`export') saving(`saving')
    }

    // --- STORE IN e() ---------------------------------------------------------
    tempname b_cupfm V_cupfm

    // Build coefficient row vector for e(b): CupFM is the primary estimator
    matrix `b_cupfm' = _cupfm_b_cupfm

    // Build diagonal variance matrix: V_jj = (b_j / t_j)^2 = se_j^2
    matrix `V_cupfm' = J(`nv', `nv', 0)
    forvalues j = 1/`nv' {
        local bj = _cupfm_b_cupfm[1, `j']
        local tj = _cupfm_t_cupfm[1, `j']
        if abs(`tj') > 0.0001 {
            matrix `V_cupfm'[`j', `j'] = (`bj'/`tj')^2
        }
        else {
            matrix `V_cupfm'[`j', `j'] = 0
        }
    }

    // Name b and V with IDENTICAL column names - required by ereturn post
    matrix colnames `b_cupfm' = `indepvars'
    matrix rownames `V_cupfm' = `indepvars'
    matrix colnames `V_cupfm' = `indepvars'

    // Post primary result - clears previous e()
    ereturn post `b_cupfm' `V_cupfm', esample(`touse')

    // Additional coefficient matrices
    ereturn matrix beta_lsdv   = _cupfm_b_lsdv
    ereturn matrix beta_baifm  = _cupfm_b_baifm
    ereturn matrix beta_cupfm  = _cupfm_b_cupfm
    ereturn matrix beta_cupfm2 = _cupfm_b_cupfm2
    ereturn matrix beta_cupbc  = _cupfm_b_cupbc

    ereturn matrix tstat_lsdv   = _cupfm_t_lsdv
    ereturn matrix tstat_baifm  = _cupfm_t_baifm
    ereturn matrix tstat_cupfm  = _cupfm_t_cupfm
    ereturn matrix tstat_cupfm2 = _cupfm_t_cupfm2
    ereturn matrix tstat_cupbc  = _cupfm_t_cupbc

    // Factor/rotation matrices
    ereturn matrix F_hat    = _cupfm_f
    ereturn matrix Lambda   = _cupfm_lambda
    ereturn matrix Aik      = _cupfm_aik
    ereturn matrix Omega    = _cupfm_omega
    ereturn matrix Omega_bc = _cupfm_omega_bc

    // Scalars
    ereturn scalar Nobs        = `N_obs'
    ereturn scalar Ng          = `N_g'
    ereturn scalar Tperiods    = `T'
    ereturn scalar nfactors    = `r_final'
    ereturn scalar bw          = `bandwidth'
    ereturn scalar maxiter     = `maxiter'
    ereturn scalar niter       = `niter_cupfm'
    ereturn scalar cvar_cupfm  = `cvar_cupfm'
    ereturn scalar cvar_cupbc  = `cvar_cupbc'
    ereturn scalar level       = `level'
    ereturn scalar converged   = `converged'

    // Strings
    ereturn local cmd       "cupfm"
    ereturn local cmdline   "cupfm `0'"
    ereturn local depvar    "`depvar'"
    ereturn local indepvars "`indepvars'"
    ereturn local panelvar  "`panelvar'"
    ereturn local timevar   "`timevar'"
    ereturn local kernel    "`kernel'"
    ereturn local estimator "CupFM"
    ereturn local papers    "Bai, Kao & Ng (2009, JoE); Bai & Kao (2005, SSRN)"
    ereturn local vcetype   "Mixed normal (BKN 2009 Theorem 3)"

    di as text "  {it:Note: Primary result e(b), e(V) correspond to CupFM.}"
    di as text "  {it:All 5 estimators available in e(beta_*) and e(tstat_*).}"
    di ""
end