Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
80,542 characters · 15 sections · 42 citation commands
Robust Estimation of Regression Models with Potentially Endogenous Outliers via a Modern Optimization Lens
\thispagestyle{empty} \setcounter{page}{1}
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 huber1964robust, the least median of squares (LMS) estimator siegel1982robust, and the least trimmed squares (LTS) estimator rousseeuw:1984:lms_lts,rousseeuw1985multivariate, among others. More recently, lee2012regularization propose regularization of the \(L_1\)-norm of case-specific parameters to achieve robust estimation against outliers and she2011outlier show that any regularized estimator can be formulated as an equivalent iterative thresholding procedure of outlier detection. Refer to yu2017robust for a comprehensive survey and detailed comparison.
The traditional literature has focused on the breakdown point and efficiency of a robust estimator huber1983notion. 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 she2011outlier, 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:nphard. Following bertsimas/king/mazumder:2016:bestsubsetmio,thompson2022robust, 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\footnote{gurobi 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.}, 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/mazumder:2020:fastsubset. 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.
{\it Notations}. For a generic vector $a = \left( a_1, a_2, \cdots, a_N \right)^{\prime}$, $\left\Vert a \right\Vert = \left( a^{\prime} a \right)^{1/2}$, $ \left\Vert a \right\Vert_1 = \sum_{i=1}^N \left\vert a_i \right\vert $, \(\left\lVert a \right\rVert_\infty = \max_i \vert a_i \vert \) and \(\left\lVert a \right\rVert _{0} = \sum_{i=1}^{N} \mathbbm{1} \left\{ a_i \neq 0 \right\} \). For matrix $A$, we define the Frobenius matrix norm $\left\Vert A \right\Vert$ as $\left\Vert A \right\Vert = \left( \operatorname{tr}\left( A^{\prime}A \right) \right)^{1/2}$. Generically, $[N] = \left\{ 1, 2, \cdots, N \right\}$ for positive integer $N$. For \(a \in \mathbb{R} \), \( \lfloor a \rfloor \) is the integer part of the real number \(a\). \(\iota_N = \left( 1,1,\cdots, 1 \right)\in \mathbb{R} ^N \) denotes the one vector.
The rest of the paper is organized as follows. Section (ref) presents the model setup and motivating examples. Section (ref) summarizes the existing \(L_1\)-regularized methods for robust estimation. Section (ref) delineates the \(L_0\)-regularized robust methods and the detailed algorithms. The performance of the algorithms is evaluated through Monte Carlo experiments in (ref). Section (ref) illustrates the proposed method in an empirical application of stock return forecasting. Section (ref) concludes.
This section sets up the linear regression models with potentially endogenous outliers. Suppose we observe samples $(y_i,x_i')'$ 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 (ref) 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')'$ 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 huber1964robust. These two estimators can be understood as special cases of the $L_1$-regularized estimator.
Denote $Y = (y_1,...,y_N)', X = \left((1,x_1'),...,(1,x_N')\right)', \beta = \left( \beta_0, \beta_1^\prime \right)^\prime, \alpha = (\alpha_1,...\alpha_N)', U = (u_1,...,u_N)'$. Then, we can present the model ((ref)) 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 (ref), the tuning parameter \(\psi \) plays the role of a soft thresholding parameter. When \(\psi\) is fixed, (ref) becomes the objective function of Huber regression
where
is the Huber loss function hastie/tibshirani/wainwright2015. As illustrated in Figure (ref), 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\vert y_i - x_i^{\prime} \beta \right\vert \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 (ref).
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 (ref), 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 (ref), 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 (ref).
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/king/mazumder:2016:bestsubsetmio,thompson2022robust, (ref) 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.\footnote{ (ref) can also be formulated as a MIO problem using big-\(M\) method,
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 (ref) is a well-posed MIO problem, and provable global optimality can be achieved by optimization solvers such as gurobi gurobi. 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 (ref), directly solving (ref) 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 (ref).
The same as best subset selection studied in bertsimas/king/mazumder:2016:bestsubsetmio, $ L_0 $-regularization robust regression (ref) cannot be solved in polynomial time, i.e. it is NP-hard natarajan:1995:nphard. The computational cost of solving for the exact solution increases steeply as the sample size increases. Heuristic algorithms that can find approximate solutions are useful for parameter tuning and providing warm-starts for the solvers. In this section, we follow bertsimas/king/mazumder:2016:bestsubsetmio,hazimeh/mazumder:2020:fastsubset,thompson2022robust,mazumder/radchenko/dedieu:2022:subsetselection to provide heuristic algorithms based on iterative hard-thresholding, neighborhood search, local combinatorial search for (ref).
Iterative hard-thresholding (IHT) based on project gradient descent is a generic heuristic algorithm for general nonlinear programming with sparsity constraints beck/eldar:2013:sparsityconstraint and it is successfully applied to best subset selection in bertsimas/king/mazumder:2016:bestsubsetmio. For problem (ref), we propose the following iterative hard-thresholding algorithm, which is a simplified version of the projected block-coordinate gradient descent in thompson2022robust.
Define the hard-thresholding operator for $ c\in \mathbb{R}^{N} $ as
As shown in bertsimas/king/mazumder:2016:bestsubsetmio, if $ \hat{\alpha} \in H_k\left( c \right)$, then $ \hat{\alpha} $ retains the $ k $ largest (in absolute value) elements of $ c $ and sets the rest to zero, i.e.
if $ \left\{ (1), (2), \cdots, (N) \right\}$ is an ordering of $ \mathcal{N} $ such that $ \left\vert c_{(1)} \right\vert \geq \left\vert c_{(2)} \right\vert \geq \cdots \geq \left\vert c_{(N)} \right\vert$.
The iterative hard-thresholding algorithm starts with an initial robust estimator of $ \beta $ and iteratively updates $ \alpha $ and $ \beta $ by applying the hard-thresholding operator to the residuals and running least square estimation based on the support of updated $ \alpha $. The details are summarized in Algorithm (ref).
Despite the iterative hard-thresholding algorithm is fast and easy to implement, it is only guaranteed to converge to a coordinate-wise optimal solution. In this section, we propose a local combinatorial search algorithm, similar to hazimeh/mazumder:2020:fastsubset, to further refine the solution.
Note that (ref) is equivalent to
For an estimate $ \hat{\alpha} $, denote $ \hat{\mathcal{I}} = \left\{ i\in \mathcal{N}: \hat{\alpha}_i = 0 \right\}$ and $ \hat{\mathcal{O}} = \left\{ i\in \mathcal{N}: \hat{\alpha}_i \neq 0 \right\}$. Following hazimeh/mazumder:2020:fastsubset, $ \hat{\alpha} $ is said to be swap-inescapable of order $ l $ if arbitrarily swapping up to $ l $ observations between $ \hat{\mathcal{I}} $ and $ \hat{\mathcal{O}} $ and then optimizing over the new support cannot improve the objective value. Formally,
If solution $ \hat{\alpha} $ is swap-inescapable of order $ k $, then $ \hat{\alpha} $, associated with the OLS estimates $ \hat{\beta} $ using observations in $ \hat{\mathcal{I}} $, is the exact global optimal solution to (ref). By solving the local combinatorial search problem in (ref) for small $ l < k $, say $ l = 1$ or $2$, we improve the IHT algorithm and verify the local combinatorial exactness.
The problem in (ref) can be formulated as a mixed-integer optimization (MIO) problem given by
The constraints \(\sum_{i\in\hat{\mathcal{I}}} \gamma_i \leq l\) and \(\sum_{i\in\hat{\mathcal{O}}} \gamma_i \geq k - l\) restrict the number of swaps between \(\hat{\mathcal{I}} \) and \(\hat{\mathcal{O}} \) up to \(l\). (ref) has a much smaller search space than the original problem (ref) when $ l $ is small.\footnote{ When $ l = 1 $, we can solve the problem in a brute force way by running least square estimation on all possible 1-1 swaps between $ \hat{\mathcal{I}} $ and $ \hat{\mathcal{O}} $ and comparing the $ k\left( N - k \right) $ resulting objective values without invoking the MIO solver. } The heuristics combining IHT and local combinatorial search are summarized in Algorithm (ref).
Algorithm (ref) is guaranteed to provide a solution that is locally optimal in the sense of swap-inescapability of order $ l $. However, the quality of the solution depends on the initial inputs due to the noncovexity. As noted in thompson2022robust and mazumder/radchenko/dedieu:2022:subsetselection, a neighborhood search procedure can serve as a useful systematic way of perturbing the initial inputs to improve the solution quality of the heuristic algorithm. Let $ [K] = \{1, 2, \cdots, K\}$ with $ K \leq \lfloor N/2\rfloor $ be the set of candidate sparsity parameters.
Denote $ b_{k, l}\left( \beta \right) $ and $ a_{k, l}\left( \beta \right) $ as the outputs of Algorithm (ref) initialized with $ \beta $, $ k $ and $ l $. The neighborhood search procedure is summarized in Algorithm (ref). As a by-product of the procedure, we obtain a solution for each $ k\in \left[ K \right] $ which is useful for sensitivity analysis and parameter tuning, which is discussed in Section (ref).
In both \(L_1\)- and \(L_0\)- regularized estimation methods, the selection of tuning parameters, \(\psi \) in (ref) and \(k\) in (ref), plays a critical role. We propose the BIC-type information criteria,
For computational efficiency, we use Algorithm (ref) to generate solutions for a grid of candidate tuning parameters and select the one that minimizes the BIC,
where $K \leq \lfloor N/2\rfloor $ is the maximum potential number of outliers.
The widely used theoretical property to evaluate robust estimation methods is the finite sample breakdown point, proposed in hampel1971general and huber1983notion. Suppose the original sample is $\left( Y, X \right)$ and the contaminated sample is $\left( \tilde{Y}_{(k_0)}, \tilde{X}_{(k_0)} \right)$ with $k_0$ observations in the original sample being arbitrarily replaced by outliers. The finite sample breakdown point of an estimator $T$ is defined as
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.
Following the setting in Section (ref), we consider the following data generating processes (DGPs).
DGP 1 (Exogenous Outliers). Consider the linear regression model with outliers,
where $\gamma_i$ is the indicator for outliers and $ \alpha_i $ is the outlier fixed effect. We generate the regressors by $x_{i, 1} = \left(v_{i,1}^2 + v_{i,2}^2 - 2\right)/2$ and $ x_{i, 2} = x_{i, 1} + v_{i, 3} $, where $v_{i, j} \sim \text{i.i.d.}N(0, 1)$ for $ j = 1, 2, 3 $, and the error term by $ u_i \sim \text{i.i.d.}N(0, 1)$. Let $p\in \left( 0, 1 \right)$ denote the fraction of outliers. $k_0 = \left\lfloor pN \right\rfloor$ and $ \gamma_i =
. $ Let $\alpha_i \sim i.i.d.N(\mu_\alpha, \sigma_\alpha^2)$ be exogenous shocks to outliers where $\left( \mu_\alpha,\sigma_\alpha \right) \in \left\{ (0, 5), (5, 5), (10, 10) \right\}$. The true coefficients are $ \beta_0 = 0.5 $, $ \beta_1 = \left( 1, 1 \right)^\prime $.
\noindentDGP 2 (Endogenous Outliers). This DGP deviates from DGP 1 by allowing the outlier fixed effects to be correlated with the regressors. Let $ \alpha_i = \rho \left( v_{i, 1} + v_{i, 2}+ v_{i, 3} \right) $ be a linear combination of the innovations to make the outlier fixed effects correlated with regressors and create endogeneity. The parameter $ \rho \in \left\{ 2, 5, 10 \right\} $ to control the degree of correlation. The rest components are the same as DGP 1.
{DGP 3} (Predictive Regression with Endogenous Shocks). Consider a linear predictive regression model as studied in kostakis2015robust,koo2020high,lee2022lasso. The dependent variable is generated as
where $\beta _{0} = 0.3$, $\beta _{1} = 1$, \(\phi = \left( 1, -1 \right) \) and \(\eta _{N} = \left( 1 / \sqrt{N}, -1 / \sqrt{N} \right) \). The vector of the stacked innovation $\xi_{i}=\left(z_{i},\underset{2\times1}{v_{i}^{\prime}},\underset{2\times1}{e_{i}^{\prime}},u_{i}\right)^{\prime}$ follows a VAR(1) process $\xi_{i}=\Phi\xi_{i-1}+\varepsilon_{i}$, where $\varepsilon_{i}\sim iid\;N\left(0,\Sigma_{\varepsilon}\right)$ in which $\Phi$ and $\Sigma_{\varepsilon}$ are empirically estimated from the welch2008comprehensive data as in lee2022lasso. $x_{i}^{c}\in\mathbb{R}^{2}$ is a vector I(1) process with cointegration rank 1 based on the VECM, $\Delta x_{i}^{c}=\Gamma^{\prime}\Lambda x_{i-1}^{c}+v_{i},$ where $\Lambda=
$ and $\Gamma=
$ are the cointegrating matrix and the loading matrix, respectively. $\left(x_{i,l}\right)_{l=1}^{2}$ are random walks generated by $x_{i,l}=x_{i-1,l}+e_{i,l}$, $l=1,2$. Let \(p \in (0,1)\) be the fraction of outliers and define \(k_0 = \lfloor pN \rfloor\). We impose two outlier periods. Let \(c_1 = \lfloor 0.25 N \rfloor\) and \(c_2 = \lfloor 0.75 N \rfloor\) be the positions in the sample where the outliers start. The outlier indicator \(\gamma_i = 1\) for \(i = c_l + 1, c_l + 2, \cdots, c_l + \lfloor k_0 / 2 \rfloor\), \(l = 1, 2\) and 0 otherwise. The outlier shifts \(\alpha _{i} = \rho \left( z_i + v_{i} ^{\prime} \iota_2 \right) \) where \(\rho \in \left\{ 2,5,10 \right\} \).
To evaluate each heuristic algorithm in Section (ref) as compared to solving the full-scale mixed integer optimization for (ref), we will report the measure Equal MIO, which is the frequency among replications the estimates obtained by the algorithm are the same as those of the MIO solution of (ref). For each implementation invoking the MIO solver, we report relative optimality gap, which is defined as \[ \text{Relative Optimality Gap} = \frac{f^{P} - f^{D}}{f^{D}}, \] where \(f^{P}\) is the attained upper bound of the objective value of the MIO solution and \(f^{D}\) is the dual lower bound of the objective value delivered by the solver.\footnote{The value can be read from the solver output. Details of the definition can be found at \href{https://www.gurobi.com/documentation/current/refman/mipgap2.html}{https://www.gurobi.com/documentation/current/refman/mipgap2.html}.} Relative optimality gap measures the progress of optimality verification of the solver. The solution is verified to be globally optimal if this gap is less than a prespecified tolerance threshold \(10^{-4}\).
For the estimation and out-of-sample prediction performance, we report the bias, root mean squared error (RMSE) and prediction errors. Generically, bias and RMSE are calculated by $R^{-1}\sum_{r=1}^R \left( \hat{ \beta}^{(r)} - \beta \right)$ and $\sqrt{ R^{-1}\sum_{r= 1}^R \left( \hat{\beta}^{(r)} -\beta \right)^2}$, respectively, for true parameter $\beta$, its estimate $\hat{\beta} ^{(r)}$ across $R$ replications. To calculate the prediction error, we generate test data following the same DGP without outliers, \(\left\{ \widetilde{y} _{i}, \widetilde{x}^{\prime}_{i} \right\}_{i=1}^{N_t} \), \(N_{t} = 1000\). The prediction error is calculated by \[ \text{Prediction Error} = R^{-1} \sum_{r=1}^R \left( \frac{1}{N_t} \sum_{i=1}^{N_t} \left( \widetilde{y} _{i}^{(r)} - \widetilde{x} _{i}^{(r)\prime} \hat{\beta}^{(r)} \right)^2 \right). \]
In all experiments, we consider $ p\in \left\{ 5\%, 10\%, 20\% \right\} $ and $ N \in \left\{ 100, 200, 400 \right\} $. All experiments are carried out on a Linux machine with an Intel i9-13900K CPU and Guorbi optimization solver at version 11.0. The time limit for the solver is set to 5 minutes.
We begin by evaluating the performance of the heuristic algorithms presented in Section (ref) in comparison to solving the full-scale mixed integer optimization problem (ref). We set the parameters \(\left( \mu _\alpha , \sigma _\alpha \right) = \left( 5,5 \right) \) for DGP 1 and \(\rho = 5 \) for DGP 2 and 3 and run \(100\) replications for each setup. The parameter \(k\) is set to the true value \(k_0\) for all algorithms. For MIO problems (ref) and (ref), \(M_\alpha = \tau \left\lVert \hat{\alpha} \right\rVert_\infty \) and \(M_{\alpha ,1} = \tau \left\lVert \hat{\alpha} \right\rVert_1 \) where \(\hat{\alpha} \) is the initial estimator and \(\tau = 1.5\).
For each simulated dataset, we implement Algorithm (ref), the iterative hard-thresholding (IHT) method. Subsequently, using the IHT estimate as the initial estimator, we apply the local combinatorial search (LCS) algorithm with local exactness levels \( l = 1\) and \( l = 2\), referred to as LCS-1 and LCS-2, respectively. Additionally, we solve the full-scale mixed integer optimization problem (ref) using the IHT estimate as the warm-start. In Table (ref), we report the average CPU runtime (in seconds) for each algorithm, the frequency of IHT, LCS-1 and LCS-2 estimates are equal to MIO and the average relative optimality gap.
The key findings presented in Table 1 highlight the performance differences between the heuristic algorithms and the mixed integer optimization (MIO) approach. Notably, within the 5-minute time limit, the MIO method fails to complete the optimality verification for cases where the scale is larger than \(N = 200\) and \(p = 0.1\). In contrast, the iterative hard-thresholding (IHT) and the local combinatorial search (LCS-1) with \(l=1\) are significantly faster. Specifically, IHT runs in milliseconds and LCS-1 achieves optimality in seconds for most cases. LCS-2 takes less than a minute except when \(N = 400\) and \(p = 0.1\) or \(0.2\). In these exceptional cases, the relative optimality gap of LCS-2 remains below 1%, indicating near convergence. Furthermore, the initial estimator, IHT, provides solutions of the same quality to MIO in 50% to 90% of the cases. When advancing to LCS-1, the frequency of obtaining solutions equivalent to MIO increases to approximately 90% while with a significantly lower computational cost. LCS-2 consistently matches the MIO solutions in nearly all replications. These findings underscore the effectiveness and computation efficiency of the heuristic algorithms. Particularly, the local combinatorial search algorithm balances the solution stability and computation costs, which is important when dealing with large sample sizes and many outliers.
In this section, we compare the performance of \(L_0\) and \(L_1\)-regularized estimation in different scenarios. For the \(L_0\) and \(L_1\)-regularized estimation, the tuning parameters \(k\) and \(\psi\) are selected using the information criteria (ref). When choosing \(k\), we rely on the neighborhood search algorithm, as detailed in Algorithm (ref) to generate the estimates for each candidate \(k \in \left[ K \right] \). Specifically, in the implementation of Algorithm (ref), we set \(K = 2k_0\) and the local exactness level \(l = 1\). We report LCS-2 with \(\hat{k}\) selected by BIC as the \(L_0\)-regularized estimation results. In addition, we include ordinary least square (OLS) and least absolute deviation (LAD) as benchmarks. Bias and root mean squared error (RMSE) are reported for the first slope parameter \(\beta_1\) in all data generating processes (DGPs). For each setup, we run 1000 replications.
The key findings from the Monte Carlo simulation reported in Table (ref) to (ref) reveal several important insights. In DGP 1, where the outlier fixed effects are exogenous and there is no endogeneity issue, the bias of the \(L_0\) method is slightly larger than that of the alternatives. However, the \(L_0\) method achieves the smallest RMSE and prediction error in most cases. Ordinary Least Squares (OLS) is not robust to outliers, and their RMSE and prediction error increase significantly as the magnitude of the outlier fixed effects becomes larger. The estimation accuracy of the \(L_1\) method is unstable when the sample size and the outlier fraction are small, but it outperforms Least Absolute Deviations (LAD) and OLS as \(N\) and \(p\) increase. In DGP 2 and 3, where endogenous outlier fixed effects are introduced, the \(L_1\)-regularized method, LAD, and OLS exhibit severe estimation bias. In contrast, the \(L_0\)-regularized method is free from estimation bias and achieves the smallest RMSE and prediction error. The accuracy gap widens as \(\rho\), the parameter controlling the degree of endogeneity, and the fraction of outliers increase. These findings demonstrate that the \(L_0\) method is more robust to the presence of endogenous outliers and provides more accurate estimation.
Linear predictive regression models have been extensively studied for stock return forecasting. For instance, koo2020high,lee2022lasso and others have applied linear predictive regression to the welch2008comprehensive dataset to investigate stock return predictability. Financial time series data often exhibit instability, particularly during periods such as the financial crisis. Including these periods in estimation and forecasting can lead to varying results. As an illustration, we apply data-driven \(L_0\) and \(L_1\)-regularized robust estimation methods to the welch2008comprehensive dataset\footnote{Retrieved from \url{http://www.hec.unil.ch/agoyal/}} and evaluate the out-of-sample stock return prediction performance.
We use monthly data from January 1990 to December 2023, which covers the dot-com bubble period, the 2007-09 financial crisis, and the Covid-19 pandemic. The dependent variable, excess return, is defined as the difference between the continuously compounded return on the S&P 500 index and the three-month Treasury bill rate, computed by \[ \text{ExReturn}_{i} = \log\left(\text{index}_{i}/\text{index}_{i-1}\right) - \log\left(1 + \text{tbl}_{i}/12\right). \] Twelve financial and macroeconomic variables\footnote{ The predictors are Dividend Price Ratio, Dividend Yield, Earning Price Ratio, Term Spread, Default Yield Spread, Default Return Spread, \texttt{Book-to-Market Ratio}, \texttt{Treasury Bill Rates}, \texttt{Long-Term Return}, \texttt{Net Equity Expansion}, \texttt{Stock Variance}, \texttt{Inflation}. Table (ref) in the appendix summarizes the detailed description of each variable. }, denoted by \(x_{i}\), are included in the model as predictors. We apply the \(L_0\) and \(L_1\)-regularized methods to estimate the model, \[ \text{ExReturn}_{i+1} = \beta_0 + x_{i}^{\prime}\beta_1 + \alpha_{i} + u_{i+1}, \] for \(i \in \left[ N \right]\), and construct one-month-ahead forecasts \(\hat{f}_{N+1} = \hat{\beta}_0 + \hat{x}_{N}^{\prime}\hat{\beta}_1\) recursively with a 10-year rolling window for each month from January 2000 to December 2023. As in Monte Carlo experiments, LAD and OLS are included as benchmarks.
Table (ref) compares the mean prediction squared error (MPSE) among different methods across various forecasting periods. We consider the forecast period from January 2000 to December 2023, and six subperiods based on the start and end dates of the dot-com bubble burst period (Mar. 2000 to Nov. 2000), the financial crisis (Dec. 2007 to Jun. 2009) and the Covid-19 shocks (Feb. 2020 to Apr. 2020). The results demonstrate that \(L_0\)-regularized method outperforms alternative methods with a margin in all subperiods except for the financial crisis period, which illuminates the forecast accuracy gain from the robustness to potentially endogenous outliers.
Additionally, Figure (ref) illustrates the outliers detected (\(\hat{\alpha}_{i} \neq 0\)) for each rolling window using \(L_0\) and \(L_1\)-regularized methods. In the grid plot, each row corresponds to a different forecast period, and each column corresponds to a period that may be included in the estimation rolling window. Within each row, the highlighted cells represent the estimation window for the corresponding forecast target period. A cell is labeled red if the corresponding period is detected as an outlier in the rolling window by both methods, blue if both methods detect this period as an inlier, purple if only the \(L_0\)-regularized method detects this period as an outlier, and yellow if only the \(L_1\)-regularized method detects this period as an outlier. As shown in the figure, the \(L_0\) method consistently detects periods around the dotcom bubble, financial crisis, and Covid-19 shocks across rolling windows. In general, the \(L_1\) method detects a similar pattern while also labeling some outlier periods outside the concentrated outlier regions.
This paper addresses the robust estimation of linear regression models in the presence of potentially endogenous outliers. We demonstrate that existing $L_1$-regularized estimation methods exhibit significant bias when outliers are endogenous and develop $L_0$-regularized estimation methods to overcome this issue. We propose systematic heuristic algorithms, notably an iterative hard-thresholding algorithm and a local combinatorial search refinement, to efficiently solve the combinatorial optimization problem of the \(L_0\)-regularized estimation. The properties of the \(L_0\) and \(L_1\)-regularized methods are examined through Monte Carlo simulations. We illustrate the practical value of our method with an empirical application to stock return forecasting.