Robust Estimation of Regression Models with Potentially Endogenous Outliers via a Modern Optimization Lens

Robust estimation of linear regression models in the presence of outliers has been extensively studied in econometrics and statistics. It is well-known that the classical ordinary least square (OLS) estimator is extremely sensitive to outliers. Various robust estimation methods have been proposed, for example, Huber’s M-estimation (Huber, 1964), the least median of squares (LMS) estimator (Siegel, 1982), and the least trimmed squares (LTS) estimator (Rousseeuw, 1984, 1985), among others. More recently, Lee et al. (2012) propose regularization of the $L_{1}$-norm of case-specific parameters to achieve robust estimation against outliers and She and Owen (2011) show that any regularized estimator can be formulated as an equivalent iterative thresholding procedure of outlier detection. Refer to Yu and Yao (2017) for a comprehensive survey and detailed comparison.

The traditional literature has focused on the breakdown point and efficiency of a robust estimator (Huber and Donoho, 1983). The finite sample breakdown point of an estimator measures the proportion of outliers that can be arbitrarily contaminated before the estimation error goes to infinity. On the other hand, efficiency measures the relative estimation efficiency of the robust estimator compared to OLS in the ideal scenario where the error term is normally distributed and there are no outliers. However, there is limited work investigating how robust estimators perform when the outliers are endogenous in the sense that the noises can be arbitrarily correlated with observed regressors. In this scenario, an outlier not only brings in a contamination noise of large magnitude but also introduces model misspecification against the classical linear regression assumptions. The endogeneity issues from the outliers can lead to severe estimation bias.

In this paper, we first examine the widely adopted $L_{1}$-regularized estimator, which includes the Huber estimator and least absolute deviation as special cases, and demonstrate that it is subject to significant bias when outliers are endogenous. As shown in She and Owen (2011), the $L_{1}$-regularized estimation is equivalent to an iterative soft-thresholding procedure with a data-driven thresholding parameter, and hence it does not completely eliminate the detected outlier.

Motivated by this finding, we turn to the $L_{0}$-regularized approach. By restricting the cardinality of the set of outliers, the $L_{0}$-regularized estimator searches for the best subset of observations and the regression coefficients that minimize the least squares objective function. This problem is equivalent to solving the least trimmed squares (LTS) estimator, which is well-known to be an NP-hard combinatorial optimization problem (Natarajan, 1995). Following Bertsimas et al. (2016); Thompson (2022), the $L_{0}$-regularized estimator can be formulated as a mixed integer optimization (MIO) problem. With the continuous advancement in modern optimization software, such as gurobi¹¹1gurobi is a commercial optimization solver that specializes in mixed integer optimization. In the past decade, gurobi has been up to 75-fold faster in computation speed for mixed integer optimization problems independent of hardware advancement (Gurobi Optimization, LLC, 2024)., and computational infrastructure, it is now tractable to solve the $L_{0}$-regularized estimation problem with a verifiable global optimal solution for datasets with up to hundreds of observations. However, the optimization routine heavily relies on the initial values because of its non-convexity nature. In addition, the implementation is not scalable, and the computational burden is increasingly heavy as the sample size and fraction of outliers grow.

To address the computational challenges, we propose a scalable and efficient heuristic algorithm that combines the iterative hard-thresholding (IHT) algorithm and a local combinatorial search refinement inspired by Hazimeh and Mazumder (2020). Given the initial solution by the IHT algorithm, the local combinatorial search algorithm checks whether swapping a few observations, say one or two, between the estimated set of inliers and outliers, can improve the objective value by solving a small-scale mixed integer optimization problem. The heuristic algorithm is more computationally efficient and can improve the initial IHT solution to be as good as the solution from solving the original $L_{0}$-regularized MIO problem in most cases, according to the Monte Carlo experiments.

Our contributions are twofold. First, we document that existing $L_{1}$-regularized methods are biased in the presence of endogenous outliers, whereas the $L_{0}$-regularized approach does not suffer from this bias through Monte Carlo simulations. This finding extends the understanding of the properties of robust estimation methods to a new scenario. Second, we propose systematic heuristic algorithms that provide stable and high-quality solutions while being computationally efficient. To the best of our knowledge, this is the first work to apply the idea of local combinatorial search to the context of robust estimation and outlier detection.

To illustrate the practical value of our method, we apply it to an empirical application in stock return forecasting. The results demonstrate that our $L_{0}$-regularized estimator outperforms $L_{1}$-regularized alternatives in terms of out-of-sample prediction errors.

The rest of the paper is organized as follows. Section 2 presents the model setup and motivating examples. Section 3 summarizes the existing $L_{1}$-regularized methods for robust estimation. Section 4 delineates the $L_{0}$-regularized robust methods and the detailed algorithms. The performance of the algorithms is evaluated through Monte Carlo experiments in 5. Section 6 illustrates the proposed method in an empirical application of stock return forecasting. Section 7 concludes.

This section sets up the linear regression models with potentially endogenous outliers. Suppose we observe samples $(y_{i},x_{i}^{\prime})^{\prime}$ for individual unit $i\in\left[N\right]$. Let $\mathcal{O}\subset\left[N\right]$ and $\mathcal{I}=\left[N\right]\setminus\mathcal{O}$. Consider the following linear regression model:

where $\mathbb{E}(u_{i}|x_{i})=0$ for all $i\in\left[N\right]$. The parameter $\alpha_{i}$ represents the conditional mean shift such that

We refer to the parameters $\alpha_{i}$ as outlier fixed effects, which are allowed to be arbitrarily correlated with the regressors $x_{i}$. An alternative formulation of model (2.1) is:

where $\gamma_{i}=\mathbbm{1}\left\{i\in\mathcal{O}\right\}$ is the outlier dummy and $\widetilde{\alpha}$ is the latent outlier fixed effect.

The existence of $\alpha_{i}$ introduces endogeneity to a subset of observations $\mathcal{O}$. The classical ordinary least squares (OLS) estimator is expected to be biased in the presence of endogenous outliers. The following motivating examples demonstrate possible sources of endogenous outliers.

The objective is to accurately estimate $\beta=(\beta_{0},\beta_{1}^{\prime})^{\prime}$ in the presence of outliers in the samples.

The most widely used robust estimators of $\beta$ in the presence of outlier observations are Huber’s M-estimation and the least absolute deviation (LAD) estimation method (Huber, 1964). These two estimators can be understood as special cases of the $L_{1}$-regularized estimator.

Denote $Y=(y_{1},...,y_{N})^{\prime},X=\left((1,x_{1}^{\prime}),...,(1,x_{N}^{\prime})\right)^{\prime},\beta=\left(\beta_{0},\beta_{1}^{\prime}\right)^{\prime},\alpha=(\alpha_{1},...\alpha_{N})^{\prime},U=(u_{1},...,u_{N})^{\prime}$. Then, we can present the model (2.1) in matrix form as

Let the least squares loss function be

and consider the following $L_{1}$-regularized optimization problem,

Given $\beta$, $\hat{\alpha}^{\psi}(\beta)\coloneq\operatorname*{arg\,min}_{\alpha}Q_{N}^{\psi}(\beta,\alpha).$ The closed-form solution is

Then the profile $L_{1}$-regularized objective function becomes

In (3.3), the tuning parameter $\psi$ plays the role of a soft thresholding parameter. When $\psi$ is fixed, (3.4) becomes the objective function of Huber regression

where

is the Huber loss function (Hastie et al., 2015, p. 26 - 27). As illustrated in Figure 1, the cutoffs $-\psi$ and $\psi$ divide the Huber loss function into two regimes: it is a squared loss if $|t|\leq\psi$ and a linear loss in absolute deviation if $|t|>\psi$. If we let $\psi=\psi_{N}$ so that $\min_{i\in\left\{1,2,\cdots,N\right\}}\left|y_{i}-x_{i}^{\prime}\beta\right|\geq\psi_{N}>0$ given $N$ and the sample $\left(Y,X\right)$ , then we have

as $\psi_{N}\to 0$, which shrinks the squared loss region in Figure 1.

The $L_{1}$-regularized estimator unifies the Huber and LAD estimators and offers greater flexibility by allowing the regularization parameter to be chosen in a data-driven manner. However, as a soft-thresholding algorithm, as shown in (3.3), it does not completely eliminate the detected outliers from the estimation procedure. This can lead to estimation biases when the outliers are endogenous.

As demonstrated in the Monte Carlo experiments in Section 5, the $L_{1}$-regularized estimator performs well in terms of estimation bias, root mean squared error (RMSE), and prediction error in DGP 1, where the outlier fixed effects are generated as exogenous random variables. However, in DGP 2 and DGP 3, where the outlier fixed effects are correlated with the regressors, the $L_{1}$-regularized estimator for the slope coefficients exhibits significant biases. This finding motivates the consideration of the hard-thresholding-based $L_{0}$-regularized estimator, as discussed in Section 4.

We consider the $L_{0}$-reguarization on the outlier fixed effects in the least square estimation,

where $\left\lVert\alpha\right\rVert_{0}=\sum_{i=1}^{N}\mathbbm{1}\left\{\alpha_{i}\neq 0\right\}$ and the tuning parameter $k\in\left[N\right]$ controls the exact sparsity of the outlier fixed effects $\alpha$.

As in Bertsimas et al. (2016); Thompson (2022), (4.1) can be formulated as the following mixed integer optimization (MIO) problem,

in which we introduce a binary variable $\gamma_{i}$ to model whether an observation is detected as an outlier. $\sum_{i=1}^{N}\gamma_{i}\leq k$ corresponds to the cardinality constraint $\left\lVert\alpha\right\rVert_{0}\leq k$. If $\gamma_{i}=1$, i.e. $i$ is labelled as an outlier, then $\alpha_{i}$ should be nonzero; on the other hand if $\gamma_{i}=0$, i.e. $i$ is labelled as an inlier, then $\alpha_{i}$ should be exactly 0. This intuition can be translated into the constraint $\left(1-\gamma_{i}\right)\alpha_{i}=0$, which can be further modeled via integer optimization using the Special Ordered Set of Type 1 (SOS-1), that is a set contains at most one non-zero variable.²²2 (4.1) can also be formulated as a MIO problem using big-$M$ method, $$\min_{\beta,\alpha,\gamma}L_{N}\left(\beta,\alpha\right)\text{ s.t. }\gamma_{i}\in\left\{0,1\right\},\,\sum_{i=1}^{N}\gamma_{i}\leq k,\,-\gamma_{i}M_{\alpha}\leq\alpha_{i}\leq\gamma_{i}M_{\alpha},$$ where $M_{\alpha}$ is provides a sufficiently large bound for $\alpha$. In practice, we set $M_{\alpha}=\tau\left\lVert\widehat{\alpha}\right\rVert_{\infty}$ with $\tau>1$ and $\widehat{\alpha}$ is the initial estimator form the heuristic algorithm. $M_{\alpha}$ and $M_{\alpha,1}$ are user-defined bound parameters that can help to tighten the parameter space and improve the computation performance.

The optimization problem (4.4) is a well-posed MIO problem, and provable global optimality can be achieved by optimization solvers such as gurobi (Gurobi Optimization, LLC, 2024). The computational efficiency and solution quality are highly dependent on the initial values provided to the solver. Furthermore, as demonstrated in the Monte Carlo experiments in Table 1, directly solving (4.4) does not achieve provable optimality within the prespecified 5-minute timeframe when the sample size $N$ and the number of outliers $k$ are large. In the following subsection, we propose systematic approximate algorithms to solve (4.1).

In this section, we examine the numerical performance of the proposed $L_{0}$-regularized estimation procedure. We compare its coefficient estimation accuracy and prediction error to those of the $L_{1}$-regularized estimation and the classical methods, LAD and OLS.