*! version 1.0.0 18jul2026

Title

xtoutliers methods -- Equations, derivations, and code®paper map

This page documents the mathematics behind the three commands and maps every
computational step to the equation it implements. It is the companion to the
_COMPATIBILITY_MAP.md shipped with the package.

1. VSOM (xtvsom)

Source: Ismadyaliana, Setiawan & Purnomo, MethodsX 13 (2024) 102900;
Gumedze (2008); Thompson (1985).

Null model and standardized residual. The within (FE) or pooled
design gives residuals eit and variance
{&sigma}^2 = e'e/(nT-p). With leverage
vit = diagonal of the hat matrix (for FE, plus the
1/Ti mean-projection term),

t^2it = eit^2 / [ {&sigma}^2 (1 - vit) ]. (paper, boxed p.4)

Cutoff. Under the null there are no outliers, so the empirical
distribution of maxi t^2it is generated by a
parametric bootstrap: draw y* = X{&beta}^ + {&epsilon}*,
{&epsilon}* {&sim} N(0,{&sigma}^2), refit, recompute
max t^2. Because (I-H)X = 0, the standardized residual is
scale free, so {&sigma} cancels and the bootstrap needs no variance
estimate. The cutoff is the (1-{&alpha}) quantile of the bootstrap maxima.

Variance-shift estimates (Thompson 1985 REML, paper p.4, boxed), for
a flagged observation with t^2 > 1:

{&psi}it = (nT-p)(t^2-1) / [ (nT-p-t^2)(1-vit) ],

{&sigma}^2it = (nT-p-t^2)/(nT-p-1) {&times} {&sigma}^2.

Key equivalence (VSOM = WLS). The variance model is
var(y) = {&sigma}^2 (D {&Psi} D' + I) = {&sigma}^2 P.
Because D is a 0/1 indicator matrix selecting the outliers and
{&Psi} = diag({&psi}), the matrix P is diagonal, with
Pii = 1+{&psi}i for an outlier and 1 otherwise. Hence the
feasible-GLS refit

{&beta}^V = (X'P^-1X)^-1X'P^-1y

is exactly weighted least squares with weights
wi = 1/(1+{&psi}i). This is what xtvsom
computes.

Residual shrinkage. The fitted value carries the shift term
{&theta}^it = {&psi}it d'P^-1(y-X{&beta}^V).
Writing rit = yit-xit'{&beta}^V
(the WLS raw residual),

yit - xit'{&beta}^V - {&theta}^it = rit - [{&psi}/(1+{&psi})] rit = rit/(1+{&psi}it) = wit rit.

So the reported VSOM residual is wit rit, which shrinks
toward zero at outliers - reproducing the residual columns of the source
paper (e.g. 0.5129 -> 0.0760, giving w{&asymp}0.148,
{&psi}{&asymp}5.8).

Simultaneous (2SLS). For the simultaneous model the second-stage
design is X** = [Y^, Xexog, cons], with Y^ the
first-stage fitted endogenous regressors (paper Eq. 10-13). Everything above
then applies with X** in place of X.

2. Robust S-estimator and WLE (xtrobust)

Source: Jaseem & Mohammad, Al-Nahrain J. Science 27(4) (2024) 40-46.

Transforms. FE: within demeaning (paper Eq. 5-6). RE:
quasi-demeaning (Eq. 7)
{&lambda}i = 1 - [{&sigma}{&epsilon}^2/({&sigma}{&epsilon}^2 + Ti{&sigma}v^2)]^(1/2),
with {&sigma}{&epsilon}^2, {&sigma}v^2
from a preliminary xtreg,re.

S-estimator (Eq. 10-13). Tukey biweight
{&rho}(u) = u^2/2 - u^4/(2c^2) + u^6/(6c^4)
for |u|{&le}c, else c^2/6, with c=1.547. The scale
{&sigma}S solves mean {&rho}(ri/{&sigma}) = K,
K = 0.199 = (c^2/6)/2, giving a 50% breakdown point. The
estimator {&beta}^ = argmin {&sigma}S is minimised by
fast-S: nsamp elemental subsets, each refined by concentration steps,
retaining the smallest scale, then full IRLS.

WLE (Eq. 14-17). With residual density
f*(ri) (Gaussian kernel, bandwidth h) and smoothed model
density m*(ri) = N(0,{&sigma}v^2+h^2),
the Pearson residual is
{&delta}i = f*(ri)/m*(ri) - 1, the RAF
is A({&delta}) = 2[({&delta}+1)^(1/2) - 1], and the weight is
wi = min{1, [A({&delta}i)+1]+/({&delta}i+1)}.
Iterating {&beta} = (X'WX)^-1X'Wy solves the weighted
likelihood equation (Eq. 14).

3. Loss functions (xtlossf)

Source: Coleman & Bryan, arXiv:2509.07014v2 (2025).

Part I (nonnegative). From symmetry, monotonicity in the difference,
and decrease in the base, the simplest loss is Cobb-Douglas; a Lie symmetry
fixes the difference exponent to 1:

L = |F-B| B^q (Eq. I.2), signed S = (F-B) B^q (Eq. I.7).

q=0 is the absolute difference, q=-1 the absolute relative
difference; q=-1/2 is optimal for comparing two unbiased estimate sets with
variance proportional to the level (Eq. I.13-I.16). The time-invariant form
uses exponent tq+t-1 (Eq. I.6), t rescaled so the last date is 1.
Calibration from a criteria table solves log {&epsilon} = -q log B + K
(Eq. I.19-I.20), C = exp(K) (fit).

Part II (mixed sign). With {&Sigma} = |F|+|B|,

L = |F-B| {&Sigma}^q (Eq. II.5), signed S = (F-B) {&Sigma}^q (Eq. II.9),

with L(0,0)=0 (a continuity point when p+q>0, Thm II.1) and
q=-1 the harmonic-mean average absolute relative difference (Eq. II.6). An
outlier is declared when L > C; C is empirical (a quantile or a Tukey
IQR fence) - there is no distributional significance test, by construction.

4. Compatibility

Every step above is traced to a numbered equation, and the equivalences (VSOM =
WLS; residual = w r; K = c^2/6/2; smoothed model density =
N(0,{&sigma}^2+h^2)) are proved in
_COMPATIBILITY_MAP.md. Where a paper provides no quantity (RE-VSOM, a
significance test for a loss, a CI for the S scale) the command reports none and
the limitation is documented rather than invented.

Author

Dr Merwan Roudane

merwanroudane920@gmail.com

github.com/merwanroudane