*! xthpool 1.0.0  11jul2026
*! Hausman poolability test for cointegrated panels (Westerlund & Hess, 2011)
*! Author: Merwan Roudane  (merwanroudane920@gmail.com)
*! https://github.com/merwanroudane
*
* Implements the maximum-Hausman poolability test of
*   Westerlund, J. and Hess, W. (2011) "A New Poolability Test for Cointegrated
*   Panels", Journal of Applied Econometrics 26: 56-88. doi:10.1002/jae.1143
*
* Step -> equation map (full derivation in xthpool_methods.sthlp):
*   [A] within (FE) demeaning of y_it, x_it        eq (1),(6); x°=x-mean_i(x)
*   [B] individual & pooled LS estimators          Sec 2.2 (beta_i, beta_pool)
*   [C] PCA of residual matrix ee' -> f_hat,lam,u  Sec 2.2 two-step, eq p.63
*   [D] IC1 factor-number selection                eq IC1(r), Prop 1
*   [E] Newey-West LR covariance of z=(f,u,dx)     Sec 2.2, Bartlett K(j)
*   [F] bias terms U_fvi,U_uvi,U_i; beta_i+,beta_pool+   p.60-61
*   [G] individual Hausman H_i (T^2 form)          p.61, Lemma A.1
*   [H] Hmax, Zmax normalization, Gumbel p-value   Thm 1, Remarks 1-3
*   [I] defactoring on observable g_t              Corollary 1, eq (4),(5)
*   [J] iterative sequential-drop scheme           Sec 3.2, crit = alpha^j pctile

program define xthpool, rclass
    version 14.0

    syntax varlist(min=2 numeric ts) [if] [in] , ///
        [ Factors(integer -1)     ///
          Rmax(integer 5)         ///
          Bandwidth(integer -1)   ///
          DEFACTor(varlist numeric ts) ///
          AConst(real 2)          ///
          Level(cilevel)          ///
          ITERate                 ///
          GRAPH                   ///
          name(string) ]

    // ---- parse dependent / independent variables --------------------------
    gettoken depvar indepvars : varlist
    local k : word count `indepvars'
    if (`k' < 1) {
        di as error "at least one regressor is required"
        exit 198
    }

    // ---- panel structure --------------------------------------------------
    qui xtset
    if ("`r(panelvar)'" == "") {
        di as error "data are not xtset; use {bf:xtset panelvar timevar} first"
        exit 459
    }
    local ivar "`r(panelvar)'"
    local tvar "`r(timevar)'"
    if ("`tvar'" == "") {
        di as error "a time variable is required; use {bf:xtset panelvar timevar}"
        exit 459
    }

    marksample touse
    markout `touse' `depvar' `indepvars' `defactor'
    qui replace `touse' = 0 if missing(`ivar') | missing(`tvar')

    // ---- require a balanced panel with no gaps ----------------------------
    tempvar nobs Tmin Tmax
    qui bysort `touse' `ivar' (`tvar'): gen long `nobs' = _N if `touse'
    sort `ivar' `tvar'
    qui su `nobs' if `touse', meanonly
    local Tbal = r(max)
    qui su `nobs' if `touse'
    if (r(min) != r(max)) {
        di as error "xthpool requires a balanced panel (each unit the same T, no gaps)"
        exit 459
    }

    local kg : word count `defactor'
    local dofact = ("`defactor'" != "")

    // ---- Newey-West Bartlett HAC bandwidth: default M = floor(4 (T/100)^(2/9))
    //      (the Newey-West rule); a user-supplied bandwidth() overrides it.
    local bw = `bandwidth'
    if (`bw' < 0) local bw = floor(4*(`Tbal'/100)^(2/9))
    if (`bw' < 1) local bw = 1

    // ---- run the engine ----------------------------------------------------
    // sort so that panelsetup blocks line up
    sort `ivar' `tvar'
    local doiter = ("`iterate'" != "")
    tempname RES HI IDV ITR BPLUS
    mata: _xthpool_engine( ///
        "`depvar'", "`indepvars'", "`defactor'", "`ivar'", "`tvar'", "`touse'", ///
        `factors', `rmax', `bw', `dofact', `aconst', `doiter', ///
        "`RES'", "`HI'", "`IDV'", "`ITR'", "`BPLUS'")

    // pull scalars back from the 1 x 9 results matrix
    tempname Hmax Zmax pval Nact Tt kk rr
    scalar `Hmax' = `RES'[1,1]
    scalar `Zmax' = `RES'[1,2]
    scalar `pval' = `RES'[1,3]
    local imax   = `RES'[1,4]
    local Nact   = `RES'[1,5]
    local Tt     = `RES'[1,6]
    local kk     = `RES'[1,7]
    local rr     = `RES'[1,8]
    local bwu    = `RES'[1,9]

    // ------------------------------------------------------------------ table
    di ""
    di as text "Hausman poolability test for cointegrated panels"
    di as text "Westerlund & Hess (2011, J. Appl. Econ.)"
    di as text "{hline 64}"
    di as text "Panel variable : " as result "`ivar'"
    di as text "Time variable  : " as result "`tvar'"
    di as text "Dependent var. : " as result "`depvar'"
    di as text "Regressors     : " as result "`indepvars'"
    if (`dofact') {
        di as text "Defactored on  : " as result "`defactor'"
    }
    di as text "{hline 64}"
    di as text "N (units) = " as result %5.0f `Nact' ///
       as text "   T = " as result %5.0f `Tt' ///
       as text "   m = " as result %2.0f `kk'
    di as text "Common factors r = " as result %2.0f `rr' ///
       as text "   NW bandwidth = " as result %3.0f `bwu' ///
       as text "   a_N = " as result %4.2f `aconst'
    di as text "{hline 64}"
    di as text %-22s "Statistic" _col(24) %12s "Value" _col(40) %12s "p-value"
    di as text "{hline 64}"
    local star ""
    if (`pval' < 0.10) local star "*"
    if (`pval' < 0.05) local star "**"
    if (`pval' < 0.01) local star "***"
    di as text %-22s "H_max" _col(24) as result %12.4f `Hmax' _col(46) " "
    di as text %-22s "Z_max (normalized)" _col(24) as result %12.4f `Zmax'
    di as text %-22s "Gumbel p-value" _col(24) as result %12.4f `pval' ///
       "  " as result "`star'"
    di as text "{hline 64}"
    di as text "H0: the slope is common across all units (poolable panel)"
    di as text "Most extreme unit (argmax H_i): " as result "`ivar' = `imax'"
    di as text "Significance: * 10%, ** 5%, *** 1%. p-value from Gumbel(0,1)."

    // ------------------------------------------------------- iterative scheme
    if ("`iterate'" != "") {
        di ""
        di as text "Iterative (sequential-drop) poolability scheme"
        di as text "{hline 64}"
        di as text %-4s "Step" _col(6) %-12s "Unit" _col(20) %10s "H_max" ///
           _col(32) %10s "Z_max" _col(44) %10s "p (raw)" _col(56) %10s "p (adj)"
        di as text "{hline 64}"
        local nsteps = rowsof(`ITR')
        forvalues s = 1/`nsteps' {
            local u  = `ITR'[`s',1]
            local hm = `ITR'[`s',2]
            local zm = `ITR'[`s',3]
            local pr = `ITR'[`s',4]
            local pa = `ITR'[`s',5]
            di as text %-4.0f `s' _col(6) as result %-12.0g `u' ///
               _col(20) %10.4f `hm' _col(32) %10.4f `zm' ///
               _col(44) %10.4f `pr' _col(56) %10.4f `pa'
        }
        di as text "{hline 64}"
        di as text "p (adj) maintains the overall level across steps"
        di as text "(critical value at step j = upper alpha^j Gumbel percentile)."
        return matrix itertable = `ITR'
    }

    // ------------------------------------------------------------- graph(s)
    if ("`graph'" != "") {
        if ("`name'" == "") local name xthpool
        _xthpool_graph "`HI'" "`IDV'" `Zmax' `Hmax' `aconst' `Nact' `kk' "`name'"
    }

    // --------------------------------------------------------- stored results
    return scalar Hmax  = `Hmax'
    return scalar Zmax  = `Zmax'
    return scalar p     = `pval'
    return scalar N     = `Nact'
    return scalar T     = `Tt'
    return scalar m     = `kk'
    return scalar r     = `rr'
    return scalar bw    = `bwu'
    return scalar aN    = `aconst'
    return scalar imax  = `imax'
    return local  depvar    "`depvar'"
    return local  indepvars "`indepvars'"
    return local  ivar      "`ivar'"
    return local  tvar      "`tvar'"
    return local  cmd       "xthpool"
    return matrix Hi     = `HI'
    return matrix unit   = `IDV'
    return matrix bplus  = `BPLUS'
end

*-------------------------------------------------------------------------------
* Graph helper: sorted individual H_i against the (per-unit) critical level
*-------------------------------------------------------------------------------
program define _xthpool_graph
    args HI IDV Zmax Hmax aN Nn kk nm
    // critical H that would make the maximum reject at 5%:
    // H_crit = a_N * z_.05 + b_N,  z_.05 = -ln(-ln(.95)) = 2.9702
    local zc = -ln(-ln(0.95))
    local bN = invchi2(`kk', 1 - 1/`Nn')
    local hc = `aN'*`zc' + `bN'
    preserve
    clear
    qui svmat double `HI', name(hi)
    qui svmat double `IDV', name(unit)
    qui gen long _ord = _n
    gsort -hi1
    qui gen _rank = _n
    capture {
        twoway (bar hi1 _rank, barwidth(0.7) color(navy%70)) ///
               (function y = `hc', range(_rank) lcolor(cranberry) lpattern(dash)), ///
            graphregion(color(white)) plotregion(color(white)) ///
            ytitle("Individual Hausman statistic H{sub:i}") ///
            xtitle("Units ordered by H{sub:i}") ///
            title("Poolability: individual Hausman statistics", size(medium)) ///
            legend(order(1 "H{sub:i}" 2 "5% rejection threshold") ///
                   rows(1) region(lstyle(none))) ///
            name("`nm'", replace)
    }
    if (_rc) {
        di as text "(graph skipped: `=_rc')"
    }
    restore
end

*-------------------------------------------------------------------------------
* Mata engine
*-------------------------------------------------------------------------------
version 14.0
mata:
mata set matastrict off

// --- Newey-West two-sided long-run covariance of the columns of Z ----------
// returns Delta (=Omega+Lam+Lam') and Lam (one-sided lagged sum) by reference
void _xthp_nw(real matrix Z, real scalar bw, real scalar T,
              real matrix Delta, real matrix Lam)
{
    real scalar Tz, j, kk2
    real matrix Om, A2
    real scalar kf

    Tz  = rows(Z)
    kk2 = cols(Z)
    Om  = (Z' * Z) / T
    Lam = J(kk2, kk2, 0)
    for (j = 1; j <= bw; j++) {
        if (j < Tz) {
            kf = 1 - j/(1 + bw)
            A2 = Z[(1)::(Tz-j), .]' * Z[(1+j)::(Tz), .]
            Lam = Lam + kf * (A2 / T)
        }
    }
    Delta = Om + Lam + Lam'
}

// --- core test on a given set of active columns of the residual matrix ------
// returns rowvector (Hmax, Zmax, praw, imax_pos, na, r_used) and, by ref,
// the full per-unit Hi vector (Hi) and the argmax position (pos)
real rowvector _xthp_core(real matrix ehat, pointer(real matrix) rowvector pXcz,
    pointer(real matrix) rowvector pDx, pointer(real matrix) rowvector pSdxxc,
    pointer(real matrix) rowvector pMinv, real matrix bi, real matrix xcy,
    pointer(real matrix) rowvector pMi,
    real colvector active, real scalar T, real scalar k,
    real scalar bw, real scalar aN, real scalar factors, real scalar rmax,
    real matrix Hi, real matrix Bout)
{
    real scalar na, r, rr, j, i, Vr, ic, best, bN, Hmax, Zmax, praw, Devi
    real matrix ea, A, Vfull, fhat, lam, uhat, uu, fh, lm
    real rowvector Lvals
    real matrix Z, Delta, Lam, Gam, fz, uz, dxi, Dvv, Dvvinv, Dfv, Duv, Dff, Duu
    real matrix Gfv, Guv, Gvv, mid, Ufv, Uuv, Ui, Upool, Mpool, Mpoolinv
    real matrix Bplus
    real colvector lami, biplus, bpool, bpoolplus, diff, xcysum, Devec
    real scalar fidx1, fidx2, uidx, vidx1, vidx2
    real scalar imaxpos
    real colvector idxo

    na = rows(active)
    ea = ehat[., active]                       // T x na  (individual FE resid)

    // ----- principal components (eigen of ea*ea', done once) ---------------
    A = ea * ea'
    symeigensystem(A, Vfull, Lvals)            // Vfull cols = eigenvectors
    // sort columns by eigenvalue descending so [.,1..r] are the r largest
    idxo  = order(Lvals', -1)
    Vfull = Vfull[., idxo]

    // choose number of factors
    if (factors >= 0) {
        r = factors
    }
    else {
        best = .
        r = 0
        for (rr = 0; rr <= rmax; rr++) {
            if (rr == 0) {
                uu = ea
            }
            else {
                fh = sqrt(T) :* Vfull[., (1)::(rr)]
                lm = (1/T) * (fh' * ea)
                uu = ea - fh * lm
            }
            Vr = sum(uu :^ 2) / (na * T)
            ic = log(Vr) + rr * log(na*T/(na+T)) * (na+T)/(na*T)
            if (rr == 0) {
                best = ic
                r = 0
            }
            else {
                if (ic < best) {
                    best = ic
                    r = rr
                }
            }
        }
    }
    if (r > na - 1) r = na - 1
    if (r < 0) r = 0

    if (r > 0) {
        fhat = sqrt(T) :* Vfull[., (1)::(r)]    // T x r
        lam  = (1/T) * (fhat' * ea)             // r x na
        uhat = ea - fhat * lam                  // T x na
    }
    else {
        fhat = J(T, 0, 0)
        lam  = J(0, na, 0)
        uhat = ea
    }

    // block indices in z = ( f(1..r), u(r+1), v(r+2..r+1+k) )
    fidx1 = 1
    fidx2 = r
    uidx  = r + 1
    vidx1 = r + 2
    vidx2 = r + 1 + k

    // ----- pass 1: per-unit bias terms, conditional LR var, unit Dvv --------
    Bplus    = J(na, k, 0)
    Devec    = J(na, 1, 0)
    Upool = J(1, k, 0)
    Mpool = J(k, k, 0)
    xcysum = J(k, 1, 0)

    for (j = 1; j <= na; j++) {
        i   = active[j]
        dxi = *pDx[i]                            // Tz x k  (first-diff of x)
        // align factor/idiosyncratic to t = 2..T
        if (r > 0) {
            fz = fhat[(2)::(T), .]               // Tz x r
        }
        else {
            fz = J(T-1, 0, 0)
        }
        uz = uhat[(2)::(T), j]                    // Tz x 1
        Z  = (fz, uz, dxi)                        // Tz x (r+1+k)

        _xthp_nw(Z, bw, T, Delta, Lam)

        // one-sided long-run covariance Gamma = Omega + Lambda = Delta - Lambda'
        // (the correction terms in Westerlund-Hess are one-sided; the driftless
        //  limit on p.79 forces contemporaneous + lagged, not lagged alone)
        Gam    = Delta - Lam'

        Dvv    = Delta[(vidx1)::(vidx2), (vidx1)::(vidx2)]   // two-sided, for scaling
        Dvvinv = invsym(Dvv)
        Duu    = Delta[uidx, uidx]
        Duv    = Delta[uidx, (vidx1)::(vidx2)]    // 1 x k
        Gvv    = Gam[(vidx1)::(vidx2), (vidx1)::(vidx2)]     // one-sided
        Guv    = Gam[uidx, (vidx1)::(vidx2)]      // 1 x k

        mid = (*pSdxxc[i]) - T * Gvv              // k x k

        // idiosyncratic bias piece (always present)
        Uuv = Duv * Dvvinv * mid + T * Guv        // 1 x k
        Ui  = Uuv

        // conditional LR variance of e given v (uses two-sided Delta)
        Devi = Duu - Duv * Dvvinv * Duv'          // scalar (u.v part)

        if (r > 0) {
            Dfv = Delta[(fidx1)::(fidx2), (vidx1)::(vidx2)]   // r x k
            Dff = Delta[(fidx1)::(fidx2), (fidx1)::(fidx2)]   // r x r
            Gfv = Gam[(fidx1)::(fidx2), (vidx1)::(vidx2)]     // r x k, one-sided
            Ufv = Dfv * Dvvinv * mid + T * Gfv               // r x k
            lami = lam[., j]                                  // r x 1
            Ui   = lami' * Ufv + Uuv                          // 1 x k
            // conditional var adds factor part
            Devi = lami' * (Dff - Dfv*Dvvinv*Dfv') * lami + Devi
        }

        biplus = (bi[i, .])' - (Ui * (*pMinv[i]))'           // k x 1
        Bplus[j, .]    = biplus'
        Devec[j]       = Devi

        Upool  = Upool + Ui
        Mpool  = Mpool + (*pMi[i])
        xcysum = xcysum + xcy[., i]
    }

    // ----- pooled bias-adjusted estimator ----------------------------------
    Mpoolinv  = invsym(Mpool)
    bpool     = Mpoolinv * xcysum                 // k x 1
    bpoolplus = bpool' - Upool * Mpoolinv          // 1 x k

    // ----- pass 2: individual Hausman statistics ---------------------------
    // Individual Hausman statistic.  Var(beta_i+) = Delta_e.vi * M_i^-1 with the
    // REALIZED information M_i = Sum x° x°' (the paper's 1/6 Delta_vv is only the
    // probability limit of the normalized moment M_i/T^2, E[int B° B°'] = 1/6
    // Delta_vv, used to state the closed-form asymptotic variance -- Theorem 1 /
    // Remark 2 require H_i to be asymptotically chi-square(m)).  Plugging the
    // realized information gives the pivotal cointegration Wald form:
    //     H_i = (beta_i+ - beta_pool+)' M_i (beta_i+ - beta_pool+) / Delta_e.vi
    Hi = J(na, 1, 0)
    for (j = 1; j <= na; j++) {
        i    = active[j]
        diff = (Bplus[j, .] - bpoolplus)'          // k x 1
        Hi[j] = (diff' * (*pMi[i]) * diff) / Devec[j]
    }

    // ----- Hmax, normalization, Gumbel p-value -----------------------------
    Hmax = max(Hi)
    imaxpos = 1
    for (j = 1; j <= na; j++) {
        if (Hi[j] >= Hmax) imaxpos = j
    }
    bN   = invchi2(k, 1 - 1/na)
    Zmax = (Hmax - bN) / aN
    praw = 1 - exp(-exp(-Zmax))

    Bout = Bplus
    return((Hmax, Zmax, praw, imaxpos, na, r))
}

// --- top-level entry -------------------------------------------------------
void _xthpool_engine(string scalar yv, string scalar xv, string scalar gv,
    string scalar idv, string scalar tv, string scalar tousev,
    real scalar factors, real scalar rmax, real scalar bw, real scalar dofact,
    real scalar aN, real scalar iterate,
    string scalar RESn, string scalar HIn, string scalar IDVn,
    string scalar ITRn, string scalar BPLUSn)
{
    real matrix Y, X, G, info, ehat, bi, xcy, Bplus, itout
    real colvector id, tt, idvals, active
    real scalar N, T, k, kg, i, r1, r2, ni, T_i
    real scalar Hmax, Zmax, praw, imaxpos, na, rused, step, uid
    pointer(real matrix) rowvector pXcz, pDx, pSdxxc, pMinv, pMi
    real matrix Yi, Xi, Gi, Xc, xcz, dxi, Mi, Miinv, res, out, Bout
    real colvector yci, meanx, meany, bihat, Hi
    real rowvector cr

    Y  = st_data(., yv, tousev)
    X  = st_data(., xv, tousev)
    id = st_data(., idv, tousev)
    tt = st_data(., tv, tousev)
    k  = cols(X)
    kg = 0
    if (dofact) {
        G  = st_data(., gv, tousev)
        kg = cols(G)
    }

    info   = panelsetup(id, 1)
    N      = rows(info)
    idvals = J(N, 1, .)
    T      = info[1,2] - info[1,1] + 1

    ehat = J(T, N, 0)
    bi   = J(N, k, 0)
    xcy  = J(k, N, 0)
    pXcz   = J(1, N, NULL)
    pDx    = J(1, N, NULL)
    pSdxxc = J(1, N, NULL)
    pMinv  = J(1, N, NULL)
    pMi    = J(1, N, NULL)

    // ---- per-unit setup ([A]-[C] pieces that do not depend on active set) --
    for (i = 1; i <= N; i++) {
        r1 = info[i,1]
        r2 = info[i,2]
        idvals[i] = id[r1]
        Yi = Y[(r1)::(r2), .]
        Xi = X[(r1)::(r2), .]

        // [I] defactoring on observable common factors g_t
        if (dofact) {
            Gi = G[(r1)::(r2), .]
            Xi = Xi - Gi * (invsym(Gi'Gi) * (Gi'Xi))
        }

        // [A] within demeaning
        meany = mean(Yi)
        meanx = mean(Xi)
        yci = Yi :- meany
        Xc  = Xi :- meanx

        // [B] individual within estimator
        Mi    = Xc' * Xc
        Miinv = invsym(Mi)
        bihat = Miinv * (Xc' * yci)               // k x 1
        bi[i, .]  = bihat'
        xcy[., i] = Xc' * yci                       // k x 1

        // [C] individual FE residual (for PCA)
        ehat[., i] = yci - Xc * bihat

        // first difference of the (defactored) regressor, aligned t=2..T
        dxi = Xi[(2)::(T), .] - Xi[(1)::(T-1), .]   // Tz x k
        xcz = Xc[(2)::(T), .]                        // Tz x k (x° for t=2..T)

        // store pointers to FRESH temporaries (":+ 0" forces a distinct copy);
        // &namedvar would alias the reused loop variable -> all units would
        // collapse to the last one (Mata pointer gotcha).
        pMi[i]    = &(Mi :+ 0)
        pMinv[i]  = &(Miinv :+ 0)
        pDx[i]    = &(dxi :+ 0)
        pXcz[i]   = &(xcz :+ 0)
        pSdxxc[i] = &(dxi' * xcz)                    // k x k  Sum dx (x°)'
    }

    // ---- base test on the full panel --------------------------------------
    active = (1::N)
    cr = _xthp_core(ehat, pXcz, pDx, pSdxxc, pMinv, bi, xcy, pMi,
                    active, T, k, bw, aN, factors, rmax, Hi, Bout)
    Hmax = cr[1]; Zmax = cr[2]; praw = cr[3]
    imaxpos = cr[4]; na = cr[5]; rused = cr[6]
    uid = idvals[active[imaxpos]]

    st_matrix(RESn, (Hmax, Zmax, praw, uid, na, T, k, rused, bw))
    st_matrix(HIn, Hi)
    st_matrix(IDVn, idvals[active])
    st_matrix(BPLUSn, (idvals[active], Bout))

    // ---- iterative sequential-drop scheme [J] -----------------------------
    if (iterate) {
        itout = J(0, 5, .)
        active = (1::N)
        step   = 1
        while (rows(active) >= 2) {
            cr = _xthp_core(ehat, pXcz, pDx, pSdxxc, pMinv, bi, xcy, pMi,
                            active, T, k, bw, aN, factors, rmax, Hi, Bout)
            Hmax = cr[1]; Zmax = cr[2]; praw = cr[3]; imaxpos = cr[4]
            uid  = idvals[active[imaxpos]]
            // adjusted p to maintain overall level: p_adj = p_raw^(1/step)
            itout = itout \ (uid, Hmax, Zmax, praw, praw^(1/step))
            // drop the maximizing unit
            active = select(active, (1::rows(active)) :!= imaxpos)
            step = step + 1
        }
        st_matrix(ITRn, itout)
    }
    else {
        st_matrix(ITRn, J(1, 5, 0))
    }
}

end