*! xtmixedroot 1.0.0  12jul2026
* Author: Merwan Roudane (merwanroudane920@gmail.com)
* https://github.com/merwanroudane
*
* Estimates theta, the fraction of I(1) (unit-root) units in a mixed panel,
* and tests H0: theta = theta0 for any theta0 in (0,1].
*
* Methods implemented (see help xtmixedroot_methods for the full
* step -> equation map):
*   - Ng, S. (2008) "A Simple Test for Nonstationarity in Mixed Panels",
*     Journal of Business & Economic Statistics 26(1), 113-127.
*     Estimators A (heterogeneous AR(p) dynamics), B (cross-sectional
*     correlation via factor proxies) and C (incidental trends), with
*     HAC (Newey-West) inference on theta-hat.
*   - Westerlund, J. (2016) "A simple test for nonstationarity in mixed
*     panels: A further investigation", Journal of Statistical Planning
*     and Inference 173, 1-30. Bias-adjusted statistics tau*_{1,T}
*     (valid for any T >= 2), tau*_{1,NT}, tau*_1 and tau*_{theta0}.
*
* Abbreviated step -> equation map (full map in the methods help):
*   V_t (cross-sec variance)          Ng eq. before Lemma 1; Westerlund Sec 3.1
*   theta-hat = mean of dV_t          Ng Thm 1; Westerlund's feasible t=2..T sum
*   AR(p) per unit, sigma_i, phi_ij   Ng Estimator A steps 1-3, eq (8)
*   D_i = prod_{k>=2}(phi_i1-phi_ik)  Ng p.117 (D_i = A^-1 B construction)
*   rescale D_i*y/sigma_i             Ng Estimator A step 4
*   factor proxies dybar_t, lags      Ng Estimator B, eq (12); Pesaran proxy
*   trend + dV on (1, 2t-1)           Ng Estimator C, eq (13)
*   HAC se of theta-hat               Ng Estimator A step 5 (Bartlett kernel)
*   theta* = theta-hat/s2e            Westerlund Sec 3.1
*   theta*_BA = theta* + theta0/T     Westerlund bias adjustment (Remark 4)
*   sigma2_theta,T / sigma2_theta,NT  Westerlund Sec 3.1.1 / 3.1.2
*   tau*_theta0                       Westerlund Theorem 2

program define xtmixedroot, rclass
    version 14.0
    syntax varname(numeric) [if] [in] [, ESTimator(string) Lags(integer 0)  ///
        MAXLags(integer 4) Factors(integer 1) PC Theta0(real 1)             ///
        HAC(integer 2) DEMean CLASSify GRaph LIST(integer 15) LEVel(cilevel) ]

    * ------------------------------------------------ panel setup
    capture quietly xtset
    if _rc {
        di as err "data must be xtset; use {bf:xtset panelvar timevar}"
        exit 459
    }
    local ivar "`r(panelvar)'"
    local tvar "`r(timevar)'"
    if ("`tvar'" == "") {
        di as err "time variable not set; use {bf:xtset panelvar timevar}"
        exit 459
    }

    marksample touse
    markout `touse' `ivar' `tvar'
    quietly count if `touse'
    if (r(N) == 0) error 2000

    * ------------------------------------------------ options
    if ("`estimator'" == "") local estimator "a"
    local estimator = lower("`estimator'")
    if (!inlist("`estimator'", "a", "b", "c")) {
        di as err "estimator() must be a, b or c"
        exit 198
    }
    local estcode = cond("`estimator'" == "a", 1, cond("`estimator'" == "b", 2, 3))
    if (`theta0' <= 0 | `theta0' > 1) {
        di as err "theta0() must lie in (0,1]; theta = 0 is not testable " ///
            "(the limiting distribution collapses; test theta0(.01) instead)"
        exit 198
    }
    if (`lags' < 0) {
        di as err "lags() must be a positive integer"
        exit 198
    }
    if (`maxlags' < 1) {
        di as err "maxlags() must be at least 1"
        exit 198
    }
    if (`hac' < 0) {
        di as err "hac() must be a nonnegative integer"
        exit 198
    }
    if (`factors' < 0) {
        di as err "factors() must be a nonnegative integer"
        exit 198
    }
    local usepc = ("`pc'" != "")
    local dem   = ("`demean'" != "")

    * ------------------------------------------------ engine
    tempname RES VMAT UMAT
    mata: _xmr_engine("`varlist'", "`ivar'", "`tvar'", "`touse'",        ///
        `estcode', `lags', `maxlags', `factors', `usepc', `theta0',      ///
        `hac', `dem', "`RES'", "`VMAT'", "`UMAT'")

    * ------------------------------------------------ unpack results
    local theta   = `RES'[1,1]
    local se      = `RES'[1,2]
    local tng     = `RES'[1,3]
    local png     = `RES'[1,4]
    local tlow    = `RES'[1,5]
    local plow    = `RES'[1,6]
    local tall    = `RES'[1,7]
    local pall    = `RES'[1,8]
    local beta    = `RES'[1,9]
    local sebeta  = `RES'[1,10]
    local thetaw  = `RES'[1,11]
    local s2e     = `RES'[1,12]
    local kappa   = `RES'[1,13]
    local s2lam   = `RES'[1,14]
    local thstar  = `RES'[1,15]
    local thba    = `RES'[1,16]
    local tau1T   = `RES'[1,17]
    local ptau1T  = `RES'[1,18]
    local tau1NT  = `RES'[1,19]
    local ptau1NT = `RES'[1,20]
    local tau1    = `RES'[1,21]
    local ptau1   = `RES'[1,22]
    local tauth   = `RES'[1,23]
    local ptauth  = `RES'[1,24]
    local tauthNT = `RES'[1,25]
    local ptauthNT = `RES'[1,26]
    local NN      = `RES'[1,27]
    local TT      = `RES'[1,28]
    local N1hat   = `RES'[1,29]
    local ngdone  = `RES'[1,30]
    local pbar    = `RES'[1,31]

    * ------------------------------------------------ display
    local estlbl = cond(`estcode' == 1, "A (heterogeneous AR(p) dynamics)", ///
        cond(`estcode' == 2, "B (A + cross-sectional correlation control)", ///
        "C (A + incidental linear trends)"))
    if (`lags' > 0) local laglbl "fixed p = `lags'"
    else            local laglbl = "BIC, 1..`maxlags' (mean p = " + string(`pbar', "%4.2f") + ")"
    if (`estcode' == 2) {
        local flbl = cond(`usepc', "first PC of D.`varlist'", "cross-sectional mean of D.`varlist'")
        local flbl "`flbl', lags 0..`factors'"
    }

    di
    di as text "{bf:Fraction of nonstationary (I(1)) units in a mixed panel}"
    di as text "Ng (2008, JBES 26:113-127); Westerlund (2016, JSPI 173:1-30)"
    di
    di as text "Panel variable: " as res "`ivar'" as text "    Time variable: " ///
        as res "`tvar'" as text "    N = " as res `NN' as text "    T = " as res `TT'
    di as text "Series: " as res "`varlist'" ///
        as text cond(`dem', "  (cross-sectionally demeaned for variance estimation)", "")
    if (`ngdone') {
        di as text "Estimator: " as res "`estlbl'"
        di as text "Lag order: " as res "`laglbl'" as text "    HAC lags (Bartlett): " as res "`hac'"
        if (`estcode' == 2) di as text "Factor proxy: " as res "`flbl'"
    }

    if (`ngdone') {
        local alpha2 = (100 - `level')/200
        local zc = invnormal(1 - `alpha2')
        local lo = `theta' - `zc'*`se'
        local hi = `theta' + `zc'*`se'
        di
        di as text "{hline 78}"
        di as text %-42s "Ng (2008) point estimate" _col(46) %9s "Estimate" ///
            _col(57) %9s "Std.Err." _col(68) "[`level'% CI]"
        di as text "{hline 78}"
        di as text "  theta  (fraction of I(1) units)" _col(46) as res %9.4f `theta' ///
            _col(57) %9.4f `se' _col(67) %6.3f `lo' as text "," as res %6.3f `hi'
        if (`estcode' == 3) {
            di as text "  var(lambda_i)  (incidental-trend var.)" _col(46) as res %9.4f `beta' ///
                _col(57) %9.4f `sebeta'
        }
        di as text "  implied number of I(1) units" _col(46) as res %9.0f `N1hat' ///
            as text "  (= [theta x N], integer part)"
        di as text "{hline 78}"
        di as text "Hypothesis tests based on theta-hat (HAC t statistics, asymptotically N(0,1))"
        di as text "{hline 78}"
        local s1 = cond(`png' < .01, "***", cond(`png' < .05, "**", cond(`png' < .1, "*", "")))
        local s2 = cond(`plow' < .01, "***", cond(`plow' < .05, "**", cond(`plow' < .1, "*", "")))
        local s3 = cond(`pall' < .01, "***", cond(`pall' < .05, "**", cond(`pall' < .1, "*", "")))
        local side1 = cond(`theta0' == 1, "H1: theta < 1     ", "H1: theta != `theta0'")
        di as text "  H0: theta = " %5.3f `theta0' "   `side1'" _col(46) as res %9.3f `tng' ///
            _col(57) as text "p = " as res %6.4f `png' as text "  `s1'"
        di as text "  H0: theta = 0.010   H1: theta > 0.01  " _col(46) as res %9.3f `tlow' ///
            _col(57) as text "p = " as res %6.4f `plow' as text "  `s2'  (any I(1) units?)"
        di as text "  H0: theta = 1.000   H1: theta < 1     " _col(46) as res %9.3f `tall' ///
            _col(57) as text "p = " as res %6.4f `pall' as text "  `s3'  (all units I(1)?)"
        di as text "{hline 78}"
    }
    else {
        di
        di as text "{hline 78}"
        di as text "Note: T = " as res `TT' as text " is too small for the Ng (2008) estimator" ///
            " (requires T >= 25;"
        di as text "      T >= 100 recommended). Only the Westerlund (2016) fixed-T statistics"
        di as text "      are reported below."
        di as text "{hline 78}"
    }

    di
    di as text "Westerlund (2016) bias-adjusted statistics"
    di as text "{hline 78}"
    di as text "  theta-hat (feasible, raw data)" _col(46) as res %9.4f `thetaw'
    di as text "  theta*  = theta-hat / s2_eps" _col(46) as res %9.4f `thstar'
    di as text "  theta*_BA = theta* + theta0/T" _col(46) as res %9.4f `thba'
    di as text "  s2_eps = " as res %7.4f `s2e' as text "   kappa_eps = " as res %7.4f `kappa' ///
        as text "   s2_lambda = " as res %7.4f `s2lam'
    if (`theta0' < 1) {
        di as text "  (s2_eps and kappa_eps from unit-wise AR(1) residuals; Westerlund Sec. 3.2)"
    }
    di as text "{hline 78}"
    if (`theta0' == 1) {
        local w1 = cond(`ptau1T' < .01, "***", cond(`ptau1T' < .05, "**", cond(`ptau1T' < .1, "*", "")))
        local w2 = cond(`ptau1NT' < .01, "***", cond(`ptau1NT' < .05, "**", cond(`ptau1NT' < .1, "*", "")))
        local w3 = cond(`ptau1' < .01, "***", cond(`ptau1' < .05, "**", cond(`ptau1' < .1, "*", "")))
        di as text "  Tests of H0: theta = 1 vs H1: theta < 1 (left-tailed)"
        di as text "  tau*_1,T   (any T >= 2; homog. intercepts)" _col(46) as res %9.3f `tau1T' ///
            _col(57) as text "p = " as res %6.4f `ptau1T' as text "  `w1'"
        di as text "  tau*_1,NT  (heterog. intercepts; large T) " _col(46) as res %9.3f `tau1NT' ///
            _col(57) as text "p = " as res %6.4f `ptau1NT' as text "  `w2'"
        di as text "  tau*_1     (large N and T)                " _col(46) as res %9.3f `tau1' ///
            _col(57) as text "p = " as res %6.4f `ptau1' as text "  `w3'"
    }
    else {
        local w1 = cond(`ptauth' < .01, "***", cond(`ptauth' < .05, "**", cond(`ptauth' < .1, "*", "")))
        local w2 = cond(`ptauthNT' < .01, "***", cond(`ptauthNT' < .05, "**", cond(`ptauthNT' < .1, "*", "")))
        di as text "  Tests of H0: theta = " %5.3f `theta0' " vs H1: theta != " %5.3f `theta0' ///
            " (two-sided; large N,T)"
        di as text "  tau*_theta0    (Theorem 2, simplified)   " _col(46) as res %9.3f `tauth' ///
            _col(57) as text "p = " as res %6.4f `ptauth' as text "  `w1'"
        di as text "  tau*_theta0,NT (sigma_theta,NT variance) " _col(46) as res %9.3f `tauthNT' ///
            _col(57) as text "p = " as res %6.4f `ptauthNT' as text "  `w2'"
    }
    di as text "{hline 78}"
    di as text "  * p<0.10, ** p<0.05, *** p<0.01." ///
        "  See {helpb xtmixedroot_methods:help xtmixedroot methods}."
    if (`ngdone' & `TT' < 100) {
        di as text "  Note: T < 100; Ng (2008) recommends T >= 100 for theta-hat. The"
        di as text "  Westerlund fixed-T statistics remain valid for any T >= 2."
    }
    if (`thstar' <= 0 & `thstar' < .) {
        di as text "  Note: theta* <= 0: the cross-sectional variance is non-increasing,"
        di as text "  which suggests theta ~ 0 (no I(1) units). tau*_theta0 is unavailable"
        di as text "  in this case; use the one-sided test of H0: theta = 0.01."
    }

    * ------------------------------------------------ classification
    if ("`classify'" != "" & `ngdone') {
        local nl = min(`list', `NN')
        di
        di as text "Classification of units (Ng 2008, Sec. 4): the [theta x N] = " ///
            as res `N1hat' as text " units"
        di as text "with the largest estimated dominant root |phi_i1| are classified I(1)."
        di as text "{hline 60}"
        di as text %12s "`ivar'" _col(16) %10s "|phi_i1|" _col(30) %10s "sigma_i" ///
            _col(44) %4s "p_i" _col(52) "class"
        di as text "{hline 60}"
        forvalues k = 1/`nl' {
            local uid   = `UMAT'[`k',1]
            local uphi  = `UMAT'[`k',2]
            local usig  = `UMAT'[`k',3]
            local up    = `UMAT'[`k',5]
            local ucls  = `UMAT'[`k',6]
            local clab = cond(`ucls' == 1, "I(1)", "I(0)")
            di as res %12.0g `uid' _col(16) %10.4f `uphi' _col(30) %10.4f `usig' ///
                _col(44) %4.0f `up' _col(52) "`clab'"
        }
        if (`nl' < `NN') di as text "  ... (" `NN' - `nl' " more units; full list in r(units))"
        di as text "{hline 60}"
    }

    * ------------------------------------------------ graph
    if ("`graph'" != "") {
        preserve
        quietly {
            clear
            set obs `TT'
            svmat double `VMAT', name(__xmrv)
            * __xmrv1 = time, __xmrv2 = raw V_t, __xmrv3 = rescaled V_t
            label variable __xmrv2 "cross-sectional variance"
            label variable __xmrv1 "`tvar'"
        }
        if (`ngdone') {
            quietly {
                gen double __idx = _n
                local b2 = cond(`estcode' == 3, `beta', 0)
                gen double __fit = __xmrv3[1] + `theta'*(__idx - 1) + `b2'*(__idx^2 - 1)
                label variable __xmrv3 "rescaled variance"
                label variable __fit "implied trend (slope = theta-hat)"
            }
            twoway (line __xmrv2 __xmrv1, lcolor(navy) lwidth(medthick)),        ///
                graphregion(color(white)) name(__xmr_g1, replace) nodraw        ///
                title("Raw data", size(medsmall)) ytitle("V{subscript:t}")      ///
                xtitle("`tvar'")
            twoway (line __xmrv3 __xmrv1, lcolor(navy) lwidth(medthick))         ///
                (line __fit __xmrv1, lcolor(cranberry) lpattern(dash)),          ///
                graphregion(color(white)) name(__xmr_g2, replace) nodraw        ///
                title("Rescaled data", size(medsmall))                          ///
                ytitle("V{subscript:t} (rescaled)") xtitle("`tvar'")            ///
                legend(order(1 "variance" 2 "slope = theta-hat") rows(1) size(small))
            graph combine __xmr_g1 __xmr_g2, rows(1) graphregion(color(white))  ///
                title("Cross-sectional variance of `varlist'", size(medium))    ///
                note("The variance of a mixed panel trends upward at rate theta (Ng 2008, Lemma 1)", size(vsmall))
            capture graph drop __xmr_g1 __xmr_g2
        }
        else {
            twoway (line __xmrv2 __xmrv1, lcolor(navy) lwidth(medthick)),        ///
                graphregion(color(white))                                        ///
                title("Cross-sectional variance of `varlist'", size(medium))    ///
                ytitle("V{subscript:t}") xtitle("`tvar'")
        }
        restore
    }

    * ------------------------------------------------ stored results
    return scalar N        = `NN'
    return scalar T        = `TT'
    return scalar theta0   = `theta0'
    return scalar level    = `level'
    if (`ngdone') {
        return scalar theta    = `theta'
        return scalar se_theta = `se'
        return scalar t        = `tng'
        return scalar p        = `png'
        return scalar t_low    = `tlow'
        return scalar p_low    = `plow'
        return scalar t_all    = `tall'
        return scalar p_all    = `pall'
        return scalar N1hat    = `N1hat'
        return scalar p_mean   = `pbar'
        if (`estcode' == 3) {
            return scalar varlambda    = `beta'
            return scalar se_varlambda = `sebeta'
        }
    }
    return scalar theta_w    = `thetaw'
    return scalar theta_star = `thstar'
    return scalar theta_ba   = `thba'
    return scalar sigma2e    = `s2e'
    return scalar kappa      = `kappa'
    return scalar sigma2lam  = `s2lam'
    if (`theta0' == 1) {
        return scalar tau1T    = `tau1T'
        return scalar p_tau1T  = `ptau1T'
        return scalar tau1NT   = `tau1NT'
        return scalar p_tau1NT = `ptau1NT'
        return scalar tau1     = `tau1'
        return scalar p_tau1   = `ptau1'
    }
    else {
        return scalar tautheta0     = `tauth'
        return scalar p_tautheta0   = `ptauth'
        return scalar tautheta0NT   = `tauthNT'
        return scalar p_tautheta0NT = `ptauthNT'
    }
    return local estimator "`estimator'"
    return local varname   "`varlist'"
    return local ivar      "`ivar'"
    return local tvar      "`tvar'"
    return local cmd       "xtmixedroot"

    matrix colnames `VMAT' = time V Vscaled
    return matrix V = `VMAT'
    if (`ngdone') {
        matrix colnames `UMAT' = id phi1 sigma D p i1
        return matrix units = `UMAT'
    }
end

* =====================================================================
* Mata engine
* =====================================================================
mata:

void _xmr_engine(string scalar yv, string scalar iv, string scalar tv,
                 string scalar touse, real scalar estcode, real scalar pfix,
                 real scalar pmax, real scalar qf, real scalar usepc,
                 real scalar theta0, real scalar M, real scalar dem,
                 string scalar resn, string scalar vmn, string scalar umn)
{
    real matrix data, Y, Yd, DY, DYr, X, XX, XXi, A, EV, UM, VM, Vb, meat
    real colvector V, Vs, dV, eta, lam, phi1, Dvec, sig, pvec, cls, times, ids
    real colvector yy, e, b, ftv, f, u, tt, chk, dt, mods, idx2, ordp
    real rowvector evals, xt
    real scalar N, T, i, t, s, k, l, p, mv, nd, s2e1, kap1, s2e0, kap0, s2lam
    real scalar thetaw, ths1, ths0, ba1, ba0, s2T, s2NT, s2NT0, sqN
    real scalar tau1T, p1T, tau1NT, p1NT, tau1, p1, tauth, pth, tauthNT, pthNT
    real scalar theta, se, tng, png, tlow, plow, tall, pall, beta, sebeta
    real scalar ngdone, pbar, N1, minT, Td, g0, gs, om2, den, sse, s4, ml
    real scalar r0, r0c, n, sig2, bic, bestbic, bestp, w, D2, sf, mf
    real scalar thsdisp, badisp, s2edisp, kapdisp
    complex rowvector L
    complex scalar lam1, Dc

    // ---------- data assembly ----------
    data = st_data(., (iv, tv, yv), touse)
    data = sort(data, (1,2))
    ids   = uniqrows(data[.,1])
    times = uniqrows(data[.,2])
    N = rows(ids)
    T = rows(times)
    if (N < 5) {
        errprintf("xtmixedroot: at least 5 cross-sectional units are required\n")
        exit(2001)
    }
    if (T < 2) {
        errprintf("xtmixedroot: at least 2 time periods are required\n")
        exit(2001)
    }
    if (rows(data) != N*T) {
        errprintf("xtmixedroot: the panel must be balanced with no gaps over the estimation sample\n")
        exit(459)
    }
    chk = J(N,1,1) # times
    if (max(abs(data[.,2] - chk)) > 0) {
        errprintf("xtmixedroot: the panel must be balanced with no gaps over the estimation sample\n")
        exit(459)
    }
    if (T > 2) {
        dt = times[|2 \ T|] - times[|1 \ T-1|]
        if (max(dt) - min(dt) > 1e-8) {
            errprintf("xtmixedroot: the time variable must be equally spaced\n")
            exit(459)
        }
    }
    Y = colshape(data[.,3], T)'

    // ---------- cross-sectional variance of the raw data ----------
    V = J(T,1,.)
    for (t=1; t<=T; t=t+1) {
        mv = sum(Y[t,.])/N
        V[t] = sum((Y[t,.] :- mv):^2)/N
    }

    // ---------- Westerlund (2016) block ----------
    thetaw = (V[T] - V[1])/T
    Yd = Y
    if (dem) {
        Yd = Y :- (rowsum(Y):/N)
    }
    DY = Yd[|2,1 \ T,N|] - Yd[|1,1 \ T-1,N|]
    nd = N*(T-1)
    s2e1 = quadsum(DY:^2)/nd
    kap1 = .
    if (s2e1 > 0) {
        kap1 = quadsum(DY:^4)/(s2e1*s2e1*nd)
    }

    // per-unit AR(1) regressions: intercepts (sigma2_lambda) and
    // residual-based error moments (used when theta0 < 1)
    lam = J(N,1,.)
    s2e0 = .
    kap0 = .
    s2lam = .
    if (T >= 4) {
        sse = 0
        s4  = 0
        for (i=1; i<=N; i=i+1) {
            yy = Yd[|2,i \ T,i|]
            X  = J(T-1,1,1), Yd[|1,i \ T-1,i|]
            XX = quadcross(X,X)
            b  = invsym(XX)*quadcross(X,yy)
            e  = yy - X*b
            lam[i] = b[1]
            sse = sse + quadsum(e:^2)
            s4  = s4  + quadsum(e:^4)
        }
        s2e0 = sse/nd
        if (s2e0 > 0) {
            kap0 = s4/(s2e0*s2e0*nd)
        }
        ml = sum(lam)/N
        s2lam = sum((lam :- ml):^2)/N
    }

    sqN = sqrt(N)
    ths1 = .
    tau1T = .
    p1T = .
    tau1NT = .
    p1NT = .
    tau1 = .
    p1 = .
    if (s2e1 > 0) {
        ths1 = thetaw/s2e1
        ba1  = ths1 + 1/T
        s2T  = 2*(T*T - 1)/(T*T) + (T - 1)/(T*T)*(kap1 - 3)
        if (s2T > 0) {
            tau1T = sqN*(ba1 - 1)/sqrt(s2T)
            p1T   = normal(tau1T)
        }
        if (s2lam < .) {
            s2NT = 2*(T*T - 1)/(T*T) + (T - 1)/(T*T)*(4*s2lam/s2e1 + kap1 - 3)
            if (s2NT > 0) {
                tau1NT = sqN*(ba1 - 1)/sqrt(s2NT)
                p1NT   = normal(tau1NT)
            }
        }
        tau1 = sqN*(ba1 - 1)/sqrt(2)
        p1   = normal(tau1)
    }

    tauth = .
    pth = .
    tauthNT = .
    pthNT = .
    ths0 = .
    if (theta0 < 1 & s2e0 < . & s2e0 > 0) {
        ths0 = thetaw/s2e0
        ba0  = ths0 + theta0/T
        den  = 2*ths0
        if (den > 1e-12) {
            tauth = sqN*(ba0 - theta0)/sqrt(den)
            pth   = 2*normal(-abs(tauth))
        }
        if (kap0 < . & s2lam < .) {
            s2NT0 = 2*(T*T - 1)/(T*T) + (T - 1)/(T*T)*(4*s2lam/s2e0 + kap0 - 3)
            if (s2NT0 > 0) {
                tauthNT = sqN*(ba0 - theta0)/(sqrt(theta0)*sqrt(s2NT0))
                pthNT   = 2*normal(-abs(tauthNT))
            }
        }
    }

    // ---------- Ng (2008) block ----------
    ngdone = 0
    theta = .
    se = .
    tng = .
    png = .
    tlow = .
    plow = .
    tall = .
    pall = .
    beta = .
    sebeta = .
    pbar = .
    N1 = .
    Vs   = J(T,1,.)
    phi1 = J(N,1,.)
    Dvec = J(N,1,.)
    sig  = J(N,1,.)
    pvec = J(N,1,.)
    cls  = J(N,1,.)
    minT = 25

    if (T >= minT) {
        ngdone = 1

        // factor proxy for Estimator B (Ng eq (12); Pesaran-type proxy)
        ftv = J(T,1,.)
        if (estcode == 2) {
            DYr = Y[|2,1 \ T,N|] - Y[|1,1 \ T-1,N|]
            if (usepc) {
                // first principal component of the differenced panel
                A = DYr*DYr'/N
                symeigensystem(A, EV = ., evals = .)
                f = EV[.,1]
                mf = sum(f)/(T-1)
                sf = sqrt(sum((f :- mf):^2)/(T-1))
                if (sf > 0) {
                    f = f/sf
                }
            }
            else {
                f = rowsum(DYr)/N
            }
            ftv[|2 \ T|] = f
        }

        // per-unit AR(p) fits (Ng Estimator A steps 1-4)
        pbar = 0
        for (i=1; i<=N; i=i+1) {
            // lag order: fixed or per-unit BIC on the common sample
            if (pfix > 0) {
                p = pfix
            }
            else {
                r0c = pmax + 1
                if (estcode == 2 & r0c < qf + 2) {
                    r0c = qf + 2
                }
                bestbic = .
                bestp = 1
                for (k=1; k<=pmax; k=k+1) {
                    sig2 = _xmr_arfit(Y[.,i], ftv, k, r0c, T, estcode, qf, b = .)
                    if (sig2 <= 0) {
                        continue
                    }
                    n = T - r0c + 1
                    bic = ln(sig2) + (k + 1 + (estcode == 2)*(qf + 1) + (estcode == 3))*ln(n)/n
                    if (bic < bestbic) {
                        bestbic = bic
                        bestp = k
                    }
                }
                p = bestp
            }
            pvec[i] = p
            pbar = pbar + p/N

            // final fit on the unit's own full sample
            r0 = p + 1
            if (estcode == 2 & r0 < qf + 2) {
                r0 = qf + 2
            }
            sig2 = _xmr_arfit(Y[.,i], ftv, p, r0, T, estcode, qf, b = .)
            // Ng step 2 verbatim: sigma_i^2 = (1/T) sum of squared residuals
            // (_xmr_arfit returns SSR/n with n = T - r0 + 1 residuals)
            sig2 = sig2*(T - r0 + 1)/T
            if (sig2 < 1e-12) {
                sig2 = 1e-12
            }
            sig[i] = sqrt(sig2)

            // dominant root and rescaling constant D_i
            // D_i = prod_{k=2}^{p} (phi_i1 - phi_ik): the partial-fractions
            // denominator defined by Ng's D_i = A^-1 B system (the printed
            // table rows for p = 3,4 contain typos; see methods help)
            if (p == 1) {
                phi1[i] = abs(b[1])
                Dvec[i] = 1
            }
            else {
                A = b[|1 \ p|]' \ (I(p-1), J(p-1,1,0))
                L = eigenvalues(A)
                mods = abs(L)'
                idx2 = order(-mods, 1)
                lam1 = L[idx2[1]]
                phi1[i] = mods[idx2[1]]
                Dc = 1
                for (k=2; k<=p; k=k+1) {
                    Dc = Dc*(lam1 - L[idx2[k]])
                }
                D2 = abs(Dc)
                if (D2 < 1e-6) {
                    D2 = 1e-6
                }
                Dvec[i] = D2
            }
        }

        // rescale (Estimator A/B: D_i*y/sigma_i; Estimator C: y/sigma_i,
        // exactly as printed in Ng's Estimator C step 4)
        for (i=1; i<=N; i=i+1) {
            if (estcode == 3) {
                w = 1/sig[i]
            }
            else {
                w = Dvec[i]/sig[i]
            }
            Y[.,i] = w*Y[.,i]
        }
        for (t=1; t<=T; t=t+1) {
            mv = sum(Y[t,.])/N
            Vs[t] = sum((Y[t,.] :- mv):^2)/N
        }
        dV = Vs[|2 \ T|] - Vs[|1 \ T-1|]
        Td = T - 1

        if (estcode <= 2) {
            // theta-hat = time average of dV (Ng Estimator A step 5);
            // HAC (Bartlett) standard error of the sample mean
            theta = sum(dV)/Td
            eta = dV :- theta
            g0 = quadcross(eta, eta)/Td
            om2 = g0
            for (s=1; s<=M; s=s+1) {
                if (s < Td) {
                    gs = quadcross(eta[|s+1 \ Td|], eta[|1 \ Td-s|])/Td
                    om2 = om2 + 2*(1 - s/(M + 1))*gs
                }
            }
            if (om2 < 0) {
                om2 = g0
            }
            se = sqrt(om2/Td)
        }
        else {
            // Estimator C: dV_t = theta + beta*(t^2-(t-1)^2) + eta_t (Ng eq (13))
            tt = (2::T)
            X = J(Td,1,1), (2*tt :- 1)
            XX = quadcross(X,X)
            XXi = invsym(XX)
            b = XXi*quadcross(X,dV)
            theta = b[1]
            beta  = b[2]
            u = dV - X*b
            meat = J(2,2,0)
            for (t=1; t<=Td; t=t+1) {
                xt = X[t,.]
                meat = meat + (u[t]*u[t])*(xt'xt)
            }
            for (s=1; s<=M; s=s+1) {
                if (s < Td) {
                    w = 1 - s/(M + 1)
                    for (t=s+1; t<=Td; t=t+1) {
                        A = X[t,.]'X[t-s,.]
                        meat = meat + w*u[t]*u[t-s]*(A + A')
                    }
                }
            }
            Vb = XXi*meat*XXi
            se = sqrt(Vb[1,1])
            sebeta = sqrt(Vb[2,2])
        }

        // hypothesis tests (Ng Sec. 4, hypotheses A, B, C)
        if (se > 0) {
            tng = (theta - theta0)/se
            if (theta0 == 1) {
                png = normal(tng)
            }
            else {
                png = 2*normal(-abs(tng))
            }
            tlow = (theta - 0.01)/se
            plow = 1 - normal(tlow)
            tall = (theta - 1)/se
            pall = normal(tall)
        }

        // classification by the dominant root (Ng Sec. 4)
        // Ng Sec. 4 verbatim: N1-hat = [theta-hat x N], the integer part
        N1 = floor(theta*N)
        if (N1 < 0) {
            N1 = 0
        }
        if (N1 > N) {
            N1 = N
        }
        cls = J(N,1,0)
        ordp = order(-phi1, 1)
        for (k=1; k<=N1; k=k+1) {
            cls[ordp[k]] = 1
        }
    }

    // ---------- outputs ----------
    // theta*/theta*_BA reported at the requested theta0: dy-based moments
    // under theta0 = 1, AR(1)-residual-based moments under theta0 < 1
    thsdisp = ths1
    s2edisp = s2e1
    kapdisp = kap1
    if (theta0 < 1 & ths0 < .) {
        thsdisp = ths0
        s2edisp = s2e0
        kapdisp = kap0
    }
    badisp = .
    if (thsdisp < .) {
        badisp = thsdisp + theta0/T
    }
    UM = ids, phi1, sig, Dvec, pvec, cls
    if (ngdone) {
        UM = UM[order(-UM[.,2], 1), .]
    }
    VM = times, V, Vs
    st_matrix(vmn, VM)
    st_matrix(umn, UM)
    st_matrix(resn, (theta, se, tng, png, tlow, plow, tall, pall, beta,
        sebeta, thetaw, s2edisp, kapdisp, s2lam, thsdisp, badisp,
        tau1T, p1T, tau1NT, p1NT, tau1, p1, tauth, pth, tauthNT, pthNT,
        N, T, N1, ngdone, pbar))
}

// OLS AR(p) fit for one unit: y_t on (const, y_{t-1..t-p} [, f_t..f_{t-q}] [, t])
// over rows r0..T; returns the residual variance and (via b) the p lag
// coefficients in b[1..p]
real scalar _xmr_arfit(real colvector y, real colvector ftv, real scalar p,
                       real scalar r0, real scalar T, real scalar estcode,
                       real scalar qf, real colvector b)
{
    real matrix X, XX
    real colvector yy, e, bb
    real scalar n, l, sig2

    n = T - r0 + 1
    if (n < p + 4) {
        b = J(p,1,.)
        return(-1)
    }
    yy = y[|r0 \ T|]
    X = J(n, 0, .)
    for (l=1; l<=p; l=l+1) {
        X = X, y[|r0-l \ T-l|]
    }
    if (estcode == 2) {
        for (l=0; l<=qf; l=l+1) {
            X = X, ftv[|r0-l \ T-l|]
        }
    }
    if (estcode == 3) {
        X = X, (r0::T)
    }
    X = X, J(n,1,1)
    XX = quadcross(X,X)
    bb = invsym(XX)*quadcross(X,yy)
    e = yy - X*bb
    sig2 = quadsum(e:^2)/n
    b = bb[|1 \ p|]
    return(sig2)
}

end