*! version 1.0.0 12jul2026

Title

xtmixedroot methods -- Formulas, assumptions and the
code-to-paper map for xtmixedroot

The model

For units i = 1,...,N and periods t = 1,...,T,

y_it = lambda_i + u_it, u_it = alpha_i u_it-1 + e_it,

where alpha_i is in [0,1], e_it is iid over i and t, and lambda_i is a
unit-specific intercept independent of alpha_i and e_it (Ng's A1-A3).
A fraction theta = N1/N of the units has alpha_i = 1 ("I(1) units"); the
rest have alpha_i < 1 ("I(0) units"). The object of interest is theta,
with null hypotheses H0: theta = theta_0 for theta_0 in (0,1].
Extensions: unit-specific AR(p) dynamics and heteroskedastic errors
(Estimator A), a common stationary factor y_it = lambda_i F_t + u_it with
beta(L)F_t = v_t (Estimator B, Ng's A4-A5), and incidental trends
y_it = lambda_i t + u_it (Estimator C).

Why the cross-sectional variance identifies theta (Ng, Lemma 1)

Let V_t = (1/N) sum_i (y_it - ybar_t)^2 be the cross-sectional variance.
For an I(0) unit the variance of y_it settles at a constant; for an I(1)
unit it grows like t. Decomposing the mixture variance into within- and
between-group parts, Ng shows that for t beyond a short transient

V_t ~ constant + theta * t, so Delta V_t -> theta.

The between-group term (Ybar1_t - Ybar0_t)^2 is a squared random walk whose
drift and variance vanish at rates 1/N1 and t/N1^2, which is why V_t is
first-differenced (a levels regression of V_t on t would be spurious for
fixed N). Time-averaging Delta V_t gives Ng's Theorem 1:

sqrt(N) (theta-hat - theta) -> N(0, 2*theta), as N -> inf, then T -> inf.

The rate is sqrt(N), not sqrt(NT), because var(Delta V_t) grows linearly in
t, so averaging over time does not accelerate convergence.

Estimator A (heterogeneous dynamics, heteroskedastic errors)

For each unit i:

1. OLS on y_it = a_0 + a_1 y_it-1 + ... + a_p y_it-p + e_it (Ng eq.
(8)). p is fixed by lags(#) or selected unit-by-unit by the BIC over
1..maxlags on a common estimation sample, then the chosen model is refit on
the unit's full sample.

2. sigma_i^2-hat = (1/T) x the sum of squared residuals (Ng step 2 verbatim).

3. phi_i1-hat,...,phi_ip-hat = the inverse roots of the estimated AR
polynomial, computed as the eigenvalues of its companion matrix and ordered
by modulus (phi_i1 is the dominant root).

4. D_i-hat = product over k = 2..p of (phi_i1 - phi_ik). This is the
partial-fractions denominator attached to the dominant root: for a series
with phi_i1 = 1, var(y_it) = (sigma_i^2/D_i^2) t + lower-order terms, so the
rescaled series D_i y_it / sigma_i has variance slope 1 and the
cross-sectional variance of the rescaled panel has slope theta.
xtmixedroot uses |D_i| (only D_i^2 matters) and floors it at 1e-6
against near-repeated roots. For p = 1, D_i = 1.

5. V-hat_t = cross-sectional variance of the rescaled panel;
theta-hat = mean of Delta V-hat_t over the T-1 available increments;
eta-hat_t = Delta V-hat_t - theta-hat. The variance of theta-hat is
estimated by the Bartlett/Newey-West long-run variance of eta-hat_t with
truncation M (hac(#)), divided by the number of increments:

se(theta-hat)^2 = [gamma(0) + 2 sum_ (1 - s/(M+1)) gamma(s)] / (T-1).

The studentized statistic t = (theta-hat - theta_0)/se is treated as N(0,1).
Ng does not state Estimator A as a theorem (nuisance parameters are
estimated) but verifies by simulation that the normal approximation is
accurate; xtmixedroot's shipped validation reproduces her Table 1.

Estimator B (cross-sectional correlation)

Identical to A except that step 1 augments the regression with proxies for
the common stationary factor (Ng eq. (12)):

y_it = a_0 + a_1 y_it-1 + ... + a_p y_it-p + l_0 F-hat_t + ... + l_q F-hat_t-q + e_it,

where F-hat_t is the cross-sectional average of the differenced data
(Pesaran-type proxy; default, with q = factors(#) lags) or the first
principal component of the differenced panel (pc). Under a stationary
factor, Delta F_t^2 acts as a stationary error component in Delta V_t, so
theta-hat remains consistent even without the augmentation - but
augmenting sharpens the estimates of the dynamic parameters and removes the
finite-sample downward bias Ng documents in her Table 3.

Estimator C (incidental trends)

With y_it = lambda_i t + u_it the cross-sectional variance acquires a
quadratic term var(lambda_i) t^2. Step 1 adds a linear trend to the AR
regression; step 4 rescales by 1/sigma_i-hat (exactly as printed in Ng's
Estimator C; see the reconciliation notes below); and theta is estimated as
the intercept of the OLS regression (Ng eq. (13))

Delta V-hat_t = theta + beta [t^2 - (t-1)^2] + eta_t,

with Newey-West HAC standard errors for both coefficients.
beta-hat estimates var(lambda_i) (returned in r(varlambda)).
Ng's Theorem 2 gives sqrt(N/T) (theta-hat - theta) -> N(0,
32 theta var(lambda_i)/15): the rate deteriorates from sqrt(N) to
sqrt(N/T), so consistency requires T/N -> 0 and very large panels
(her Table 4 uses T = 300-600, N = 200-400).

Westerlund (2016) fixed-T statistics

Westerlund replaces Ng's sequential asymptotics with a finite-sample
expansion valid for any T >= 2 and derives the exact fixed-T bias and
variance. With theta-hat_w = (1/T) sum_ Delta V_t computed on the
raw data (V_T - V_1, telescoped):

theta* = theta-hat_w / sigma_eps^2-hat, theta*_BA = theta* + theta_0/T,

since E(theta* - theta_0) = -theta_0/T exactly under the null. The
statistics, all asymptotically N(0,1) as N -> inf:

tau*_1,T = sqrt(N)(theta*_BA - 1)/sigma-hat_theta,T with
sigma-hat_theta,T^2 = 2(T^2-1)/T^2 + (T-1)/T^2 (kappa-hat - 3). Valid for
any T >= 2 under homogeneous intercepts (sigma_lambda^2 = 0); a common time
effect is handled by demean. Left-tailed test of theta = 1.

tau*_1,NT = same numerator with sigma-hat_theta,NT^2 =
2(T^2-1)/T^2 + (T-1)/T^2 [4 sigma-hat_lambda^2/sigma-hat_eps^2 +
kappa-hat - 3], where sigma-hat_lambda^2 is the cross-sectional variance of
the intercepts from unit-wise regressions of y_it on (1, y_it-1). Requires
larger T for sigma-hat_lambda^2 to be consistent.

tau*_1 = sqrt(N)(theta*_BA - 1)/sqrt(2), the large-(N,T) form.

tau*_theta0 = sqrt(N)(theta*_BA - theta_0)/sqrt(2 theta*) for
theta_0 in (0,1) (Westerlund Theorem 2), two-sided; requires N,T -> inf
jointly with sqrt(N)/T -> 0. The sigma_theta,NT-variance variant is also
reported. Under theta_0 < 1, sigma_eps^2 and kappa are estimated from the
unit-wise AR(1) residuals rather than from Delta y (Delta y is a valid
proxy for e only under theta = 1).

Moment estimators: sigma-hat_eps^2 and kappa-hat are the pooled second and
fourth moments of Delta y_it (t = 2..T), normalized by the number of
increments N(T-1); see the reconciliation notes.

Local power and identification (what a rejection means)

Westerlund models the stationary units as alpha_i = exp(c_i/N^eta), c_i <= 0.
Key results: (i) under fixed alternatives (eta = 0) power goes to 1 as
N -> inf and theta is estimable; (ii) under local alternatives (eta > 0)
theta-hat converges to 1 regardless of theta - the fraction is not
identified, and only H0: theta = 1 is testable; (iii) with T fixed the
test has nontrivial local power exactly at eta = 1/2; with T = N^gamma
growing, at eta = gamma + 1/2, so a larger T buys power against closer
alternatives. Practical reading: in a short panel, failing to reject
theta = 1 - or estimating theta-hat near 1 - may only mean the stationary
units are very persistent, and having a large T is necessary but not
sufficient to pin down theta (the deviation from the null must also not be
too small).

Step-by-equation map (code block -> paper)

code step source

V_t = (1/N) sum (y_it - ybar_t)^2 Ng Sec. 2 (definition of V_t,N); Westerlund Sec. 3.1 (V_t)

theta-hat_w = (V_T - V_1)/T Westerlund's feasible estimator, sum over t = 2..T

per-unit OLS AR(p), BIC lags Ng Estimator A step 1, eq. (8); BIC per Ng Sec. 4 and applications

sigma_i^2-hat = mean sq. residual Ng Estimator A step 2

phi-hat = companion eigenvalues Ng Estimator A step 3 ("reciprocal of the roots of alpha_i(L) = 0")

D_i = prod_k>=2 (phi_i1 - phi_ik) Ng p. 117, the D_i = A^-1 B construction (see reconciliations)

rescale D_i y/sigma_i; V-hat_t Ng Estimator A step 4

theta-hat = mean Delta V-hat_t; HAC se Ng Estimator A step 5, Bartlett kernel K(s,M) = 1 - s/(M+1)

factor proxy dybar_t, lags 0..q Ng Estimator B step 1, eq. (12); Pesaran cross-sectional average

first-PC proxy (pc) Ng Sec. 3.2 ("first principal component of Delta y_it")

trend in AR; Delta V on (1, 2t-1) Ng Estimator C steps 1 and 5, eq. (13)

t tests at theta_0, .01, 1 Ng Sec. 4, hypotheses A, B, C

classification by |phi_i1| rank Ng Sec. 4 (informal proposal; ~60 percent correct rate)

s2_eps, kappa from Delta y Westerlund Sec. 3.1.1

theta*_BA = theta* + theta_0/T Westerlund Remark 4 / Theorem 1

sigma_theta,T^2; tau*_1,T Westerlund Sec. 3.1.1

sigma_theta,NT^2; tau*_1,NT; tau*_1 Westerlund Sec. 3.1.2

residual-based s2_eps under theta_0<1 Westerlund Sec. 3.2

tau*_theta0, two-sided Westerlund Theorem 2 and following display

demean option Westerlund Remark 1 and footnote 3

Reconciliations with the printed papers

Three places where the printed text is internally inconsistent were resolved
in favor of the papers' own defining constructions and simulation tables;
each choice is verified by the shipped Monte Carlo validation.

1. sigma_eps^2-hat denominator (Westerlund). The text prints
sigma-hat_eps^2 = sum_i sum_ (Delta y_it)^2 / NT, but the sum has
N(T-1) terms and the /NT normalization would shift the mean of tau*_1,T by
sqrt(N)/T (e.g. +8.9 at T = 2, N = 320), contradicting the near-zero means
in his Table 1. xtmixedroot divides by N(T-1), which makes
E(theta*_BA) = 1 exact under the null; the shipped validation reproduces his
Table 1 size (empirical 3.0 percent vs his 3.5 at T = 4, N = 160).

2. D_i table (Ng, p. 117). The rows for p = 3 and p = 4 of the
convenience table contain typographical errors (a sign on the phi_2 phi_3
term at p = 3; a dropped product in the last term at p = 4). Solving Ng's
own defining system D_i = A^-1 B gives the dominant-root partial-fractions
weight A_i1/D_i = phi_i1^(p-1) / prod_(phi_i1 - phi_ik), which matches
the printed table exactly at p = 1, 2 and matches A^-1 B numerically for
all p. xtmixedroot therefore uses D_i = prod_(phi_i1 -
phi_ik).

3. Estimator C rescaling (Ng, p. 123). Step 4 of Estimator C
rescales by 1/sigma_i only (no D_i), although step 3 still computes D_i-hat.
xtmixedroot follows the printed step 4 verbatim.

Two further implementation notes. (i) Ng's step 5 prints the studentized
statistic as (theta-hat - theta)/omega with omega^2 the long-run variance of
eta_t; the variance of a sample mean is that long-run variance divided by
the number of increments, and only this normalization yields a t statistic
with unit variance - confirmed by the validation, which reproduces
N var(theta-hat)/theta = 3.2 against Ng's Table 1 value 2.996. (ii) V_0 is
unobservable, so the time average runs over the T-1 observable increments
(Westerlund's feasible convention); for the Ng block the average divides by
T-1, which is the exact sample mean of the increments.

Validation

The ancillary file xtmixedroot_example.do regenerates the papers' own
Monte Carlo designs and checks every code path. Headline results (seeds
fixed in the file):

quantity (design) xtmixedroot vs paper

mean theta-hat, Ng Model 1a (N=30, T=100, theta=.5) .515 vs .517

N var(theta-hat)/theta, same design 3.22 vs 2.996

5 percent two-sided rejection of the true null .125 vs .114

size of tau*_1,T (theta=1, T=4, N=160) .030 vs .035

est. A vs B means under a common factor (N=60, T=100) .421/.452 vs .459/.481

References

Ng, S. 2008. A simple test for nonstationarity in mixed panels.
Journal of Business & Economic Statistics 26(1): 113-127.
https://doi.org/10.1198/073500106000000675

Westerlund, J. 2016. A simple test for nonstationarity in mixed
panels: A further investigation. Journal of Statistical Planning and
Inference 173: 1-30.
https://doi.org/10.1016/j.jspi.2016.01.004

Author

Merwan Roudane

Email: merwanroudane920@gmail.com

GitHub: https://github.com/merwanroudane