{smcl}
{* *! version 1.0.0  17sep2026}{...}
{vieweralsosee "ltscoint" "help ltscoint"}{...}
{viewerjumpto "The model" "ltscoint_methods##model"}{...}
{viewerjumpto "The LTS estimator" "ltscoint_methods##lts"}{...}
{viewerjumpto "Assumptions and the oracle property" "ltscoint_methods##theory"}{...}
{viewerjumpto "Step-to-equation map" "ltscoint_methods##map"}{...}
{viewerjumpto "Inference after LTS" "ltscoint_methods##inference"}{...}
{viewerjumpto "Determining h" "ltscoint_methods##h"}{...}
{viewerjumpto "Decisions the paper leaves open" "ltscoint_methods##decisions"}{...}
{viewerjumpto "Validation" "ltscoint_methods##validation"}{...}
{title:Title}

{p2colset 5 26 28 2}{...}
{p2col:{bf:ltscoint methods} {hline 2}}Methods and formulas behind {helpb ltscoint}{p_end}
{p2colreset}{...}

{pstd}Equation numbers refer to Berenguer-Rico and Nielsen (2026), {it:Oxford Bulletin of Economics and
Statistics} 88, 690-711, unless stated otherwise.


{marker model}{...}
{title:The model}

{pstd}
For a scalar y_t and a (p-1)-vector z_t, with x_t = (y_t, z_t')', the ADL in equilibrium-correction form is

{p 12 12 2}(2.1)  D.y_t = omega' D.z_t + alpha (y_t-1 - kappa' z_t-1 - nu_c) + sum_{j=1}^{k-1} gamma_j' D.x_t-j + sigma eps_t,{p_end}

{pstd}or with a linear trend{p_end}

{p 12 12 2}(2.3)  D.y_t = omega' D.z_t + alpha (y_t-1 - kappa' z_t-1 - nu_l t) + sum_j gamma_j' D.x_t-j + mu_c + sigma eps_t.{p_end}

{pstd}
In vector form D.y_t = x_t' beta + sigma eps_t (2.2) with x_t = (D.z_t', y_t-1, z_t-1', D.x_t-1', ...,
D.x_t-k+1', 1)' [and t] and beta = (omega', alpha, -alpha kappa', gamma_1', ..., gamma_k-1', -alpha nu_c)'.
The z_t are generated by a VAR (3.2) with cointegrating rank 1, beta' = (1, -kappa'), and the weak exogeneity
restriction alpha = (alpha, 0)' (3.3); the ADL is the conditional equation of that VAR (Johansen 1992).
{cmd:ltscoint} builds x_t in exactly this order ({cmd:_ltsc_build}), so that the coefficient on L.y is
alpha and the coefficients on L.z_j are psi_j = -alpha kappa_j.


{marker lts}{...}
{title:The LTS estimator}

{pstd}
Given an index set zeta with #zeta = h, (2.4) defines beta_zeta = (sum_{t in zeta} x_t x_t')^-1
sum_{t in zeta} x_t D.y_t and sigma_zeta{c 94}2 = h{c 94}-1 sum_{t in zeta} (D.y_t - x_t' beta_zeta){c 94}2.
The LTS estimators are (2.5)

{p 12 12 2}zeta-hat = argmin_zeta sigma_zeta{c 94}2,   beta-hat = beta_{zeta-hat},   sigma-hat{c 94}2 = sigma{c 94}2_{zeta-hat}.{p_end}

{pstd}
Two algorithms are implemented in Mata.

{phang2}{bf:Exact enumeration} ({cmd:_ltsc_exact}, option {opt exact}): every set of T-h outliers is visited in
lexicographic order, OLS is computed on its complement and the smallest residual sum of squares is kept. This is
(2.5) literally and is used as the reference in small samples.{p_end}

{phang2}{bf:FAST-LTS} ({cmd:_ltsc_fastlts}; Rousseeuw and Van Driessen 2000, 2006). For s = 1,...,nsamp: draw a
random elemental subset of K observations (extended by further random observations if singular), fit OLS,
compute squared residuals for all T observations, take the h smallest as the new subset and refit
(a {it:concentration step}, {cmd:_ltsc_csteps}); each C-step cannot increase the h-trimmed sum of squares
(RvD, Theorem 1). After {it:csteps} steps keep the {it:nkeep} best subsets, iterate each to convergence
(the subset stops changing), and return the best. The full-sample OLS fit is always one of the starts, so h = T
returns OLS exactly. Because every C-step selects rows of the (D.y_t, x_t) design, no time-contiguity is imposed
on zeta, in line with Section 2.2 ("the indices in zeta need not be consecutive").{p_end}

{pstd}
Once zeta-hat is found, the command runs {helpb regress} of D.y on x_t and one impulse dummy per outlier over
all T observations. This regression has the same beta-hat and residuals on zeta-hat as OLS on the retained
observations, gives the oracle variance of Theorem 3 with T-K-(T-h) = h-K degrees of freedom, and its
log-likelihood is the one reported in the paper's (7.2). It also lets every official postestimation tool see
the right sample.


{marker theory}{...}
{title:Assumptions and the oracle property}

{pstd}
Good ADL errors are i.i.d. N(0,1) for t in zeta_T (4.1) and independent of the past (4.3); outlier errors satisfy
|eps_t| >= sqrt(2 log h) (4.2), i.e. they exceed the almost-sure maximum of h normal draws. Outliers may also
enter the VAR errors of z_t (Section 4.2), in G episodes of non-decreasing length, so that a good episode starts
from an initial value that aggregates the previous outliers (Section 3.2, Theorem 1).

{phang2}{bf:Theorem 2 (boundedness)} needs h_VAR > 2T/3 (4.8) and polynomially bounded outliers (4.6): up to a
third of the observations may be outlying. Remark B.1 shows 1/3 is binding in an autoregression because past
outliers propagate into the regressors. The command therefore refuses h <= 2T/3 unless {opt force}.{p_end}

{phang2}{bf:Theorem 3 (consistent selection and oracle expansion)} strengthens this to outliers of the same
logarithmic order as the largest good errors (4.9), a non-concentration condition on the regressors (4.11), and
h_VAR = T - o(sqrt(T)/log T) (4.12): the proportion of outliers vanishes although their number may diverge.
Then zeta-hat misclassifies a vanishing fraction, and sqrt(h)(sigma-hat{c 94}2 - sigma{c 94}2_zeta_T) and the
normalised beta-hat - beta_zeta_T are o_P(1). LTS inherits the asymptotic distribution of infeasible OLS on the
true good observations; Lemma B.7 shows the propagation of outliers into the regressors is negligible, and Lemma
B.8 gives normal limits for the I(0) directions and Dickey-Fuller-type limits for the I(1) directions.{p_end}

{pstd}
The command prints a note when T-h > sqrt(T), the regime where the paper's Tables 1-2 show inference starting
to deteriorate (most visibly for concentrated outliers and alpha = -0.2).


{marker map}{...}
{title:Step-to-equation map}

{col 5}{it:Code block}{col 34}{it:Paper}{col 60}{it:What it computes}
{col 5}{hline 90}
{col 5}_ltsc_build{col 34}(2.2), (2.3){col 60}x_t in paper order; D., L., L#D. operators
{col 5}_ltsc_fastlts / _ltsc_exact{col 34}(2.4), (2.5); RvD 2006{col 60}zeta-hat, beta-hat, RSS
{col 5}regress + dummies{col 34}Theorem 3; (7.2){col 60}oracle V, sigma-hat, log-likelihood
{col 5}sigma() option{col 34}(2.4) vs (7.2) vs suppl.{col 60}divisor h, h-K or T-K
{col 5}varsigma{c 94}2{col 34}ET (2026) eq. 2.4{col 60}SLTS variance = V / varsigma{c 94}4
{col 5}_ltsc_longrun{col 34}Sec. 6.1; suppl. p. 2{col 60}kappa = -psi/alpha, nu = -mu/alpha, D = dkappa/dtheta
{col 5}_ltsc_misspec{col 34}Sec. 7, (7.1)-(7.2){col 60}F_ar(1-2), F_arch(1), chi2_normal, F_hetero, T_h
{col 5}LR block{col 34}Sec. 5.3 (1), 6.3; HJNR 1998{col 60}T log(RSS_r/RSS_u); Doornik (2003) Tables 12-13
{col 5}_ltsc_ecmcv{col 34}Sec. 5.3 (1); E&M 2002{col 60}t_alpha with response-surface critical values
{col 5}_ltsc_wetest{col 34}Sec. 7 "Testing weak exogeneity"{col 60}D.z_j on (y - kappa'z)_t-1 + dummies
{col 5}_ltsc_profile{col 34}Table 7; BJN 2023 (6.1)-(6.2){col 60}T_h and IC_h over h
{col 5}_ltsc_dash{col 34}Figures 3-4{col 60}fit, scaled residuals, correlogram, QQ
{col 5}{hline 90}

{pstd}
{bf:Long-run block.} With theta = (omega, alpha, psi, mu) and kappa = -psi/alpha, D_kappa = (0, psi/alpha{c 94}2,
-1/alpha, 0)' and se(kappa-hat) = (D' V D){c 94}1/2 (supplement, p. 2); the same with mu for nu. Under weak
exogeneity, (kappa-hat - kappa)/se is asymptotically N(0,1) (Section 5.3, item 2; Johansen 1992).

{pstd}
{bf:Tests of no cointegration.} Under alpha = 0 both psi = -alpha kappa and the restricted deterministic term vanish
(Section 6.3), so the restricted regression drops y_t-1, z_t-1 and the constant (case H_c) or the trend (case
H_l, constant kept). LR = T log(RSS_r/RSS_u) is the likelihood-ratio statistic of Harbo et al. (1998) for a
partial system with p1 = 1 endogenous and p2 = dim z weakly exogenous variables; its limit is a functional of
Brownian motions whose quantiles, mean and variance Doornik (2003) tabulates from the Gamma approximation of
Doornik (1998). The p-value is 1 - Gamma(LR; shape = mean{c 94}2/var, scale = var/mean). For the paper's
example (H_l, p2 = 1) the 80% quantile is 11.18 (the paper quotes 11.0 from HJNR's own Table 2). The ECM
t-statistic t_alpha = alpha-hat/se is the PcGive unit-root test; Ericsson and MacKinnon (2002) give its
finite-sample quantiles q(T_a) = theta_inf + theta_1/T_a + theta_2/T_a{c 94}2 + theta_3/T_a{c 94}3 with
k = 1 + dim z and T_a = h - K (their Tables 3 and 4, constant and constant+trend).

{pstd}
{bf:Mis-specification tests} follow PcGive's definitions as used in (7.1)-(7.2): F_ar(1-2) is the F-form of the
Godfrey (1978) auxiliary regression of the residual on x_t, the dummies and two lagged residuals (missing lags set
to zero); F_arch(1) is the F-test of e_t{c 94}2 on e_t-1{c 94}2; chi{c 94}2_normal(2) is Doornik-Hansen (2008)
on the retained residuals; F_hetero regresses e_t{c 94}2 on the regressors and their squares (White 1980 without
cross-products). Degrees of freedom match the paper: e.g. F_ar(2,53) and F_hetero(12,49) for T = 67, K = 7,
T-h = 5. The paper notes ([64], [65]) that only some of these tests have theory covering outlier removal;
the truncation-corrected normality test of Berenguer-Rico and Nielsen (2026, Econometrics and Statistics) is
not implemented.


{marker inference}{...}
{title:Inference after LTS: oracle versus SLTS}

{pstd}
Theorem 3 implies that the variance of beta-hat is estimated by sigma-hat{c 94}2 (sum_{t in zeta-hat} x_t x_t')^-1
with no correction: this is {cmd:e(V)} by default. The robust-statistics convention (SLTS) instead assumes
uncontaminated normal errors and truncation, giving sigma-hat{c 94}2 -> sigma{c 94}2 varsigma{c 94}2 and an
efficiency factor, so that Var_SLTS = V / varsigma{c 94}4 with varsigma{c 94}2 = {lambda - 2 c phi(c)}/lambda,
c = Phi{c 94}-1((1+lambda)/2), lambda = h/T (Berenguer-Rico and Nielsen 2026, ET, eq. 2.4; the supplement's
SLTS rows). Both are stored; {opt slts} posts the latter. The paper's supplement shows SLTS is conservative
under the LTS model.


{marker h}{...}
{title:Determining h}

{pstd}
The paper takes h as given and, in Section 7, chooses it by inspecting the cumulant statistic over h (Table 7) and
by impulse-indicator saturation. {cmd:ltscoint hsearch} computes, for each h in [hmin, hmax], the LTS fit and

{p 12 12 2}T_h = h k3_h{c 94}2/6 + h k4_h{c 94}2/24,   k3_h, k4_h = skewness and excess kurtosis of the h retained residuals{p_end}
{p 12 12 2}IC_h = log(sigma_h{c 94}2) + 2 log log(T) (T-h)/T,{p_end}

{pstd}
(Berenguer-Rico, Johansen and Nielsen 2023, eqs 6.2 and 6.1). T_h is asymptotically chi{c 94}2(2) at the true h,
diverges at a logarithmic rate for h above it, and diverges for h below it because the retained residuals are then
truncated-normal; hence its minimiser estimates h. The published paper minimises T_h (h = 62 in their vintage);
the working-paper version used PcGive's chi{c 94}2_normal and a search restricted to h >= T - sqrt(T)/2 (h = 63).
The search default hmin = floor(2T/3)+1 respects boundedness. When {it:h} is not supplied, {cmd:ltscoint}
selects the T_h minimiser ({opt hsel(cumulant)}) or the IC_h minimiser ({opt hsel(ic)}) and reports the profile
in {cmd:e(hprofile)}. No theory for LTS inference with an estimated h is available (Section 8).


{marker decisions}{...}
{title:Decisions the paper leaves open, and what the command does}

{phang2}1. {bf:Divisor of sigma-hat{c 94}2.} (2.4) uses h; the supplement uses T-K; (7.2) implies h-K. Default
h-K ({opt sigma(df)}), matching the reported empirical results; the others are options. All are asymptotically
equivalent under (4.12).{p_end}
{phang2}2. {bf:Log-likelihood and LR.} Computed on the dummy-augmented regression over T observations, as in (7.2)
and the LR = 2(206.41-203.39) example; an h-based version would differ by the factor h/T inside the log.{p_end}
{phang2}3. {bf:Raw versus reweighted LTS.} The paper's estimator is the raw h-subset OLS; {cmd:robustbase::ltsReg}
returns reweighted coefficients by default. The command returns raw only.{p_end}
{phang2}4. {bf:No consistency factor on sigma.} Under the LTS model none is needed (BJN 2023); SLTS is offered as an
option, not applied silently.{p_end}
{phang2}5. {bf:Deterministic cases.} Restricted constant (H_c) or restricted trend with unrestricted constant (H_l)
only, as in (2.1) and (2.3).{p_end}
{phang2}6. {bf:Singular subsets.} Elemental subsets whose cross-product is singular are extended with further random
observations (RvD 2006); if still singular they are skipped.{p_end}
{phang2}7. {bf:Data shape.} Gaps in the time series produce missing lags and are dropped from the estimation sample;
the retained/outlier indicator refers to the estimation sample only.{p_end}


{marker validation}{...}
{title:Validation}

{phang2}- Exact enumeration and FAST-LTS return identical zeta-hat and RSS on all tested samples.{p_end}
{phang2}- On the UK data with T-h = 5, beta-hat equals the raw LTS of R's {cmd:robustbase::ltsReg} (alpha = 62/67,
nsamp = 5000) and Ben Jann's {cmd:robreg lts} (bp = 6) to all printed digits; all three select the outliers 2009,
2020, 2021, 2022, 2023, the set reported in the paper.{p_end}
{phang2}- Monte Carlo on the paper's DGP (6.1)-(6.2), Table 1 design (alpha = -0.2, kappa = 1, sqrt(T)/2 outliers in 5
episodes): with T = 400 and 400 replications the 5% t-test on kappa rejects 0.050 (LTS) and 0.005 (OLS), against
the paper's 0.054 and 0.004; all planted outliers are recovered. See {cmd:ltscoint_example.do}.{p_end}
{phang2}- The Gamma p-values return 0.050 and 0.200 at the tabulated 95% and 80% quantiles; the response surface
returns -3.84 for Ericsson and MacKinnon's worked example (k = 4, constant, T_a = 39).{p_end}


{title:Author}

{pstd}Dr Merwan Roudane{break}
merwanroudane920@gmail.com{break}
{browse "https://github.com/merwanroudane"}{p_end}