A model where the least trimmed squares estimator is maximum likelihood

Abstract

The least trimmed squares (LTS) estimator is a popular robust regression estimator. It finds a subsample of h ‘good’ observations among n observations and applies least squares on that subsample. We formulate a model in which this estimator is maximum likelihood. The model has ‘outliers’ of a new type, where the outlying observations are drawn from a distribution with values outside the realized range of h ‘good’, normal observations. The LTS estimator is found to be $h1/2$ consistent and asymptotically standard normal in the location-scale case. Consistent estimation of h is discussed. The model differs from the commonly used $ϵ$-contamination models and opens the door for statistical discussion on contamination schemes, new methodological developments on tests for contamination as well as inferences based on the estimated good data.

1 Introduction

The LTS estimator suggested by Rousseeuw (1984) is a popular robust regression estimator. It is defined as follows. Consider a sample with n observations, where some are thought to be ‘good’ and others are ‘outliers’. The user specifies that there are h ‘good’ observations. The LTS estimator finds the h sub-sample with the smallest residual sum of squares. Rousseeuw (1984) developed LTS in the tradition of Huber (1964) and Hampel (1971), who were instrumental in formalizing robust statistics. Huber suggested a framework of i.i.d. errors from an $ϵ$-contaminated normal distribution. Hampel formalized robustness and breakdown points. The present work, with its focus on maximum likelihood, deviates from this tradition.

Specifically, we propose a model in which the LTS estimator is maximum likelihood, a new approach to robust statistics. In this model, we first draw h ‘good’ regression errors from a normal distribution. Conditionally on these ‘good’ errors, we draw $n−h$ ‘outlier’ errors from a distribution with support outside the range of the drawn ‘good’ errors. The model is semi-parametric, so we apply a general notion of maximum likelihood. We provide an asymptotic theory for a location-scale model and find that the LTS estimator is $h1/2$-consistent and asymptotically standard normal. More than 50% contamination can be allowed under mild regularity conditions on the distribution function for the ‘outliers’. Interestingly, the associated scale estimator does not require a consistency factor. In practice, h is typically unknown. We suggest a consistent estimator for the proportion of ‘good’ observations, $h/n$⁠.

The approach of asking in which model the LTS estimator is maximum likelihood is similar to that taken by Gauss in 1809 (Hald, 2007, Section 5.5, 7.2). In the terminology of Fisher (1922), Gauss asked in which continuous i.i.d. location model is the arithmetic mean the maximum-likelihood estimator and found the answer to be the normal model. Maximum likelihood often brings a host of attractive features. The model provides interpretation and reveals the circumstances under which an estimator works well in terms of nice distributional results and optimality properties. Further, we have a framework for testing the goodness-of-fit, which leads to the possibility of first refuting and then improving a model.

To take advantage of these attractive features of the likelihood framework, we follow Gauss and suggest a model in which the LTS estimator is maximum likelihood. The model is distinctive in that the errors are not i.i.d. Rather, the h ‘good’ errors are i.i.d. normal, whereas the $n−h$ ‘outlier’ errors are i.i.d., conditionally on the ‘good’ errors, with distributions assigning zero probability to the realized range of the ‘good’ errors. When $h=n$⁠, we have a standard i.i.d. normal model, just as the LTS estimator reduces to the ordinary least squares (OLS) estimator. The model is semi-parametric, so we use an extension of traditional likelihoods, in which we compare pairs of probability measures (Kiefer & Wolfowitz, 1956) and consider probabilities of small hyper-cubes including the data point (Fisher, 1922; Scholz, 1980).

In practice, it is of considerable interest to develop a theory for inference for LTS. Within a framework of i.i.d. $ϵ$-contaminated errors, the asymptotic theory depends on the contamination distribution and the scale estimator requires a consistency correction (Butler, 1982; Čížek, 2005; Croux & Rousseeuw, 1992; Rousseeuw, 1984; Víšek, 2006). Since the contamination distribution is unknown in practice, inference is typically done using the asymptotic distribution of the LTS estimator derived as if there is no contamination. This seems fine for an infinitesimal deviation from the central normal model (Huber & Ronchetti, 2009, Section 12). However, these approaches would lead to invalid inference in case of stronger contamination—see Section 7.1. Within the present framework, the asymptotic theory is simpler. Specifically, we derive the asymptotic properties of the LTS estimator for a location-scale version of the presented model and find that the LTS estimator has the same asymptotic theory as the infeasible OLS estimator computed from the ‘good’ data, when it is known which data are ‘good’. As the asymptotic distribution does not depend on the contamination distribution, inference is much simpler.

In regression, a major concern is that the OLS estimator may be very sensitive to inclusion or omission of particular data points, referred to as ‘bad’ leverage points. LTS provides one of the first high-breakdown point estimators suggested for regression that also avoids the problem of ‘bad’ leverage. The presented model allows ‘bad’ leverage points.

Another key feature of the LTS model presented here is a separation of ‘good’ and ‘outlying’ errors, where the ‘outliers’ are placed outside the realized range of the ‘good’, normal observations. This way, the asymptotic error in estimating the set of ‘good’ observations is so small that it does not influence the asymptotic distribution of the LTS estimator. This throws light on the discussion regarding estimators’ ability to separate overlapping populations of ‘good’ and ‘outlying’ observations (Doornik, 2016; Riani et al., 2014).

LTS is widely used in its own right and often as a starting point for algorithms such as the MM estimator (Yohai, 1987) and the Forward Search (Atkinson et al., 2010). Many variants of LTS have been developed: nonlinear regression in time series (Čížek, 2005), algorithms for fraud detection (Rousseeuw et al., 2019) and sparse regression (Alfons et al., 2013).

In all cases, the practitioner must choose h, the number of ‘good’ observations. In our reading, this remains a major issue in robust statistics. We propose a consistent estimator for the proportion of ‘good’ observations, $h/n$⁠, in a location-scale model and study its finite sample properties by simulation. We apply it in the empirical example in Section 8.

The least median of squares (LMS) estimator, also suggested by Rousseeuw (1984), is closely related to LTS. The LMS estimator seemed numerically more attractive than LTS until the advent of the fast LTS approximation to LTS (Rousseeuw & van Driessen, 2000). Nonetheless, LMS remains of considerable interest. Replacing the normal distribution in the LTS model with a uniform distribution gives a model in which LMS is maximum likelihood. In the Online supplementary material, we show that the LMS estimator is h-consistent and asymptotically Laplace in the location-scale case. This is at odds with the slow $n1/3$ consistency rate found in the context of i.i.d. models.

We start by presenting the LTS estimator in Section 2. The LTS regression model is given in Section 3. The general maximum likelihood concept and its application to LTS are presented in Section 4 with details in Appendix  A. Section 5 presents an asymptotic theory for LTS in the location-scale case with proofs in Appendix  B. Estimation of h is discussed in Section 6. Monte Carlo simulations are presented in Section 7. An empirical illustration is given in Section 8.

In a supplement, the LMS estimator is analysed and detailed derivations of some identities in the LTS proof are given. All codes for simulations and empirical illustration are available from https://users.ox.ac.uk/˜nuff0078/Discuss and as supplementary files to the paper.

Notation: Vectors are column vectors. The transpose of a vector v is denoted $v′$⁠.

2 The LTS estimator

The available data are a scalar $yi$ and a p-vector of regressors $xi$ for $i=1,…,n$⁠. We consider a regression equation $yi=β′xi+σεi$ with regression parameter β and scale σ.

The least trimmed squares (LTS) estimator suggested by Rousseeuw (1984) is defined as follows. Given a value of β, the residuals are $ri(β)=yi−β′xi$⁠. The ordered squared residuals are denoted $r(1)2(β)≤⋯≤r(n)2(β)$⁠. The user chooses an integer $h≤n$⁠. Given that choice, the sum of the h smallest residual squares is obtained. Minimizing over β gives the LTS estimator

The LTS minimization classifies the observations as ‘good’ or ‘outliers’. The set of indices of the h ‘good’ observations is an h-subset of $(1,…,n)$⁠. We let $ζ^$ denote the set of indices i of the estimated ‘good’ observations satisfying $ri2(β^)≤r(h)2(β^)$⁠. Rousseeuw and van Driessen (2000) point out that $β^$ is a minimizer over least-squares estimators, that is $β^=β^ζ^$⁠, where

In the model proposed below, the maximum-likelihood estimator for the scale is the residual variance of the estimated ‘good’ observations

We show the consistency of $σ^2$ in the location-scale case. The problem of estimating the scale is intricately linked to the choice of model for the innovations $εi$⁠. Previously, Croux and Rousseeuw (1992) considered $σ^2$ in the context of a regression model with i.i.d. innovations with distribution function $F$⁠. For that model, they found that $σ^2$ should be divided by a consistency factor defined as the conditional variance of $εi$ given that $|F(εi)−1/2|≤h/(2n)$⁠. In practice, $F$ is unknown, so the normal distribution is often chosen for the consistency factor.

3 The LTS regression model

We formulate the LTS model where h ‘good’ errors are assumed i.i.d. normal, while $n−h$ ‘outliers’ are located outside the realized range of the ‘good’ errors. This differs from the $ϵ$-contaminated normal model of Huber (1964), where i.i.d. errors are normal with probability $1−ϵ$ and otherwise drawn from a contamination distribution. The supports of the normal and the contamination distributions overlap in the $ϵ$-contamination case, whereas ‘good’ and ‘outlying’ observations are separated in the below LTS model.

LTS regression model

 

The ‘outlier’ distributions $Gj$ are parameters in the LTS model. They are constrained to be continuous at zero to avoid overlap of ‘good’ and ‘outlier’ errors. Although the regressors are deterministic in Model 3.1, randomness of $x1,…,xn$ can be easily accommodated by a conditional model where the set of distributions $Gj(z)$ are replaced by the set of conditional distributions $Gj(z|xj)$⁠. The likelihood below then becomes a conditional likelihood.

The LTS model allows leverage points, since $Gj$ can vary with j. Informally, an observation is a ‘bad’ leverage point for the OLS estimator, if that estimator is very sensitive to inclusion or omission of that observation (Rousseeuw & Leroy, 1987, Section 1.1). An example is the star data shown in Figure 1 and used in the empirical illustration in Section 8.

A consequence of the model is that the ‘good’ errors must have consecutive order statistics. Randomly, according to the choice of $Gj$⁠, some ‘outlier’ errors are to the left of the ‘good’ errors. The count of left ‘outlier’ errors is the random variable

Thus, the ordered errors satisfy

The set ζ collects the indices of the observations corresponding to $ε(δ+1),…,ε(δ+h)$⁠. The difficulty in robust regression is of course that the errors are unknown. We note, that the extreme ‘good’ errors, $ε(δ+1)$ and $ε(δ+h)$⁠, are finite in any realization. As the ‘good’ errors are normal, the extremes grow at a $(2logh)1/2$ rate for large h, see Example 5.1 below. Likewise, the ‘outlier’ errors are also finite, but located outside the realized range from $ε(δ+1)$ to $ε(δ+h)$⁠.

4 Maximum-likelihood estimation in the LTS model

The LTS model is semi-parametric. We therefore start with a general maximum likelihood concept before proceeding to analysis of the LTS model.

4.1 A general definition of maximum likelihood

Traditional parametric maximum likelihood is defined in terms of densities, which are not well-defined here. Thus, we follow the generalization proposed by Scholz (1980), which has two ingredients. First, it uses pairwise comparison of measures, as suggested by Kiefer and Wolfowitz (1956), see Johansen (1978) and Gissibl et al. (2021) for applications. This way, a dominating measure is avoided. Second, it compares probabilities of small sets that include the data point, following the informality of Fisher (1922). This way, densities are not needed. Scholz’ approach is suited to the present situation, where the candidate maximum-likelihood estimator is known and we are searching for a model.

We consider data in $Rn$ and can therefore simplify the approach of Scholz. Let $P$ be a family of probability measures on the Borel sets of $Rn$⁠. Given a (data) point $z∈Rn$ and a distance $ϵ>0$ define the hypercube $Czϵ=(z1−ϵ,z1]×⋯×(zn−ϵ,zn]$⁠, which is a Borel set.

 

Following Scholz, define $P,Q$ to be equivalent at z and write $P=zQ$ if $P≤zQ$ and $Q≤zP$⁠. We get that (i) $P=zQ$ if and only if $limϵ→0{P(Czϵ)/Q(Czϵ)}$ exists and equals 1; (ii) $P<zQ$ and $Q<zR$ imply $P<zR$ (transitivity); and (iii) $P=zP$ for all $P∈P$ (reflexivity).

 

The maximum-likelihood estimator is unique if $P<zP^$ for all $P≠P^$⁠. As for Scholz, traditional maximum likelihood is a special case. Further, the empirical distribution function is maximum-likelihood estimator in a model of i.i.d. draws from an unknown distribution.

4.2 The LTS likelihood

To obtain the $ϵ$-likelihood for the LTS regression Model 3.1, we find the probability that the random n-vector of observations $yi$ belongs to an $ϵ$-cube $Czϵ$⁠. Throughout the argument, the regressors $xi$ are kept fixed. Conditioning ‘outliers’ on ‘good’ observations, we get

For the first product for ‘good’ observations, we note that $εi=(yi−β′xi)/σ$ is standard normal and define $ziβσ=(zi−β′xi)/σ$ and $ziβσϵ=ziβσ−ϵ/σ$⁠. Then, the factors of the first product in (4.1) are $Φ(ziβσ)−Φ(ziβσϵ)$⁠, which we denote $ΔϵΦ(ziβσ)$⁠.

For the second product, let $yiβσ=(yi−β′xi)/σ=εi$ and introduce

and a similar expression $z~jβσϵ$ for $zjβσϵ$⁠, noting that $z~jβσϵ=0$ if $zjβσϵ$ falls within the range from $mini∈ζεi$ to $maxi∈ζεi$⁠. A derivation in Appendix  A shows that the factors of the second product in (4.1) are $ΔϵGj(z~jβσ)=Gj(z~jβσ)−Gj(z~jβσϵ)$⁠. With this notation, the probability (4.1) of the $ϵ$-cube is

The $ϵ$-likelihood arises from the probability (4.2) evaluated in the data point. That is, we replace the vector z by the data vector y. We define $y~jβσ$ and $ΔϵGj(y~jβσ)$ in a similar fashion as before. Thus, we get the $ϵ$-likelihood

We note that the ‘good’ errors cannot be repeated due to the assumptions that the ‘good’ observations are continuous and that the centred ‘outliers’ are continuous at zero. Thus, parameter values resulting in repetitions of the ‘good’ errors are to be ignored. As an example, suppose we have $n=10$ observations and we have chosen a β so that the residual $yi−β′xi$ takes the ordered values: 1,1,1,2,3,6,6,7,8,9. The values 1 and 6 are repetitions and cannot be ‘good’. Thus, for $h=3$⁠, we can only select ζ as the index triplet 7,8,9. All other choices are asigned a zero likelihood. The issue is essentially the same as in OLS theory, where the least-squares estimator may be useful when there are repetitions, but it cannot be maximum likelihood, as repetitions happen with probability zero under normality.

4.3 Maximum-likelihood estimation

The two products in the $ϵ$-likelihood (4.3) resemble a standard normal likelihood and a product of $n−h$ likelihoods, each with one observations from an arbitrary distribution. We will exploit those examples using profile likelihood arguments.

First, suppose $β,σ,ζ$ are given. Then, the first product in the LTS $ϵ$-likelihood (4.3) is constant. In the second product, the jth factor only depends on $Gj$⁠. Since the $Gj$ functions vary in a product space, we maximize each factor separately. The maximum value for $ΔϵGj(yjβσ)$ is unity for any $ϵ>0$ and regardless of the value of $β,σ,ζ$⁠. The maximum is attained for any distribution function that is continuous at zero and that rises from zero to unity over the interval from $yjβσϵ$ to $yjβσ$⁠. This set of distribution functions includes the distribution with a point mass at $y~jβσ$⁠, which satisfies the requirement of continuity at zero, because the ‘outliers’ are separated from the ‘good’ observations. Moreover, the point mass distribution is unique in the limit for vanishing $ϵ$⁠. Thus, maximizing (4.3) over $Gj$ gives the profile $ϵ$-likelihood

Second, suppose ζ is given. We argue that the unique maximum-likelihood estimator in the sense of Definition 4.2 is given by

We need to show that $limsupϵ→0{LGϵ(β,σ,ζ)/LGϵ(β^ζ,σ^ζ,ζ)}<1$ for all $(β,σ)≠(β^ζ,σ^ζ)$⁠. Note that $ϵ−hLGϵ(β,σ,ζ)$ converges to a standard likelihood for a regression with normal errors. Therefore, the condition is satisfied as long as $∑i∈ζxixi′$ is invertible. Thus, maximizing (4.4) over $(β,σ)$ gives a profile $ϵ$-likelihood for ζ satisfying, for $ϵ→0$⁠,

Third, the profile likelihood is maximized by choosing ζ so that $σ^ζ$ is as small as possible. We note that the maximization is subject to the constraint that none of the ‘good’ residuals are repeated. Apart from the latter constraint, this is the LTS estimator described in (2.2). We summarize.

 

5 Asymptotic theory for the location-scale model

We present an asymptotic theory of the LTS estimator in the special case of a location-scale model $yi=μ+σεi$⁠. In this model, the observations $yi$ and the unknown errors $εi$ have the same ranks. Thus, the ordering in (3.3) is observable and we only need to consider values of ζ corresponding to indices of observations of the form $y(δ+1),…,y(δ+h)$ for some $δ=0,…,n−h$⁠. Following Rousseeuw and Leroy (1987, p. 171), the LTS estimators then reduce to $μ^=μ^δ^$ and $σ^=σ^δ^$⁠, where

5.1 Sequence of data generating processes

We consider a sequence of LTS models indexed by n, so that $hn→∞$ as $n→∞$⁠. In this sequence, the ‘outliers’ have a common distribution $G$⁠, which is continuous at zero. If $hn=n$ the LTS estimator reduces to the full sample least-squares estimator with standard asymptotic theory. Here, we choose

where λ is the asymptotic proportion of ‘good’ observations. The parameters $ζn$⁠, $δn$ vary with n, while $μ,σ,G$ are constant in n.

We reparametrize $G$ in terms of separate distributions for left and right ‘outliers’. Let

Thus, as $ξj$ is $G$-distributed, we have that $ε_j=−ξj$ is $G_$-distributed when $ξj<0$ and $ε¯j=ξj$ is $G¯$-distributed when $ξj>0$⁠. This leads to

This way, the ‘outliers’, $εj$ for $j∉ζn$⁠, can be constructed through a binomial experiment as follows. Draw $n−h$ independent $Bernoulli(ρ)$ variables. If the jth variable is unity then $εj=mini∈ζεi−ε_j$⁠. If it is zero then $εj=maxi∈ζεi+ε¯j$⁠. Hence, the number of left ‘outliers’ is

by the Law of Large Numbers. In summary, the sequence of data generating processes is defined by $μ,σ,ρ,G_,G¯$⁠, and $hn,δn$⁠.

5.2 The main result

We show that the asymptotic distribution of the LTS estimator is the same as it would have been, if it had been known which observations were ‘good’. Here, we consider the case where the ‘good’ observations are normal and there are more ‘good’ observations than ‘outliers’. The result does not require any regularity conditions for the common ‘outlier’ distribution $G$⁠.

 

We expect that Theorem 5.1 would generalize to the LTS regression Model 3.1. That is, the LTS regression estimator will have the same asymptotic theory as if we knew which observations were ‘good’. The proof of such a result is, however, an open problem.

Theorem 5.1 for the LTS location-scale model differs from the known results for i.i.d. models. Butler (1982) proved a general result for the model with i.i.d. errors, showing that $n1/2(μ^−μ)$ is asymptotically normal when the errors are symmetric with a strongly unimodal distribution. Further, the asymptotic variance involves an efficiency factor that differs from unity even in the normal case. He also gives an asymptotic theory for nonsymmetric errors. Further, Víšek (2006) analysed the regression estimator for i.i.d. errors. Johansen and Nielsen (2016b, Theorem 5) provide an asymptotic expansion of the scale estimator under normality.

A corresponding result for the LMS estimator is given in the supplement. The LMS estimator is maximum likelihood in a model where the ‘good’ observations are uniform. When there are no outliers, $n=h$⁠, the LMS estimator is a Chebychev estimator. Knight (2017) provides asymptotic analysis of Chebychev estimators in regression models. Coinciding with that theory, Online Supplementary Material, Theorem C.2 shows that the LMS estimator is super-consistent at an $n−1$ rate in that model, which contrasts with the well-known $n−1/3$ rate for i.i.d. models. However, Online Supplementary Material, Theorem C.2 does require regularity conditions for the ‘outliers’ to give sufficient separation from the ‘good’ observations.

5.3 Extentions of the main result

5.3.1 Nonnormal errors

We start by relaxing the normality assumption. In the spirit of OLS asymptotics, we present sufficient conditions for the ‘good’ errors: $2+$ moments, symmetry, and exponential tails.

    

5.3.2 More ‘outliers’ than ‘good’ observations

In this section, we allow for more ‘outliers’ than ‘good’ observations. Although in the traditional breakdown point analysis there has to be more ‘good’ observations than ‘outliers’ (Rousseeuw & Leroy, 1987, Section 3.4), in practice, LTS estimators are sometimes used with more ‘outliers’ than ‘good’ observations, as in the Forward Search algorithm (Atkinson et al., 2010) for instance. We show that, within the LTS model framework, it is possible to allow for more ‘outliers’ than ‘good’ observations, but in this case, we need some regularity conditions on the ‘outliers’ to make sure the ‘good’ observations can be found.

Let the proportion of ‘good’ observations satisfy $0<λ<1$⁠. Recall, that $δn/(n−hn)→ρ=G(0)$$a.s$⁠. is the proportion of the ‘outliers’ that are to the left. The number of ‘outliers’ to the right is $n¯=n−hn−δn$⁠. Define also the proportion of left and right ‘outliers’ relative to the number of ‘good’ observations through

Regularity conditions are needed for the ‘outlier’ distribution when $ω_≥1$ or $ω¯≥1$⁠. Note that $ω_=ω¯<1$ when the proportion of ‘good’ observations is $λ>1/3$ and proportions of left and right ‘outliers’ are the same, so that $ρ=1/2$⁠. The following definition is convenient for analysing the empirical distribution function of the ‘outliers’, evaluated at a random quantile.

  

The following assumption requires that the ‘outlier’ distribution is regular and that the ‘outliers’ are more spread out than the ‘good’ observations.

  

Given the novelty (and perhaps at first surprising content) of Theorem 5.3, a discussion of this result and its relationship with Theorem 5.1 seems in order. Recall that Theorem 5.1 has no regularity conditions on the ‘outliers’. The result exploits that there are more than $50%$ ‘good’ observations and the thin normal tails of these ‘good’ observations separate them sufficiently from the ‘outliers’. In Theorem 5.3, we allow less than $50%$ ‘good’ observations but require regularity conditions on the outliers to ensure sufficient separation. More specifically, the proof of Theorem 5.3 covers two situations. First, if the proportions of left and right ‘outliers’ are equal and more than one-third of the observations are ‘good’, then $ω¯=ω_<1$⁠, so that Assumption 5.2 is nonbinding. Second, in general, Assumption 5.2(i,ii) implies that the ‘outliers’ are spread out more than the ‘good’ observations.

This has a potential parallel with the breakdown point analysis of (Rousseeuw & Leroy, 1987, Section 3.4), which shows that for a given dataset the distribution of the LTS estimator is bounded with less than $50%$ contamination of an arbitrary type. It may still be possible to obtain the same result allowing for more than $50%$ contamination that satisfies regularity conditions of the type considered in Theorem 5.3, but we leave this as an open question.

5.4 The OLS estimator in the LTS model

We show that the least-squares estimator $μ¯=n−1∑i=1nyi$ can diverge in the sequence of LTS models. This implies that the least-squares estimator is not robust within the LTS model in the sense of Hampel (1971). We assume all ‘outliers’ are to the right, so that $ρ=0$⁠.

 

6 Estimating the proportion of ‘good’ observations

The LTS regression model takes the number of good observations, h, as given. In practice, an investigator has to choose h. Estimation of h is a difficult problem, and methods for estimating h are scarce in the literature. In this section, we start by reviewing a commonly used method for choosing h and discuss its asymptotic properties. We will then propose some new methods for consistently estimating the rate $λ=h/n$ of ‘good’ observations. The best of these methods is consistent at a $logh$ rate. This is a slow rate and the method may be imprecise about the number h of ‘good’ observations even in large samples. Notwithstanding, to the best of our knowledge, this is the only consistent method available in the literature. Further research and improvements in this area are needed.

6.1 The index plot method

A common method for estimating h is to apply a high breakdown point LTS estimator selecting approximately $n/2$ observations, then compute scaled residuals for all observations and keep those observations for which the scaled residuals are less than 2.5 in absolute value. This is conveniently done using an index plot (Rousseeuw & Hubert, 1997; Rousseeuw & Leroy, 1987). More specifically, scaled residuals $ε^iς$ are plotted against the index i. Here, the consistency factor $ς=0.615$ is the conditional standard deviation of $εi$ given that $1/4<Φ(εi)<3/4$ (Croux & Rousseeuw, 1992). Bands on $±2.5$ are displayed on the plot and observations corresponding to scaled residuals $ε^iς$ beyond these bands are declared outliers. Hence, the estimator of the number of outliers is $n−h^IP=∑i=1n1(|ε^iς|>2.5).$ We note that when scaling the residuals, the consistency factor is based on the assumption that all errors are i.i.d. normal.

We can analyse the asymptotic properties of the index plot method using Theorems 4, 5, 7 in Johansen and Nielsen (2016b). This will show, that in a clean, normal sample, the index plot method, asymptotically, finds $γ^IP=(n−h^IP)/n→Pγ=P(|εi|>2.5)=2Φ(−2.5)=1.24%$ outliers. In the terminology of Castle et al. (2011); Hendry and Santos (2010), γ is the asymptotic gauge of the procedure. Thus, the index plot method will on average estimate a nonzero proportion of outliers, even in uncontaminated samples.

6.2 Methods based on the LTS model

We consider methods for consistent estimation of the proportion of ‘good’ observations in the LTS model. This is a semi-parametric model selection problem where the dimension of the parameter space increases with decreasing h, since the number of $Gj$ functions increases when the number of ‘outliers’ increases. We consider three estimation methods. In all cases, suppose $h≥h_$ where $h_>dimxi$ is a lower bound chosen by the user.

Maximum likelihood. We argue that $h^=h_$ is the maximum-likelihood estimator for h. Let $Lϵ(β,σ,ζ,Gj,h)$ be the $ϵ$-likelihood given in (4.3) for each h. Let $σ^h$ denote the maximum-likelihood estimator in Theorem 4.1 for each h. By (4.6), the profile likelihood for h is

We note that $σ^h2>0$ for all $h≥h_$⁠, so that $ϵh_−hL^hϵ/L^h_ϵ$ is bounded uniformly in h. Thus, $L^hϵ/L^h_ϵ→0$ as $ϵ→0$ whenever $h>h_$⁠, so that $h^=h_$⁠.

An information criteria. Penalizing $logσ^h2$ gives the information criteria

for a penalty f. Let $h^IC$ be the minimizer. Then $λ^IC=h^IC/n$ estimates the proportion of ‘good’ observations, $λ∘$ say. Below, we argue, for the location-scale case, that $λ^IC$ is consistent if the penalty is chosen so that $f(n)→∞$ for increasing n, but $f(n)<λ∘−1loglogn$⁠, where $0<λ∘≤1$⁠. If, for instance, it is expected that more than half of the observations are ‘good’, so that $λ∘>1/2$⁠, the penalty can be chosen as $f(n)=2loglogn$⁠. The intuition is as follows.

Consider data generating processes as in Section 5.1 and let $hn∘$ denote the number of ‘good’ observations, so that $hn∘/n→λ∘$⁠. We want to show that $λ^IC=λ∘$ in the limit. Theorem 5.1 shows that $σ^hn∘2$ is consistent for $σ2$ when $h=hn∘$⁠. Thus, we consider sequences $hn$ so that $hn/n→λ≠λ∘$ and argue that $IChn−IChn∘∼log(σ^hn2/σ2)+f(n)(λ∘−λ)$ has a positive limit, so that the minimizer $λ^IC=λ∘$ in the limit.

Suppose $λ<λ∘$⁠. In that case, one could expect that, asymptotically, the estimated set of ‘good’ observations is a subset of the good observations, so that $ζ^h⊂ζ$⁠. Further, $σ^h2$ converges to $σλ2$⁠, say, the variance of a truncated normal distribution as analysed by Butler (1982). Then, $IChn−IChn∘∼log(σλ2/σ2)+f(n)(λ∘−λ)$⁠, which is positive in the limit when $f(n)$ diverges. Thus, $λ^IC≥λ∘$ in the limit.

Suppose $λ>λ∘$⁠. Then $IChn−IChn∘∼log(σ^hn2/σ2)−f(n)(λ−λ∘)$⁠, where the penalty term is negative, so that $f(n)$ must be chosen carefully. With the normal LTS model, $σ^hn2$ must diverge. The reason is that, for $λ>λ∘$⁠, the estimated set of ‘good’ observations $ζ^h$ includes both ‘good’ observations and ‘outliers’. The ‘outliers’ diverge at a $(2loghn∘)1/2$ rate, see Example 5.1, so that $σ^hn2$ must diverge at a $2loghn∘$ rate. Noting that $hn∘∼nλ∘$⁠, we get $IChn−IChn∘∼loglogn−f(n)(λ−λ∘)$⁠. When $1/2<λ∘<λ≤1$⁠, we have $λ−λ∘<1/2$⁠, and the $loglogn$ term dominates when $f(n)≤2loglogn$⁠. Thus, $λ^IC≤λ∘$ in the limit.

Combining these arguments indicates that $λ^IC=λ∘$ in the limit. Unfortunately, in simulations that are not reported here, we find that a $2loglogn$ penalty grows so slowly in n that for some specifications the consistency is only realized for extremely large samples.

Cumulant based normality test. A more useful estimator for the proportion of ‘good’ observations λ is the minimizer $h^T$⁠, say, of the normality test statistic

based on third and fourth cumulants of the estimated ‘good’ residuals. We argue that $h^T/n$ is consistent for $λ∘$⁠. This is supported by a simulation study in Section 7.2.

The intuition of the consistency argument is similar to that for the information criteria. First, for $h=hn∘$⁠, Theorem 5.1 may be extended to show that $Th0$ is asymptotically $χ22$⁠. For $h<hn∘$⁠, the sample moments may converge to the moments of a truncated normal distribution, see Berenguer-Rico and Nielsen (2017). Thus, $Th$ would diverge as it is normalized using the normal distribution. For $h>hn∘$⁠, the estimated set of ‘good’ observations $ζ^h$ contains both ‘good’ and ‘outlying’ observations, so that $Th$ diverges at a logarithmic rate, instead of the iterated logarithmic rate for the information criteria. A formal proof is left for future work.

7 Simulations

7.1 Inference in the location-scale model

In the location-scale model $yi=μ+σεi$⁠, we study the finite sample properties of tests for $μ=0$ using four different tests statistics and six different data generating processes, which we describe below. We consider sample sizes $n=25,100,400,1600$ and let $h/n=λ=0.8$⁠.

Tests. We consider four different estimators: LTS, LMS, OLS and SLTS (scale-corrected LTS estimator). For each estimator, s say, we compute t-type statistics, $ts=(μ^s−μ)/ses$⁠, with associated asymptotic 95% quantiles $qs$⁠.

For the OLS estimator, we have $μ^OLS=n−1∑i=1nyi$ and $seOLS2=σ^OLS2/n$ with $σ^OLS2=n−1∑i=1n(yi−y¯)2$⁠. The asymptotic quantile $qOLS$ is standard normal.

For the LTS estimator, we apply the estimators $μ^LTS$ and $σ^LTS$ given in (5.1). By Theorem 5.2, we get that $seLTS2=σ^LTS2/h$⁠. The asymptotic quantile $qLTS$ is standard normal.

For the LMS estimator, we apply the estimators $μ^LMS$ and $σ^LMS$ given in Online Supplementary Material, (C.4). Online Supplementary Material, Theorem C.2 gives that $seLMS2=σ^LMS2/h2$⁠. The asymptotic quantile $qLMS$ is standard Laplace.

Finally, SLTS uses the usual approach to LTS estimation with consistency and efficiency correction factors arising from truncation in a standard normal distribution as outlined in the end of Section 2. Let $ςh/n2=∫−ccx2φ(x)dx/∫−ccφ(x)dx$ with c chosen so that $∫−ccφ(x)dx=h/n$⁠, which gives $ς0.82=0.438$⁠. Then, we have estimators $μ^SLTS=μ^LTS$ and $σ^SLTS2=σ^LTS2/ςh/n2$⁠, while $seSLTS2=σ^SLTS2/(nςh/n2)$⁠. The asymptotic quantile $qSLTS$ is standard normal.

Data generating processes (DGPs). Table 1 gives an overview of the DGPs. The first three DGPs are examples of the LTS Model 3.1 and the LMS Online Supplementary Material, Model C.1. The proportion of ‘good’ observations is $λ=80%$⁠. The ‘good’ errors $εi$ are i.i.d. normal $N(0,1)$ in DGP1 and DGP2 and i.i.d. uniform $U[−1,1]$ in DGP3. The ‘outlier’ errors are defined as in (3.1), so that $εj=(maxi∈ζεi+ξj)1(ξj>0)+(mini∈ζεi+ξj)1(ξj<0)$⁠, where $ξj−ν+1(ξj>0)+ν−1(ξj<0)$ are i.i.d. normal $N(0,1)$ and $ν+$ and $ν−$ separate ‘good’ and ‘outlying’ observations. The separators are $ν+=ν−=0$ in DGP1 and $ν+=3$⁠, $ν−=1$ in DGP2 and DGP3.

The last three DGPs are examples of $ϵ$-contamination (Huber, 1964). We draw n observations from the distribution function $0.8Φ+0.2H$⁠, where Φ is standard normal and $H$ represents contamination. In DGP4, $H=Φ$⁠, giving a standard i.i.d. normal model. In DGP5 and DGP6, $H$ is $N(0,3)$ and $N(2,1)$⁠, giving symmetric and nonsymmetric mixtures, respectively.

We have different maximum-likelihood estimators for the different models. These are LTS for DGP1 and DGP2, LMS for DGP3, and OLS for DGP4. None of the considered estimators are maximum likelihood for DGP5 and DGP6.

Table 2 reports results from $106$ repetitions. The Monte Carlo standard error is 0.001.

The OLS statistic is maximum likelihood with DGP4 and performs well. It performs equally well with the symmetric, i.i.d. DGP5. For DGP1, it is slowly diverging, possibly because the absolute sample mean is diverging with n. OLS performs poorly with the nonsymmetric DGP2, DGP3, and DGP6.

The LTS statistic is maximum likelihood with DGP1 and DGP2. The asymptotic theory also applies for DGP3. The convergence is slow for DGP1, where there is no separation. The LTS statistic does not perform well with $ϵ$-contamination in DGP4, DGP5, and DGP6.

The LMS statistic is maximum likelihood with DGP3 and perform well with that DGP, but poorly with all other DGPs. The SLTS statistic is not maximum likelihood for any of the considered models. It is calibrated to be asymptotically unbiased for DGP4 and performs well for that model, but poorly with all other DGPs.

Overall, we see that it is a good idea to apply maximum likelihood but this does require that the model specification is checked. In particular, the LTS estimator is best in DGP1–DGP3, although with some finite sample distortion with DGP1 where ‘good’ and ‘outlying’ observations are not well-separated. The LTS estimator does not work well for the $ϵ$-contaminated models, where a model dependent scale correction is needed. The usual approach of using the normal scale correction as in SLTS does not work well in general. All estimators are poor for asymmetric $ϵ$-contamination.

7.2 h estimation

Next, we study the finite sample properties of estimating h using the cumulant-based normality test statistic $Th$ in (6.2). Results are reported in Table 3, based on $103$ repetitions.

In each repetition, we compute the $Th$ statistic for each h in the range from 60% to 90% of n. The estimator of h is the minimizer of $Th$ over that range.

The data generating processes are DGP1 and DGP2 from above. These are examples of the LTS model, so that DGP1 has symmetric ‘outliers’ that are not separated from the ‘good’ observations, while DGP2 has asymmetric ‘outliers’ that are separated from the ‘good’ observations. For each DGP, we consider cases with 70% or 80% ‘good’ observations.

Table 3 reports three quantities for each of the DGPs. First, $h^bias$ and $h^sd$ are the Monte Carlo average and standard deviation of the estimation error $h^−h$⁠. Further, $h^r$ is a binary variable, which is checked if $h^$ is in the interior of the range from 60% to 90% of n for all $103$ simulations. The theory suggests that the proportion $h/n$ of ‘good’ observations is consistently estimated, whereas h is not consistently estimated. Thus, we would expect $h^bias$ and $h^sd$ to grow slower than linearly in n, but not to vanish.

For all the four set-ups, the simulations confirm that $h^/n$ is consistent for λ. In DGP1, where there is only little separation between ‘good’ and ‘outlying’ observations, the performance differs substantially between the cases $λ=0.7$ and $λ=0.8$⁠. We do not have an explanation for this difference. Nonetheless, the estimation works much better for DGP2, with its separation between ‘good’ observations and ‘outliers’, in both cases $λ=0.7$ and $λ=0.8$⁠, matching the theoretical discussion in the above sections.

We also considered the information criteria $Φh$ in (6.1), but omit the results. The performance of $Φh$ is poor, quite possibly due to the $loglogn$ rates involved.

8 Empirical illustrations

We provide two empirical illustrations using the stars data of Rousseeuw and Leroy (1987) and the state infant mortality data of Wooldridge (2015). Both analyses illustrate the estimation of h. The second case also illustrates inference in the LTS model. Therefore, we must determine empirically if the LTS model is appropriate. Indeed, there could be a contamination pattern that is consistent with the $ϵ$-contamination, with the LTS model, or with neither of those. For both illustrations, we study the source of the data to arrive at reasonable empirical models.

Throughout, we use R version 4.0.2 (R Core Team, 2020) estimating LTS using ltsReg from the Robustbase package. Before each LTS call we apply set.seed(0). R code with embedded data is available as Supplementary material.

8.1 The stars data

For this empirical illustration, we consider the data on log light intensity and log temperature for the Hertzsprung–Russell diagram of the star cluster CYG OB1 containing $n=47$ stars as reported by Rousseeuw and Leroy (1987, Table 2.3). Figure 1 shows a cross plot of the variables, where the log temperature axis is reversed. The majority of observations follow a steep band called the main sequence. Rousseeuw and Leroy (1987) refer to the four stars to the top right of Figure 1 as ‘outliers’.

By consulting the original source of the data in Humphreys (1978), see also (Vansina & De Grève, 1982, Appendix  A), we found that the 4 stars to the right of Figure 1 are of M-type (observations 11, 20, 30, 34) and they are red supergiants. Further, the fifth star from the right is of F-type (observation 7, called 44 Cyg). The next 31 stars (1 doublet) from the right are of B-type and the remaining 11 stars (1 doublet) furthest to the left are of the O-type. The doublets are not exact doublets in Humphreys’ original data, so this in itself should not be seen as evidence against a normality assumption.

We fitted the linear model $loglight=β1+β2logtemperature+σεi$⁠. Table 4, left panel, shows LTS estimates for different h values. For $h=n=47$⁠, the LTS estimator is the full sample OLS estimator. The β estimates are the ‘raw coefficients’ found by ltsReg while $σ^$ is computed directly from (2.3) without any consistency correction. Figure 1 shows LTS fits for selected values of h. It is seen that the fits rotate when h increases. Table 4 also reports the $Th$ criterion as a function of h. It is minimized for $h=42$ pointing at five ‘outliers’: The four M-stars and the F-star. Figure 1 indicates estimated ‘good’ observations and ‘outliers’ with solid and open dots.

The estimation of h by the statistic $Th$ appears a little shaky in Table 4, left panel. The lowest value is obtained for $h=42$⁠, while the value for $h=44$ is nearly as low. The slope coefficients $β^2$ change gradually for $h>42$ with no obvious choice of h for $h>42$⁠. Table 4, right panel, shows corresponding results when dropping the F star from the sample. The results are the same for $h≤42$⁠. We see that $h=42$ is now a clear minimum identifying the four M-type stars as ‘outliers’. The M stars appear to have a masking effect, where after their deletion, the F star emerges as very influential in the sense of Rousseeuw and Leroy (1987, p. 81). Perhaps for this reason, different conclusions are reached by different traditional methods. An LTS index plot points at 5 ‘outliers’, see Section 6.1, while an MM index plot (using lmrob in the R package robustbase) and LMS residuals (using lmsreg in the R package MASS) both point at 4 ‘outliers’, not detecting the F star as ‘outlying’.

Figure 2 shows kernel density plots for the scaled residuals for $h=25,42,47$⁠. The black, thin lines gives kernel densities for the full sample. The red, dashed lines gives kernel densities for the estimated ‘good’ observations. The standard normal distribution is shown with a blue, thick line. For $h=42$⁠, the red, blue, and part of the black lines coincide, which indicates the normality of the ‘good’ observations. The full sample kernel density has a probability mass in the right tail corresponding to the four giants. There is a slight discrepancy between the full sample and the ‘good’ kernel densities in the region from 2 to 4 corresponding to the F star. By construction this 43rd residual will be outside the range of the 42 ‘good’ residuals, but not by far. Kernel densities are very sensitive to the choice of kernel and bandwidth. For illustration, we chose a Gaussian kernel and bandwidth $1.5h−1/5$ to get the best match of the red and blue curves for $h=42$⁠. With that choice, we can more clearly see discrepancies between the kernel density for the ‘good’ observations and the normal curve for $h=25,47$⁠, that is, for the LTS with a high breakdown point and the OLS estimator, respectively.

8.2 State infant mortality rates

We consider 1990 data for the United States on infant mortality rates by state, including Washington DC, which has a particularly high infant mortality rate (Wooldridge, 2015, p. 299). The data are available as infmrt in the R package wooldridge and are taken from U.S. Bureau of the Census (1994). We analyse two models.

The first model follows Wooldridge. It is a linear regression of the number of deaths within the first year per 1,000 live births, infmort, on the log of per capita income, lpcinc, the log of physicians per 100,000 members of the civilian population, lphysic, and the log of population in thousands, lpopul. In Figure 3, the graph using circle symbols shows the cumulant-based criteria $Th$ as function of h. The OLS regression is obtained for $h=n=51$ and has a rather large value of $Th$—notice the gap in the $Th$-axis—indicating that the full sample model is mis-specified. Choosing $h=50$ would lead to one ‘outlier’, which is Washington DC. However, the $Th$ function is minimized at $h=45$ indicating six ‘outliers’: Delaware, Washington DC, Georgia, Texas, California, and Alaska (U.S. Bureau of the Census, 1994, Table 123). The interpretation is not obvious and could be due to a missing regressor.

The second model differes in two respects. It applies logarithms to the regressand to stabilize rates close to zero as well as for Washington DC. It also includes a regressor for the log proportion of black people in the population, lblack (U.S. Bureau of the Census, 1992, Table 255), since infant mortality is quite different for white and black infants in most states (U.S. Bureau of the Census, 1994, Table 123). In Figure 3, the graph using square symbols shows $Th$ versus h for this model. The minimizer is $h=50$ with Washington DC as the ‘outlier’. The minimum of $0.08$ is small compared with a $χ22$ distribution, so no evidence against the LTS model. We note that the $Th$ function is also quite low for h in the left side of the plot, albeit not as low as for $h=50$⁠. The estimated LTS model for $h=50$ is

The standard errors are the usual OLS standard errors as the LTS model appears to apply. We conclude that the variables lpcinc and lpopul are not significant.

Changing the regressand to be measured on the original scale introduces a second ‘outlier’, South Dakota, which has one of the lowest infant mortality rates for the black population.

9 Concluding remarks

Estimating the proportion of ‘good’ observations. We consider a few issues regarding the estimation of h and the adequacy of the model in the two empirical examples.

First, the sample sizes, $n=47$ and $n=50$⁠, are rather small. The results in Section 6.2 indicate that the estimation of $λ=h/n$ based on the $Th$ criterion is $logh$ consistent, hence, it may be imprecise in small samples.

Second, the results in Section 6.2 focus on the location-scale LTS model. The empirical illustrations above consider models with regressors. We believe the estimation of λ based on the $Th$ criterion extends to a well-specified multiple LTS regression model. However, as pointed out in Section 8.2, the omission of regressors could interfere with the estimation of λ. This aspect needs further investigation.

Third, the results in Section 6.2 refer to the LTS model, where, once ‘outliers’ have been removed, the remaining ‘good’ observations look normal. In practice, the adequacy of the model in a given data set has to be studied. A formal testing procedure is not yet available for the present model. Estimating λ by the $Th$ criterion is unlikely to be consistent in a model where ‘outliers’ are generated by $ϵ$-contamination.

Fourth, notwithstanding the above issues, the suggested procedure for selecting h seems to work well in the two empirical illustrations in this section and helped in finding satisfactory LTS models for the data.

Other models of the LTS type. New models and estimators can be generated by replacing the normal assumption in the LTS models with some other distribution. For instance, the uniform leads to the least median of squares (LMS) estimator analysed in Online Supplementary Material, Appendix C, while the Laplace distribution leads to the least trimmed sum of absolute deviations (LTA) estimator (Dodge & Jurečková, 2000; Hawkins & Olive, 1999; Hössjer, 1994, Section 2.7). Ron Butler has suggested to us that the approach in this paper could also be applied to the minimum covariance determinant (MCD) estimator for multivariate location and scale (Butler et al., 1993; Rousseeuw, 1985).

Alternative models for ‘outliers’. The maximum likelihood argument in Section 4 would also work if the ‘outlier’ errors $εj$ rather than $ξj$ have distribution $Gj$⁠. The analysis would be related to the trimmed likelihood argument of Gallegos and Ritter (2009). However, the resulting LTS estimator would have less attractive asymptotic properties in that model compared to those derived in this paper.

Inference requires a model for both ‘good’ and ‘outlying’ observations. In the presented theory, the ‘good’ and the ‘outlying’ observations are separated. The traditional approach, as advocated by Huber (1964), is to consider mixture distributions formed by mixing a reference distribution with a contamination distribution. Any subsequent inference on the regression parameter β would require a specific formulation of the $ϵ$-contaminated distribution. This is a nontrivial practical problem. Instead, there has been a focus on showing that the bias of estimators is bounded under contamination (Huber & Ronchetti, 2009, Section 4) while inference is conducted using the asymptotic distribution that assumes all observations are i.i.d. normal, implying that the ‘good’ observations are truncated normal. This gives a different distribution theory for inference compared with the one presented here and simulations in Section 7.1 indicate that this will not control size in general. It is therefore of interest to formulate models allowing ‘outliers’ under which consistency can be proved. The LTS model in this paper does so. The presented simulations show that the two inferential theories are really very different. In practice, LTS estimation should therefore be evaluated in the context of a particular model and inference should be conducted accordingly.

Alternative estimators of h. We have proposed consistent estimators for $h/n$⁠, but it would be useful to investigate their performance further. In a regression context, it may be worth considering the Forward Search algorithm (Atkinson et al., 2010). Omitted regressors and ‘outliers’ may confound each other, so a simultaneous search over these may be useful as in the Autometrics algorithm (Castle et al., 2021; Hendry & Doornik, 2014). Some asymptotic theory for these algorithms are provided in Johansen and Nielsen (2016a, 2016b).

Misspecification tests can be developed for the present model. The asymptotic theory developed here shows that standard normality tests can be applied to the set of estimated ‘good’ observations. Other tests could also be investigated, in particular those that are concerned with functional form or omission of regressors.

More ‘outliers’ than ‘good’ observations. This is allowed in the LTS model and supported by the asymptotic analysis under regularity conditions for the ‘outliers’. This lends some support to the practice of starting the Forward Search algorithm by an LTS estimator for fewer than half of the observations (Atkinson et al., 2010). Whether it makes more sense to model the ‘outliers’ or the ‘good’ observations as normally distributed in this situation must rest on a careful consideration of the data and the substantive context.

Acknowledgments

Comments from R. Butler, Y. Gao, O. Hao, C. Henning, referees, and the associate editors are gratefully acknowledged.

Data availability

All the data that we use are publicly available; the precise sources are cited in the article. The source code files in R available as Supplementary material.

Supplementary material

Supplementary material are available at Journal of the Royal Statistical Society: Series B online.

Appendix A. Derivation of the LTS likelihood

We prove that $P(zj−ϵ<yj≤zj|yifori∈ζ)=ΔϵGj(z~jβσ)=Gj(z~jβσϵ)−Gj(z~jβσ).$

First, we prove that $P0=P(yj≤zj|⋅)=G(z~jβσ)$⁠, where the dot indicates the conditioning set. Since $εj=(yj−β′xj)/σ$ and $zjβσ=(zj−β′xj)/σ$⁠, we have that $P0=P(εj≤zjβσ|⋅)$⁠.

Recall that the ‘outlier’ errors are $εj=(maxi∈ζyiβσ+ξj)1(ξj>0)+(mini∈ζyiβσ+ξj)1(ξj<0),$ see (3.1) and since $yiβσ=εi$⁠. Further, $ξj$ has distribution function $Gj$⁠, which is continuous at zero. As a consequence, $P(mini∈ζyiβσ≤εj≤maxi∈ζyiβσ|⋅)=0$ and $P(εj≤mini∈ζyiβσ)=P(ξj<0)=Gj(0)$⁠. Thus, $P0=P(εj≤zjβσ|⋅)=Gj(0)$ for $mini∈ζyiβσ≤zjβσ≤maxi∈ζyiβσ$⁠.

Let $P0=P1+P2$⁠, where $P1=P{(εj≤zjβσ)∩(ξj<0)|⋅}$ and $P2=P{(εj≤zjβσ)∩(ξj>0)|⋅}$⁠.

By (3.1), we have $(εj≤zjβσ)=(ξj≤zjβσ−mini∈ζyiβσ)$⁠, so that $P1$ can be written as $P{ξj≤min(zjβσ−mini∈ζyiβσ,0)|⋅}$⁠. Hence,

Similarly, $P2=P{(ξj≤zjβσ−maxi∈ζyiβσ)∩(ξj>0)|⋅}$ by (3.1). If $zjβσ<maxi∈ζyiβσ$⁠, then the intersection is empty. If instead $zjβσ>maxi∈ζyiβσ$⁠, then, the intersection is the set $(0<ξi≤zjβσ−maxi∈ζyiβσ)$⁠. Hence,

Note also that if $zjβσ<mini∈ζyiβσ$⁠, then $zjβσ<maxi∈ζyiβσ$⁠. And, if $zjβσ>maxi∈ζyiβσ$⁠, then $zjβσ>mini∈ζyiβσ$⁠. In combination, we have

Recall the notation $z~jβσ=(zjβσ−mini∈ζεi)1(zjβσ<mini∈ζεi)+(zjβσ−maxi∈ζεi)1(zjβσ>maxi∈ζεi)$⁠. We then get $P0=P(yj<zj|⋅)=Gj(z~jβσ)$⁠.

Second, similarly, $P(yj<zj−ϵ|⋅)=Gj(z~jβσϵ)$⁠. Combining, the desired result follows.

Appendix B. Proofs of asymptotic theory for the LTS model

We start with some preliminary extreme value results, which then allows analysis of the case with more ‘good’ observations than ‘outliers’. Then some results on empirical processes follows, which are needed for the general case.

B.1 Extreme values

For a distribution function $F$ define the quantile function $Q(ψ)=inf{c:F(c)≥ψ}$⁠.

    

B.2 Some initial properties of the LTS estimator

We consider the LTS estimator under the sequence of data generating processes defined in Section 5.1. The main challenge is to show that $δ^$ is close to the $Binomial(n−hn,ρ)$-distributed variable $δn=∑j∈ζn1(εj<mini∈ζnεi)$⁠. In the following lemmas, we condition on the sequence $δn$⁠, so that the randomness stems from ‘good’ errors $ε(δn+1),…,ε(δn+hn)$ and the magnitudes of the ‘outliers’, $ε_(δn+1−j)$ for $j≤δn$ and $ε¯(j−δn−hn)$ for $j>δn+hn$⁠. The unconditional statements in the Theorems about $δ^$ are then derived as follows. If $P(δ^−δn∈In|δn)→1$ for an interval $In$ and a sequence $δn$ then by the law of iterated expectations

due to the dominated converges theorem, because $P(δ^−δn∈In|δn)$ is bounded.

We give detailed proofs for the case $δ^>δn$⁠, so that some of the small ‘good’ observations are considered left ‘outliers’ and some of the small right ‘outliers’ are considered ‘good’. The case $δ^<δn$ is analogous, since we can multiply all observations by $−1$ and relabel left and right. When considering $δ^>δn$ we note that $δ^−δn≤n¯$⁠, the number of right ‘outliers’. Recall $ρ=G(0)$⁠. Hence, if $ρ=1$⁠, then all outliers are to the left. In particular, due to the binomial construction of $δn$ then $n¯=0$$a.s$⁠. when $ρ=1$⁠, so that the event $δ^>δn$ is a null set. Thus, when analysing $(δ^−δn)/hn$ it suffices to consider $δ^>δn$ and $ρ<1$⁠.

The next lemma concerns cases where $s^=δ^−δn$ is small. We show that the LTS estimators are close to the infeasible OLS estimators for the ‘good’ observations. We note that, by Lemma B.2 and Assumption 5.1(ii) that $E|εi|2+ω<∞$⁠, we can find a small $η>0$ so that $ε(δn+1),ε(δn+hn)=oP(hn1/2−η)$⁠. We write $ε(δ^+i)p$ for ${ε(δ^+i)}p$⁠.

  

The next Lemma is the main ingredient to showing consistency of $δ^$⁠. It is convenient to define sequences

We note that by Assumption 5.1(iii), $εi$ has a distribution function with infinite support. Hence, the extreme value $ε(δn+hn)$ diverges to infinity and $ε(δn+hn)−1=oP(1)$⁠. By Assumption 5.1(ii), Lemma B.2 applies, so that $ε(δn+hn)=oP(hn1/2−η)$ for some $η>0$⁠. Then, for large $hn$ and η small, $0<s_n<s¯n$ where $s¯n/hn=oP(1)$⁠.

  

B.3 Fewer ‘outliers’ than ‘good’ observations

 

B.4 Marked empirical processes evaluated at quantiles

We start with some preliminary results on marked empirical processes evaluated at quantiles.

Consider random variables $εi$ for $i=1,…,n$ and define the marked empirical distribution and its expectation, for $c≥0$⁠, by

For $p=0$⁠, let $Fn0=Fn$⁠. We also define the quantile function $Q(ψ)=inf{c:F(c)≥ψ}$ and the empirical quantiles $Qn(ψ)=inf{c:Fn(c)≥ψ}$⁠. The first result follows from the theory of empirical quantile processes.

  

We need the following Glivenko–Cantelli result.

  

The next result is inspired by Johansen and Nielsen (2016a, Lemma D.11).

  

B.5 More ‘outliers’ than ‘good’ observations

The next lemma is needed when there are more than half of the observations are ‘outliers’. As $σ^δ2$ is not diverging, additional regularity conditions are needed to ensure that $σ^δ2>σ^δn2$⁠.

   

B.6 The OLS estimator in the LTS model

 

References

Author notes