2.1 Huber IV Estimator of the SAR Model

Kelejian and Prucha (Reference Kelejian and Prucha1998) propose the 2SLS estimator for the SAR model. Building upon their approach, we integrate Huber’s (Reference Huber1964) robust estimation technique to propose the HIV estimator in this section. The SAR model is formulated as

(2.1) $$ \begin{align} y_{i,n}=\lambda w_{i\cdot,n}Y_{n}+x_{i,n}'\beta+\epsilon_{i,n}\equiv\lambda\sum_{j=1}^{n}w_{ij,n}y_{j,n}+x_{i,n}'\beta+\epsilon_{i,n}, \end{align} $$

or in matrix form, $Y_{n}=\lambda W_{n}Y_{n}+X_{n}\beta +\epsilon _{n}$ , where $Y_{n}\equiv (y_{1,n},\ldots ,y_{n,n})'$ is the column vector of dependent variables, $W_{n}=(w_{ij,n})$ is an $n\times n$ non-stochastic spatial weights matrix and $w_{i\cdot ,n}$ is its $i^{th}$ row, $x_{i,n}=(x_{i1},\ldots ,x_{ip_{x}})'$ is a $p_{x}$ -dimensional column vector of exogenous regressors, including the intercept term $x_{i1}=1$ , $X_{n}=(x_{1,n},\ldots ,x_{n,n})'$ , and $\epsilon _{n}=(\epsilon _{1,n},\ldots ,\epsilon _{n,n})'$ is the column vector of disturbances. Denote the parameter vector $\theta =(\lambda ,\beta ')'$ and the true parameter vector ${\theta _{0}=(\lambda _{0},\beta _{0}')'}$ .

Huber (Reference Huber1964) introduces the following loss function:

$$\begin{align*}\rho(t)=\begin{cases} \begin{array}{@{}ccl} \frac{1}{2}t^{2} & & |t|\leq M\\ M|t|-\frac{1}{2}M^{2} & & |t|>M, \end{array} & \end{cases} \end{align*}$$

where M is a critical value. When $M=\infty $ , $\rho (\cdot )$ is the squared loss function, while as $M\to 0$ , $\rho (\cdot )$ approximates the absolute value loss function. Since $w_{i\cdot ,n}Y_{n}$ is an endogenous variable, the standard Huber M-estimator based on the loss function $\rho (\cdot )$ is typically inconsistent. However, the idea of moment estimation based on the derivative of $\rho (\cdot )$ can be adapted to the SAR model, provided that suitable IVs are available.

Notice that $\rho (t)$ is a continuously differentiable even function, and its derivative is

(2.2) $$ \begin{align} \psi(t)=\frac{d\rho(t)}{dt}=\begin{cases} \begin{array}{@{}ccl} M & & t>M\\ t & & |t|\leq M\\ -M & & t<-M. \end{array} & \end{cases} \end{align} $$

Suppose that

(2.3) $$ \begin{align} \operatorname{\mathrm{E}}[\psi(\epsilon_{i,n})|X_{n}]=0, \end{align} $$

which is a standard assumption in the Huber estimation literature, and $\epsilon _{i,n}$ ’s are independent conditional on $X_{n}$ . When the distribution of $\epsilon _{i,n}$ is symmetric about 0 given $x_{i,n}$ , Eq. (2.3) holds trivially. When $\epsilon _{i,n}$ ’s are i.i.d. but not symmetrically distributed conditional on $X_{n}$ , we can also shift the intercept such that Eq. (2.3) is satisfied,Footnote ⁴ and that is why an intercept term must be included. Notice that Eq. (2.3) allows conditional heteroskedasticity, so $\epsilon _{i,n}$ and $x_{i,n}$ might not be independent.

Since $W_{n}$ is nonstochastic, Eq. (2.3) implies

$$\begin{align*}\operatorname{\mathrm{E}}[\psi(\epsilon_{i,n})x_{i,n}]=0\text{ and }\operatorname{\mathrm{E}}[\psi(\epsilon_{i,n})w_{i\cdot,n}X_{n}]=0. \end{align*}$$

Hence, $x_{i,n}$ and $w_{i\cdot ,n}X_{n}$ can be regarded as IVs.Footnote ⁵ We can also add the second-order spatial effects $(w_{i\cdot ,n}W_{n}X_{n})$ into the IV vector $q_{i,n}$ , or use $G_{i\cdot ,n}X_{n}\beta _{0}$ as an IV for the endogenous variable $w_{i\cdot ,n}Y_{n}$ in the SAR model, where $G_{n}\equiv W_{n}(I_{n}-\lambda _{0}W_{n})^{-1}$ and $G_{i\cdot ,n}$ is its i-th row.Footnote ⁶ We can also censor these IVs such that they can satisfy the theory. Let $q_{i,n}\in \mathbb {R}^{p_{q}}$ be the column IV vector in the following and ${Q_{n}\equiv (q_{1,n},\ldots ,q_{n,n})'}$ denote the corresponding $n\times p_{q}$ IV matrix, where $p_{q}$ must be greater than or equal to $p_{x}+1$ , the dimension of $\theta $ . Then, we can write the moment conditions as

(2.4) $$ \begin{align} \operatorname{\mathrm{E}}[\psi(\epsilon_{i,n})q_{i,n}]=0. \end{align} $$

Based on Eq. (2.4) and the corresponding sample moments, we can derive an estimator for $\theta $ , which we term the HIV estimator. Denote $ {{\epsilon _{i,n}(\theta )\equiv }}y_{i,n}-\lambda w_{i\cdot ,n}Y_{n}-x_{i,n}'\beta $ , $\epsilon _{n}(\theta )\equiv (\epsilon _{1,n}(\theta ),\dots ,\epsilon _{n,n}(\theta ))'$ , and $\psi (\epsilon _{n}(\theta ))=(\psi (\epsilon _{1,n}(\theta )),\dots ,\psi (\epsilon _{n,n}(\theta )))'$ . The sample moments corresponding to Eq. (2.4) are $\frac {1}{n}Q_{n}'\psi (\epsilon _{n}(\theta ))=\frac {1}{n}\sum _{i=1}^{n}\psi (\epsilon _{i,n}(\theta ))q_{i,n}$ . Given a positive-definite GMM weighting matrix $\Omega _{n}$ , we can obtain an HIV estimator $\hat {\theta }_{n}=\arg \min _{\theta \in \Theta }\hat {T}_{n}(\theta )$ , where

(2.5) $$ \begin{align} \hat{T}_{n}(\theta)&=\left[\frac{1}{n}Q_{n}'\psi(\epsilon_{n}(\theta))\right]'\Omega_{n}\left[\frac{1}{n}Q_{n}'\psi(\epsilon_{n}(\theta))\right]\nonumber\\&=\left[\frac{1}{n}\sum_{i=1}^{n}\psi(\epsilon_{i,n}(\theta))q_{i,n}\right]'\Omega_{n}\left[\frac{1}{n}\sum_{i=1}^{n}\psi(\epsilon_{i,n}(\theta))q_{i,n}\right]. \end{align} $$

In practice, we can let $\Omega _{n}\equiv \left (\frac {1}{n}\sum _{i=1}^{n}q_{i,n}q_{i,n}'\right )^{-1}$ . In this case, when $M=\infty $ , the HIV estimator is reduced to the 2SLS estimator based on $q_{i,n}$ .

If conditional on $X_{n}$ , $\epsilon _{i,n}$ ’s are i.i.d., the best HIV (BHIV) for $w_{i\cdot ,n}Y_{n}$ is $G_{i\cdot ,n}X_{n}\beta _{0}$ , which is the same as that proposed in Lee (Reference Lee2003). To implement the BHIV estimation, we need a consistent estimator for $\theta _{0}$ and replace $G_{n}$ and $\beta _{0}$ by their estimates (see Section B of the Supplementary Material for details).

2.2 Huber GMM Estimator

Building upon the 2SLS in Kelejian and Prucha (Reference Kelejian and Prucha1998), the quadratic moment conditions for a pure SAR model in Kelejian and Prucha (Reference Kelejian and Prucha1999), and the BIV estimation in Lee (Reference Lee2003) for the SAR model, Lee (Reference Lee2007) proposes a GMM estimator for an SAR model that employs both linear and quadratic moments. This GMM estimator surpasses the 2SLS estimator in efficiency and can rival the MLE when the disturbances are normally distributed and appropriate quadratic moments are chosen. Lin and Lee (Reference Lin and Lee2010) extend Lee (Reference Lee2007) to allow heteroskedasticity.

Like a 2SLS estimator, the HIV estimator proposed in the previous section employs only the first-order moments based on the conditional IV moments motivated by Eq. (2.3). To enhance estimation efficiency, drawing inspiration from Kelejian and Prucha (Reference Kelejian and Prucha1999), Lee (Reference Lee2007), and Lin and Lee (Reference Lin and Lee2010), we explore “quadratic Huber moments” in this section.

Our discussion proceeds under the same assumptions outlined in Section 2.1, especially Eq. (2.3). Let $P_{kn}=(P_{ij,kn})$ ( $k=1,\ldots ,m$ ) be n-dimensional nonstochastic square matrices with zero diagonals, i.e., $P_{ii,kn}=0$ for all $i=1,\ldots ,n$ . Then, it follows from Eq. (2.3) and the independence of $\epsilon _{i,n}$ ’s that

$$\begin{align*}\operatorname{\mathrm{E}}[\psi(\epsilon_{n})'P_{kn}\psi(\epsilon_{n})]&=\sum_{i=1}^{n}\sum_{j\neq i}P_{ij,kn}\operatorname{\mathrm{E}}[\psi(\epsilon_{i,n})\psi(\epsilon_{j,n})]\nonumber\\&=\sum_{i=1}^{n}\sum_{j\neq i}P_{ij,kn}\operatorname{\mathrm{E}}[\psi(\epsilon_{i,n})]\operatorname{\mathrm{E}}\psi(\epsilon_{j,n})]=0. \end{align*}$$

Consequently, we can construct m sample quadratic Huber moments: $\frac {1}{n}\psi (\epsilon _{n}(\theta ))' P_{kn}\psi (\epsilon _{n}(\theta ))$ . As in Section 2.1, $q_{i,n}$ denotes a column IV vector and $Q_{n}\equiv (q_{1,n},\ldots ,q_{n,n})'$ . We construct a sample moment vector as in Lee (Reference Lee2007) and Lin and Lee (Reference Lin and Lee2010):

(2.6) $$ \begin{align} g_{n}(\theta)=\frac{1}{n}\left(\begin{array}{@{}c} \psi(\epsilon_{n}(\theta))'P_{1n}\psi(\epsilon_{n}(\theta))\\ \vdots\\ \psi(\epsilon_{n}(\theta))'P_{mn}\psi(\epsilon_{n}(\theta))\\ Q_{n}'\psi(\epsilon_{n}(\theta)) \end{array}\right). \end{align} $$

Given a GMM weighting matrix $\Omega _{n}$ , we obtain an HGMM estimator

$$\begin{align*}\hat{\theta}_{HGMM}=\arg\min_{\theta}g_{n}(\theta)'\Omega_{n}g_{n}(\theta). \end{align*}$$

2.3 Selection of the Critical Value M

The Huber loss function contains a critical value M. In Section 3, we show that the two proposed estimators are consistent and asymptotically normal for any fixed constant M. Nonetheless, in finite samples, it might be beneficial to employ a data-driven approach to select a suitable M that balances robustness and efficiency. Following Bickel (Reference Bickel1975), we propose to use

(2.7) $$ \begin{align} M=c\cdot\operatorname{\mathrm{median}}\left(|{{\epsilon}}-\operatorname{\mathrm{median}}({{\epsilon}})|\right)/0.6745, \end{align} $$

where c is a positive constant and $\operatorname {\mathrm {median}}\left (| { {\epsilon }}-\operatorname {\mathrm {median}}({{\epsilon }})|\right )/0.6745$ serves as a robust estimator for the standard deviation. The factor $0.6745$ corresponds to the 75% quantile of the standard normal distribution.

Huber (Reference Huber1964) and Maronna et al. (Reference Maronna, Martin, Yohai and Salibián-Barrera2019) demonstrate through simulation studies that the performance of Huber M estimators is relatively insensitive to the choice of c, when it lies between 1 and 1.5. Any value within the interval $[1,1.5]$ can yield satisfactory results. Both MATLAB and R default to a critical value of 1.345, which corresponds to a 95% relative asymptotic efficiency for location estimators under normality.Footnote ⁷ We also investigate this issue through simulations in Section 4. The conclusion is similar. When c is a constant between 1 and 1.5, our proposed estimators outperform traditional estimators for the SAR model under long-tailed disturbances or extreme outliers, with performances being comparable across this range. Thus, it makes little practical difference whether $c=1$ , 1.2, 1.345, or 1.5. Following the practice in MATLAB and R, we determine the critical value M using Eq. (2.7) with $c=1.345$ in both our simulation studies and empirical applications.

3.1 Spatial Near-Epoch Dependence

Since Kelejian and Prucha (Reference Kelejian and Prucha2001) establish a central limit theorem (CLT) for linear-quadratic forms in the context of SAR models, the theory of linear-quadratic forms has become a cornerstone in spatial econometrics (see Kelejian and Prucha, Reference Kelejian and Prucha2001; Lee, Reference Lee2004, Reference Lee2007). However, as $\psi (\epsilon _{i,n}(\theta ))$ does not conform to a linear-quadratic form, the theory of linear-quadratic forms is not applicable to Huber estimators, except when $\theta =\theta _{0}$ . Hence, we require certain weak spatial dependence concepts such that some laws of large numbers (LLN) can be utilized to establish asymptotic properties of our estimators. The theory of the spatial near-epoch dependence (NED) developed in Jenish and Prucha (Reference Jenish and Prucha2012) is prevalent in spatial econometrics and is capable of handling nonlinear models. Thus, using this approach, we establish large sample properties of our estimators via spatial NED theory.

Suppose that there are n individuals located in $D_{n}\subset \mathbb {R}^{d}$ and $|D_{n}|=n$ . As in Jenish and Prucha (Reference Jenish and Prucha2012), the distance between any two points $x,y\in \mathbb {R}^{d}$ is defined as $d_{xy}=\max _{i=1,\ldots ,d}|x_{i}-y_{i}|$ . We index these n individuals as $1,2,\ldots ,n$ . We identify their indices and locations in this article.Footnote ⁸

Next, we list the assumptions required to establish some moment conditions and spatial NED properties for $\{y_{i,n}:i\in D_{n}\}$ and related variables, which will be applied to establish the asymptotic theories in Sections 3.2 and 3.3.

Under Assumption 1, we consider the scenario of increasing domain asymptotics. This setting is widely presumed in spatial econometrics (Qu and Lee, Reference Qu and Lee2015; Xu and Lee, Reference Xu and Lee2015a, Reference Xu and Lee2015b, Reference Xu and Lee2018; Liu and Lee, Reference Liu and Lee2019), as it is the foundation of the spatial mixing and NED theory (Jenish and Prucha, Reference Jenish and Prucha2009, Reference Jenish and Prucha2012). This setting is particularly intuitive when individuals are situated in geographical locations like countries and cities. Nonetheless, social and economic characteristics can also be employed to define coordinates, thereby constructing distances between entities. Although Assumption 1 is popular in spatial econometrics, it rules out the case of infilled asymptotics, i.e., within a bounded space, the number of individuals tends to infinity. Infill asymptotics is popular in time series, especially functional time series (see, e.g., Phillips and Jiang, Reference Phillips and Jiang2025, among many other papers). In spatial statistics, it is also widely studied, e.g., the estimation of a Gaussian random field (see Cressie, Reference Cressie1993 and Gelfand et al., Reference Gelfand, Diggle, Guttorp and Fuentes2010 for details). Under infilled asymptotics, the spatial correlations of many observations might be strong and will not tend to zero. In this case, traditional estimation methods for the SAR model are generally inconsistent.Footnote ¹⁰

Assumption 2(i) imposes some moment conditions on the exogenous variable $x_{i,n}$ and the error term $\epsilon _{i,n}$ . Specifically, we require only the $\gamma $ -th-order moment for $x_{i,n}$ and $\operatorname {\mathrm {E}}|\epsilon _{i,n}|^{\gamma }$ for some $\gamma \geq 1$ . For the QMLE and the GMM estimator, when exogenous variables are uniformly bounded, Lee (Reference Lee2004, Reference Lee2007) requires disturbances to have a finite moment of order higher than 4, i.e., $\gamma>4$ in our notation. Assumption 2(ii) allows $x_{i,n}$ to be spatially correlated. Assumption 2(iii) has been discussed below Eq. (2.3).

Assumption 3 guarantees that the von Neumann series expansion for ${(I_{n}-\lambda W_{n})^{-1}}$ exists. This assumption is one of the critical assumptions leading to the NED of $\{y_{i,n}\}$ and $\{w_{i\cdot ,n}Y_{n}\}$ . It is also widely used in the spatial econometrics literature (see, e.g., the conditions in Lemma 2 in Kelejian and Prucha, Reference Kelejian and Prucha2010).

Assumption 4 allows two setups for the spatial weights matrix $W_{n}$ . Under Assumption 4(i), two individuals have direct interactions only if they are within a fixed finite distance. Under Assumption 4(ii), two individuals are allowed to have direct interactions even if they are far away from each other, but the strength of their interactions should decrease as distance increases. Assumptions 1 and 4 imply that $\sup _{n}||W_{n}||_{\infty }<\infty $ and $\sup _{n}||W_{n}||_{1}<\infty $ . With Assumption 3, they further imply that $\sup _{n}||(I_{n}-\lambda W_{n})^{-1}||_{\infty }<\infty $ . The uniform boundedness of the norms of $W_{n}$ and $(I_{n}-\lambda W_{n})^{-1}$ is prevalent in the literature of spatial econometrics (see, e.g., Assumption 3 in Kelejian and Prucha, Reference Kelejian and Prucha2010 and Assumption 5 in Lee, Reference Lee2004).

Assumptions 1–4 imply that $\{y_{i,n}\}$ and $\{w_{i\cdot ,n}Y_{n}\}$ are $L_{\gamma }$ -NED on $\{x_{i,n},\epsilon _{i,n}:i\in D_{n}\}$ (Lemma A.8). We summarize the moment and $L_{\gamma }$ -NED properties of the regressor $z_{i,n}\equiv (w_{i\cdot ,n}Y_{n},x_{i,n}')'$ , ${{\psi }({\epsilon _{i,n}(\theta ))}}$ , and common choices for IVs in Lemma A.8 in Appendix A.

3.2 Large Sample Properties of the Huber IV Estimator

In this section, we establish the large sample properties of the HIV estimator defined in Eq. (2.5).

In this article, we can use $w_{i\cdot ,n}X_{n}$ , $w_{i\cdot ,n}W_{n}X_{n}$ , or $G_{i\cdot ,n}X_{n}\beta _{0}$ as the IV for $w_{i\cdot ,n}Y_{n}$ . When $x_{i,n}$ ’s are uniformly bounded, Assumption 6 holds for $x_{i,n}$ , $w_{i\cdot ,n}X_{n}$ , $w_{i\cdot ,n}W_{n}X_{n}$ , and $G_{i\cdot ,n}X_{n}\beta _{0}$ ; if $x_{i,n}$ ’s are unbounded, we may censor them such that $q_{i,n}$ is uniformly bounded. We assume $\{q_{i,n}:i\in D_{n}\}$ to be uniformly bounded to simplify the proof.Footnote ¹¹ Alternatively, a more involved proof would permit unbounded $q_{i,n}$ .

Assumption 7 is the standard requirement for a GMM weighting matrix.

Assumption 8 is crucial for ensuring the identifiable uniqueness of $\theta _{0}$ . It implies that $\theta _{0}$ is the unique solution to $\frac {1}{n}\sum _{i=1}^{n}\operatorname {\mathrm {E}}\left [\psi (\epsilon _{i,n}(\theta ))q_{i,n}\right ]=0$ when n is large enough. However, Assumption 8 is slightly stronger than the condition that $\theta _{0}$ is the unique solution to $\frac {1}{n}\sum _{i=1}^{n}\operatorname {\mathrm {E}}\left [\psi (\epsilon _{i,n}(\theta ))q_{i,n}\right ]=0$ . Assumption 8 excludes the possibility that $0<\left \Vert \frac {1}{n}\sum _{i=1}^{n}\operatorname {\mathrm {E}}\left [\psi (\epsilon _{i,n}(\theta ))q_{i,n}\right ]\right \Vert \to 0$ as $n\to \infty $ .

Under the above assumptions, the HIV estimator is consistent, as stated in Theorem 1.

Next, we derive the asymptotic distribution of the HIV estimator $\hat {\theta }$ . Given that the objective function $\hat {T}_{n}(\theta )$ is not differentiable with respect to $\theta $ , the conventional Taylor expansion cannot be employed to derive the asymptotic distribution of $\hat {\theta }$ . However, we can utilize the empirical process theory for derivation. We generalize Theorem 3.3 in Pakes and Pollard (Reference Pakes and Pollard1989), which studies the asymptotic distributions of GMM estimators under non-differentiable moment conditions, to accommodate heterogeneous data in Proposition B.2. The critical condition in Proposition B.2 is a condition similar to stochastic equicontinuity (SEC), specifically condition (iii) there. Within the empirical process literature, several results exist to validate the SEC condition of an empirical process. However, most of them are tailored to independent or stationary data. Data generated by spatial econometric models are usually both heterogeneous and spatially correlated. Therefore, we need to generalize those results to allow for both heterogeneity and spatial correlation. Pollard (Reference Pollard1985, p. 305) extends the result in Huber (Reference Huber1967) to establish a set of sufficient conditions for SEC in the context of independent data based on a bracketing condition. In this article, we generalize the result in Pollard (Reference Pollard1985) to account for serial or spatial correlation (see Proposition B.1 for details).

To establish the asymptotic distribution of $\hat {\theta }$ , we need some additional regularity conditions. $\dot {\psi }(\epsilon )\equiv 1(|\epsilon |\leq M)$ equals the derivative of $\psi (\epsilon )$ almost everywhere. Define $\Gamma _{n}(\theta )=\frac {1}{n}\sum _{i=1}^{n}\operatorname {\mathrm {E}}[\dot {\psi }(y_{i,n}-z_{i,n}'\theta )q_{i,n}z_{i,n}']$ , $\Gamma _{n}=\Gamma _{n}(\theta _{0})$ , and $V_{n}\equiv \operatorname {\mathrm {var}}\left \{ \frac {1}{n}\sum _{i=1}^{n}\psi (\epsilon _{i,n})q_{i,n}\right \} $ . By the independence of $\epsilon _{i,n}$ ’s and the law of iterated expectation, straightforward calculations yield that

$$\begin{align*}V_{n}=\frac{1}{n}\sum_{i=1}^{n}\operatorname{\mathrm{E}}[\psi(\epsilon_{i,n})^{2}q_{i,n}q_{i,n}']. \end{align*}$$

Assumption 10. (i) The class of functions $\{\Gamma _{n}(\theta ):\,n\in \mathbb {N}\}$ is equicontinuous at $\theta _{0}$ ; (ii) there exists a positive constant $\delta _{0}$ such that

$$\begin{align*}\liminf_{n\to\infty}\inf_{\theta_{1},\dots,\theta_{l}\in B(\theta_{0},\delta_{0})}\min\operatorname{\mathrm{sv}}\left(\left[\begin{array}{@{}c} \Gamma_{1n}(\theta_{1})\\ \vdots\\ \Gamma_{ln}(\theta_{l}) \end{array}\right]\right)>0. \end{align*}$$

$\Gamma _{in}(\theta )$ is the $i^{th}$ row of $\Gamma _{n}(\theta )$ ; and (iii) $\lim _{n\to \infty }\Gamma _{n}=\Gamma _{0}$ and $\Gamma _{0}$ has a full column rank.

Assumption 9 is a standard condition for establishing asymptotic distributions. Assumption 10 imposes some conditions on $\Gamma _{n}(\theta )$ and $\Gamma _{n}$ such that condition (ii) of Proposition B.2 holds. We need different $\theta _{i}$ for different rows because there does not exist a mean-value theorem for vector-valued functions. Under standard conditions, $\Gamma _{n}(\theta )$ is continuous at $\theta _{0}$ . The equicontinuity at $\theta _{0}$ for $\Gamma _{n}(\theta )$ means that as the sample size n increases to infinity, $\Gamma _{n}(\theta )$ does not change too much around $\theta _{0}$ . This is a reasonable requirement; otherwise, it would be infeasible to conduct any reliable inference based on finite samples.

Theorem 2. Under Assumptions 1–11, if the $\alpha $ in Assumption 4 satisfies $\alpha>\left [\frac {2}{\min (2,\gamma )}+1\right ]d$ ,

$$\begin{align*}\sqrt{n}(\hat{\theta}-\theta_{0})\xrightarrow{d}N(0,(\Gamma_{0}'\Omega\Gamma_{0})^{-1}\Gamma_{0}'\Omega V_{0}\Omega\Gamma_{0}(\Gamma_{0}'\Omega\Gamma_{0})^{-1}). \end{align*}$$

The optimal choice for $\Omega $ is $\Omega \propto V_{0}^{-1}$ , and with such an $\Omega $ , $\sqrt {n}(\hat {\theta }-\theta _{0})\xrightarrow {d}N(0,(\Gamma _{0}'V_{0}^{-1}\Gamma _{0})^{-1})$ .

A standard two-step GMM procedure can be applied to obtain an HIV estimator based on the optimal GMM weighting matrix, where $\Omega _{n}\equiv \left [\frac {1}{n}\sum _{i=1}^{n}\psi (\hat {\epsilon }_{i,n})^{2} q_{i,n}q_{i,n}'\right ]^{-1}$ for the optimal GMM estimation, where $\hat {\epsilon }_{i,n}\equiv y_{i,n}-z_{i,n}'\tilde {\theta }$ and $\tilde {\theta }$ is an initial consistent estimator. We can estimate $\Gamma _{0}$ and $V_{0}$ by their sample means $\frac {1}{n}\sum _{i=1}^{n}\dot {\psi }(\hat {\epsilon }_{i,n})q_{i,n}z_{i,n}'$ and $\frac {1}{n}\sum _{i=1}^{n}\psi (\hat {\epsilon }_{i,n})^{2}q_{i,n}q_{i,n}'$ , respectively (see Theorem 5 for their consistency).

3.3 Large Sample Properties of the Huber GMM Estimator

For the HGMM estimator, we require the linear Huber moments to satisfy the conditions for the HIV estimators. In addition, some conditions for the quadratic Huber moments are also needed. We use $W_{n}$ , $W_{n}^{2}-\operatorname {\mathrm {diag}}(W_{n}^{2})$ , $S_{n}^{-1}-\operatorname {\mathrm {diag}}(S_{n}^{-1})$ , or $G_{n}-\operatorname {\mathrm {diag}}(G_{n})$ as $P_{kn}$ in this article. Since

(3.1) $$ \begin{align} \psi(\epsilon_{n}(\theta))'P_{kn}\psi(\epsilon_{n}(\theta))=\sum_{i=1}^{n}\psi(\epsilon_{i,n}(\theta))\left[P_{i\cdot,kn}\psi(\epsilon_{n}(\theta))\right], \end{align} $$

if $\left \{ P_{i\cdot ,kn}\psi (\epsilon _{n}(\theta ))\right \} _{i\in D_{n}}$ is uniformly bounded and uniformly $L_{2}$ -NED on $\{x_{i,n}\,\epsilon _{i,n}:i\in D_{n}\}$ , by Lemma A.5, the product term $\left \{ \psi (\epsilon _{i,n}(\theta ))\left [P_{i\cdot ,kn}\psi (\epsilon _{n}(\theta ))\right ]\right \} _{i\in D_{n}}$ is also uniformly bounded and uniformly $L_{2}$ -NED. Then, we can investigate the limiting properties of the quadratic Huber moments using the NED theory. We show that $\left \{ P_{i\cdot ,kn}\psi (\epsilon _{n}(\theta ))\right \} _{i\in D_{n}}$ is uniformly $L_{2}$ -NED on $\{x_{i,n}\,\epsilon _{i,n}:i\in D_{n}\}$ in Assumption A.9. Thus, we are ready to establish asymptotic properties of the HGMM estimator.

Next, we obtain the asymptotic distribution of the HGMM estimator $\hat {\theta }_{HGMM}$ . To do so, we need some regularity conditions for $\Gamma _{n}(\theta )\equiv \frac {\partial g_{n}(\theta )}{\partial \theta '}$ , $\Gamma _{n}\equiv \Gamma _{n}(\theta _{0})$ , and $V_{n}\equiv n\operatorname {\mathrm {var}}\left (g_{n}(\theta _{0})\right )$ . We collect the expressions of these matrices in Section C of the Supplementary Material. Notice that the $V_{n}$ in Eq. (C.2) is a block diagonal matrix. This is because a linear form and a quadratic form (with all diagonal elements being zero) are uncorrelated, as suggested by Eq. (3.2) in Kelejian and Prucha (Reference Kelejian and Prucha2001).

Theorem 4. Under Assumptions 1–9, 12, and 13, if the $\alpha $ in Assumption 4(ii) satisfies $\alpha>\left [\frac {2}{\min (2,\gamma )}+1\right ]d$ , then

$$\begin{align*}\sqrt{n}(\hat{\theta}_{HGMM}-\theta_{0})\xrightarrow{d}N(0,(\Gamma_{0}'\Omega\Gamma_{0})^{-1}\Gamma_{0}'\Omega V_{0}\Omega\Gamma_{0}(\Gamma_{0}'\Omega\Gamma_{0})^{-1}). \end{align*}$$

When $\Omega \propto V_{0}^{-1}$ , $\sqrt {n}(\hat {\theta }_{HGMM}-\theta _{0})\xrightarrow {d}N(0,(\Gamma _{0}'V_{0}^{-1}\Gamma _{0})^{-1})$ .

The condition $\alpha>\left [\frac {2}{\min (2,\gamma )}+1\right ]d$ ensures that the NED coefficient satisfies the condition of Lemma A.3 in Jenish and Prucha (Reference Jenish and Prucha2012). If $|w_{ij,n}|$ decreases exponentially as $d_{ij}$ increases, this condition becomes unnecessary.

Following Lin and Lee (Reference Lin and Lee2010), we estimate $V_{0}$ in Eq. (C.2) in the Supplementary Material by replacing $\Sigma _{n}\equiv \operatorname {\mathrm {var}}(\psi (\epsilon _{n}))$ by $\hat {\Sigma }_{n}\equiv \operatorname {\mathrm {diag}}(\psi (\hat {\epsilon }_{1,n})^{2},\ldots ,\psi (\hat {\epsilon }_{n,n})^{2})$ and replacing $\Sigma _{n}^{d}$ by $\hat {\Sigma }_{n}^{d}\equiv (\psi (\hat {\epsilon }_{1,n})^{2},\ldots ,\psi (\hat {\epsilon }_{n,n})^{2})'$ , where $\hat {\epsilon }_{i,n}=y_{i,n}-z_{i,n}'\hat {\theta }$ is the residual, i.e., the estimator for $V_{0}$ is

Similarly, we estimate $\Gamma _{n}$ in Eq. (C.1) in the Supplementary Material by

Then, we can estimate the asymptotic variance of the optimal GMM estimator by $(\hat {\Gamma }_{n}'\hat {V}^{-1}\hat {\Gamma }_{n})^{-1}$ .

The conditions $\sup _{i,n,\epsilon ,x}|f_{in}(\epsilon |x)|<\infty $ and $\sup _{i,n,\epsilon ,x}|\epsilon f_{in}(\epsilon |x)|<\infty $ are satisfied by most continuous and piecewise continuous p.d.f.s in practice. As long as $f_{in}(x)$ is continuous and $\lim _{|x|\to \infty }|xf_{in}(x)|=0$ , Assumption 2(ii) holds. Even when $f_{in}(\cdot )$ is not continuous, Assumption 2(ii) is still satisfied in many cases, e.g., the p.d.f. of an exponential distribution or a uniform distribution.

3.4 Finite Sample Breakdown Point Analysis

Breakdown point analysis is a key criterion for assessing the robustness of a statistical procedure. “Informally, a breakdown of a statistical procedure means that the procedure no longer conveys useful information on the data-generating mechanism” (Genton and Lucas, Reference Genton and Lucas2003). The concept of breakdown point varies across different literature. Traditionally, the breakdown point refers to the fraction of contamination (or outliers) that results in unbounded estimators or estimators on the boundary of the parameter space (Hampel, Reference Hampel1971). This definition is well-suited for independent data, but Genton and Lucas (Reference Genton and Lucas2003) argue that it is not readily applicable to dependent data. Thus, they extend the definition in Hampel (Reference Hampel1971) to accommodate dependent data. Genton, (Reference Genton, Dutter, Filzmoser, Gather and Rousseeuw2003, p. 151) finds that when the spatial weights matrix is lower triangular, the breakdown point for the MLE of a pure SAR model without exogenous variables is zero; when $w_{ij}=1$ if and only if $|i-j|=1$ , the asymptotic breakdown point of the least median of squares estimator for the pure SAR model with additive outliers is 14.4%.

Since “the breakdown point is most useful in small sample situations” (Huber and Ronchetti, Reference Huber and Ronchetti2009, p.279), we focus on finite sample breakdown point. As in Genton and Lucas (Reference Genton and Lucas2003), we consider innovation outliers, i.e., outliers in the error terms $\epsilon _{i,n}$ .Footnote ¹² Unlike additive outliers used in Genton (Reference Genton, Dutter, Filzmoser, Gather and Rousseeuw2003), the impact of innovation outliers can be propagated through the spatial weights matrix $W_{n}$ to all dependent variables $y_{j,n}$ ’s. Assuming the presence of contamination $\zeta \bar {\epsilon }_{n}$ , where $\zeta \in [0,\infty ]$ controls the magnitude of the contamination, the model then becomes

$$\begin{align*}Y_{n}=\lambda_{0}W_{n}Y_{n}+X_{n}\beta_{0}+(\epsilon_{n}+\zeta\bar{\epsilon}_{n}). \end{align*}$$

So, $Y_{n}=S_{n}^{-1}(X_{n}\beta +\epsilon _{n}+\zeta \bar {\epsilon }_{n})$ . Suppose all entries of $W_{n}$ are nonnegative and $\lambda _{0}>0$ to simplify our analysis. We consider two cases:

(i) The network implied by $W_{n}$ is bidirectional and connected, i.e., $W_{n}$ cannot be written as a block diagonal matrix. In this case, for any $i\neq j$ , there exists a $k\in \mathbb {N}$ such that $\left (W_{n}^{k}\right )_{ij}\neq 0$ . Consequently, all entries of $S_{n}^{-1}=\sum _{l=0}^{\infty }\lambda _{0}^{l}W_{n}^{l}$ are greater than 0. Therefore, even if there exists only one outlier, e.g., $\bar {\epsilon }_{n}=(1,0,\ldots ,0)'$ , as $\zeta \to +\infty $ , all $y_{i,n}$ ’s will tend to $+\infty $ . In this case, it is impossible to obtain any reasonable estimators. Thus, the breakdown point for innovations outliers is zero. This outcome is independent of any specific estimation method and is attributed to the network’s connectedness. This conclusion parallels the zero breakdown point for the SAR model discussed in Genton (Reference Genton, Dutter, Filzmoser, Gather and Rousseeuw2003, p. 151).

(ii) $W_{n}$ is a block diagonal matrix. Suppose the network can be separated into $G>1$ subnetworks, represented by $W_{g}$ ( $g=1,\ldots ,G$ ). These subnetworks are not linked, but each of them is a bidirectional and connected network, i.e., $W_{n}=\operatorname {\mathrm {diag}}(W_{1},\ldots ,W_{G})$ . In each subnetwork $W_{g}$ ( $g=1,\ldots ,G$ ), if we introduce one extreme outlier to one $\epsilon _{i,n}$ , where i is an individual in group g, i.e., the contaminated error is $\epsilon _{i,n}+\zeta $ . Since $W_{g}$ is connected, by an argument similar to case (i), as $\zeta \to +\infty $ , all $y_{j,n}$ ’s in group g will tend to $+\infty $ . Therefore, the finite sample breakdown point must be $\leq \frac {G-1}{n}$ . Consequently, if the number of groups G is not large, the finite sample breakdown point is quite small. When G is large, the data are more similar to independent data, and it is reasonable to have a greater breakdown point.

In summary, due to network structures, when extreme contamination occurs in innovations, it will also lead to extreme values for other observations. As a result, the breakdown point is smaller than that in the typical case of independent data.

3.5 Influence Function

The influence function is another useful criterion for measuring the robustness of an estimator. In this section, we will demonstrate that, under regularity conditions, the influence functions of both the HIV estimator and the HGMM estimator are bounded.

Consider the influence functions of the HIV estimator first. Suppose that with probability $1-\delta $ , $(Y_{n},X_{n})$ is generated by an uncontaminated model, and with probability $\delta $ , $(Y_{n},X_{n})=(Y_{n}^{0},X_{n}^{0})$ , where $(Y_{n}^{0},X_{n}^{0})$ is a nonstochastic outlier. With such contaminated data, the pseudo true value of the parameter $\theta $ is defined as

(3.2) $$ \begin{align} & \theta_{\delta}=\arg\min\left\{ \frac{1-\delta}{n}\sum_{i=1}^{n}\operatorname{\mathrm{E}}\left[q_{i,n}\psi(y_{i,n}-z_{i,n}'\theta)\right]+\frac{\delta}{n}\sum_{i=1}^{n}\left(\operatorname{\mathrm{E}} q_{i,n}\right)\psi(y_{i,n}^{0}-(z_{i,n}^{0})'\theta)\right\} ' \nonumber\\& \Omega_{n}\left\{ \frac{1-\delta}{n}\sum_{i=1}^{n}\operatorname{\mathrm{E}}\left[q_{i,n}\psi(y_{i,n}-z_{i,n}'\theta)\right]+\frac{\delta}{n}\sum_{i=1}^{n}\left(\operatorname{\mathrm{E}} q_{i,n}\right)\psi(y_{i,n}^{0}-(z_{i,n}^{0})'\theta)\right\} , \end{align} $$

where the expectation $\operatorname {\mathrm {E}}$ is taken for the uncontaminated data and $z_{i,n}^{0}\equiv (w_{i\cdot ,n}Y_{n}^{0},(x_{i,n}^{0})')'$ . The influence function is defined as

$$\begin{align*}\lim_{\delta\to0}\frac{\theta_{\delta}-\theta_{0}}{\delta}. \end{align*}$$

Assume that $\lim _{\delta \to 0}\theta _{\delta }=\theta _{0}$ , which is a standard condition for calculating influence functions, and that $\psi (y_{i,n}^{0}-(z_{i,n}^{0})'\theta )$ is differentiable around $\theta _{0}$ . The first-order condition of Eq. (3.2) is

$$\begin{align*}\begin{aligned} & \left\{ \frac{1-\delta}{n}\sum_{i=1}^{n}\operatorname{\mathrm{E}}\left[\dot{\psi}(y_{i,n}-z_{i,n}\theta_{\delta})q_{i,n}z_{i,n}'\right]-\frac{\delta}{n}\sum_{i=1}^{n}\operatorname{\mathrm{E}}\left[q_{i,n}(z_{i,n}^{0})'\right]\dot{\psi}(y_{i,n}^{0}-(z_{i,n}^{0})'\theta_{\delta})\right\} '\cdot\\ & \Omega_{n}\left\{ \frac{1-\delta}{n}\sum_{i=1}^{n}\operatorname{\mathrm{E}}\left[q_{i,n}\psi(y_{i,n}-z_{i,n}\theta_{\delta})\right]+\frac{\delta}{n}\sum_{i=1}^{n}\left(\operatorname{\mathrm{E}} q_{i,n}\right)\psi(y_{i,n}^{0}-(z_{i,n}^{0})'\theta_{\delta})\right\} =0. \end{aligned} \end{align*}$$

Doing Taylor expansion for the second term at $\theta _{0}$ in the above equation, we obtain

$$\begin{align*}\begin{aligned} & \left\{ \frac{1-\delta}{n}\sum_{i=1}^{n}\operatorname{\mathrm{E}}\left[\dot{\psi}(y_{i,n}-z_{i,n}\theta_{\delta})q_{i,n}z_{i,n}'\right]-\frac{\delta}{n}\sum_{i=1}^{n}\operatorname{\mathrm{E}}\left[q_{i,n}(z_{i,n}^{0})'\right]\psi(y_{i,n}^{0}-(z_{i,n}^{0})'\theta_{\delta})\right\} '\Omega_{n}\cdot\\ & \left\{ \frac{1-\delta}{n}\sum_{i=1}^{n}\left(\operatorname{\mathrm{E}}\left[\psi(\epsilon_{i,n})q_{i,n}z_{i,n}'\right]+\operatorname{\mathrm{E}}\left[\dot{\psi}(y_{i,n}-z_{i,n}\bar{\theta})q_{i,n}z_{i,n}'\right](\theta_{0}-\theta_{\delta})\right)\right.\nonumber\\&\qquad\left.+\frac{\delta}{n}\sum_{i=1}^{n}\left(\operatorname{\mathrm{E}} q_{i,n}\right)\psi(y_{i,n}^{0}-(z_{i,n}^{0})'\theta_{\delta})\right\} =0, \end{aligned} \end{align*}$$

where $\bar {\theta }$ is between $\theta _{0}$ and $\theta _{\delta }$ .Footnote ¹³ Hence, as $\delta \to 0$ ,

$$\begin{align*}\begin{aligned}\frac{\theta_{\delta}-\theta_{0}}{\delta}\to & \left\{ \left(\frac{1}{n}\sum_{i=1}^{n}\operatorname{\mathrm{E}}\left[\dot{\psi}(\epsilon_{i,n})z_{i,n}q_{i,n}'\right]\right)\Omega_{n}\left(\frac{1}{n}\sum_{i=1}^{n}\operatorname{\mathrm{E}}\left[\dot{\psi}(\epsilon_{i,n})q_{i,n}z_{i,n}'\right]\right)\right\} ^{-1}\cdot\\ & \left(\frac{1}{n}\sum_{i=1}^{n}\operatorname{\mathrm{E}}\left[\dot{\psi}(\epsilon_{i,n})z_{i,n}q_{i,n}'\right]\right)\Omega_{n}\left[\frac{1}{n}\sum_{i=1}^{n}\left(\operatorname{\mathrm{E}} q_{i,n}\right)\psi(y_{i,n}^{0}-(z_{i,n}^{0})'\theta_{0})\right]. \end{aligned} \end{align*}$$

It follows from $|\psi (\cdot )|\leq M$ that $\frac {\theta _{\delta }-\theta _{0}}{\delta }$ is bounded as $\delta \to 0$ . In addition, under Assumptions 7 and 10, because $\sup _{i,n}\operatorname {\mathrm {E}}\left \Vert q_{i,n}\right \Vert <\infty $ , we have $\limsup _{n\to \infty }\operatorname {\mathrm {E}}\lim _{\delta \to 0}\left \Vert \frac {\theta _{\delta }-\theta _{0}}{\delta }\right \Vert <\infty $ . In other words, the bound for the influence function of the HIV estimator for each n remains bounded as $n\to \infty $ .

By similar arguments, the influence function of the HGMM estimator is also bounded.

4.1 The Effects of Different c in M

In Section 2.3, we propose to use Eq. (2.7) to determine the critical value M: $M=c\cdot \frac {\operatorname {\mathrm {median}}\left (|{{\epsilon }}-\operatorname {\mathrm {median}}({{\epsilon }})|\right )}{0.6745}$ , where $c=1.345$ , the default critical value in both R and MATLAB. We first conduct some simulations to examine the performance of the HIV estimator and the HGMM estimator under different choices of c.

We generate independent disturbances $\epsilon _{i,n}$ ’s from the following seven distributions: (1) i.i.d. $N(0,1)$ ; (2) a mixed normal distribution $0.95N(0,1)+0.05N(0,25)$ ; (3) a mixed distribution of $N(0,1)$ and the standard Cauchy distribution $0.95N(0,1)+0.05t_{1}$ ; (4) $N(0,1)$ contaminated by a single extreme error term: $\epsilon _{i}\sim N(0,1)$ , $i=2,3,\ldots ,n$ and $\epsilon _{1}\sim N(150,1)$ ; (5) $N(0,1)$ contaminated by two extreme error terms: $\epsilon _{i}\sim N(0,1)$ , $i=3,4,\ldots ,n$ and $\epsilon _{1},\epsilon _{2}\sim N(150,1)$ ; (6) Laplace distribution times $x_{i2,n}$ ; and (7) $\exp (1)$ times $x_{i2,n}^{2}$ . Setups (2)–(7) are contaminated or long-tailed distributions. Setups (2) and (3) belong to mixture distributions, and setups (4) and (5) are contaminated by extreme outliers.Footnote ¹⁵ Setups (6) and (7) are long-tailed and subject to conditional heteroskedasticity. Estimation results are summarized in Tables E.1–E.3 in the Supplementary Material.

For the HIV estimator, we consider a single tuning parameter c under the several different types of disturbances with a sample size of $n=100.$ As the default value of c is $1.345$ in R and MATLAB, we take this value as a benchmark to calculate the relative mean squared error (RelaMSE) for each tuning parameter value, i.e., $RelaMSE=\frac {MSE(\hat {\theta }_{c})}{MSE(\hat {\theta }_{1.345})}$ . The results are shown in Table E.1 in the Supplementary Material. In most cases, for contaminated or long-tailed distributions, the HIV estimator performs better when $c\in [1,1.5]$ than when $c>1.5.$ Within the range $[1,1.5]$ , the relative efficiency varies by less than 10% across different tuning parameters, with the exception of the case involving Cauchy contamination. As a comparison, from Table 3 in Section 4.2, the relative efficiency gain of the HIV estimator over 2SLS is at least 86% (in most cases more than 100%) when $n=100$ . Therefore, the efficiency difference due to varying c within [1, 1.5] is minor, compared to the efficiency gained by using Huber estimation. Notably, $c=1.345$ works substantially better than other values of c in the presence of Cauchy contamination. Hence, we use $c=1.345$ for all Huber IV estimators in subsequent experiments.

For the HGMM estimator, we adopt two tuning parameters: $c_{1}$ for quadratic moments and $c_{2}$ for linear moments, where $c_{1},c_{2}\in \{1,1.2,1.345,1.5,1.8,2,2.5,3\}$ . We explore the impact of using different pairs of tuning parameters in the following three scenarios: (1) disturbances drawn from $0.95N(0,1)+0.05N(0,25)$ , (2) $N(0,1)$ contaminated by a single outlier from $N(150,1)$ , and (3) Laplace distribution times $x_{i2,n}$ . Similarly, we use the pair $(1.345,1.345)$ as a benchmark to calculate the relative MSE for each tuning parameter combination, i.e., $RelaMSE=\frac {MSE(\hat {\theta }_{(c_{1},c_{2})})}{MSE(\hat {\theta }_{(1.345,1.345)})}$ . The results regarding $\lambda $ and $\beta _{2}$ are presented in Tables E.2 and E.3 in the Supplementary Material. In general, the relative MSEs of both $\hat {\lambda }$ and $\hat {\beta }_{2}$ tend to get larger when $c_{2}\geq 1.8$ and do not vary much for $c_{2}\in [1,1.5]$ . For a fixed $c_{2}$ , especially when $c_{2}\in [1,1.5]$ , the influence of $c_{1}$ on the relative MSE is rather small, with variations typically within 5% in most cases. In cases where there is no clear pattern, such as the relative RMSE of $\hat {\lambda }$ in scenario (2), $(c_{1},c_{2})=(1.345,1.345)$ dominates most other tuning parameter combinations. Thus, Eq. (2.7) serves as a reasonable criterion for selecting M with $(c_{1},c_{2})=(1.345,1.345)$ , and there is no need to differentiate $c_{1}$ and $c_{2}$ .

In conclusion, for both the HIV estimator and the HGMM estimator, our simulation results indicate that $c=1.345$ is a satisfactory choice for a variety of scenarios, including contaminated data, long-tailed data, and data with conditional heteroskedasticity. Thus, we maintain $c=1.345$ in our subsequent simulation studies.

4.2 Compare the Performance of Various Estimators

Next, we compare the performances of the QMLE, the 2SLS estimator in Kelejian and Prucha (Reference Kelejian and Prucha1998), the BIV estimator in Lee (Reference Lee2003), the GMM estimators in Lin and Lee (Reference Lin and Lee2010), the quantile IV (QIV) estimator in Su and Yang (Reference Su and Yang2011), the 2SGM estimator in Wagenvoort and Waldmann (Reference Wagenvoort and Waldmann2002),Footnote ¹⁶ the Huber maximum likelihood estimator (HMLE) in Tho et al. (Reference Tho, Ding, Hui, Welsh and Zou2023), and the three estimators proposed in our article: the HIV estimator, the BHIV estimator, and HGMM estimator. For both the GMM estimator and the HGMM estimator, we adopt $P_{n}=\tilde {G}_{n}-\operatorname {\mathrm {diag}}(\tilde {G}_{n})$ for the quadratic moments, where $\tilde {G}_{n}=W_{n}(I_{n}-\tilde {\lambda }W_{n})^{-1}$ , and $\tilde {\lambda }$ is a consistent estimator for $\lambda _{0}$ .Footnote ¹⁷ And we use $X_{n}$ and $G_{n}X_{n}\tilde {\beta }$ as the IVs, where $\tilde {\beta }$ is a consistent estimator for the true parameter value $\beta _{0}$ . For Huber type estimators, the critical value of M is determined as in Section 2.3 with $c=1.345$ .

We have tried three different sample sizes in the experiments: $n=100,200,$ and 400. For each setting, simulations are repeated 1,000 times to obtain these estimates’ sample bias and RMSE. We report the results on parameters $\lambda $ and $\beta _{2}$ under sample sizes $n=100$ and $n=400$ in Tables 1 and 2, respectively. Full results—including all parameters (except for the intercept coefficient $\beta _{1}$ ) and all sample sizes ( $n=100, 200, 400$ )—are provided in Tables E.4–E.6 in the Supplementary Material. Under asymmetrically distributed disturbances, the true intercept $\beta _{10}$ depends on the critical value M, and is different from the true intercept under the moment condition $\operatorname {\mathrm {E}}\epsilon _{i,n}=0$ , so it is troublesome to compare its estimates. Thus, we omit them from all reported tables—including the full results in Tables E.4–E.6 in the Supplementary Material. To obtain a better insight into the merits of robust estimators, we present the relative efficiency of each Huber-type estimator $\hat {\theta }$ compared to its corresponding traditional estimator, e.g., HIV/2SLS, BHIV/BIV, and HGMM/GMM, in Table 3. The relative efficiency of $\hat {\theta }_{1}$ relative to $\hat {\theta }_{2}$ is calculated as $MSE(\hat {\theta }_{2})/MSE(\hat {\theta }_{1})$ .

Table 1 Estimation results ( $n=100$ )

1: iid $N(0,1)$ ; 2: Mixed normal: iid $0.95N(0,1)+0.05N(0,25)$ ; 3: NC: $0.95N(0,1)+0.05t_{1}$ ; 4: Extreme contamination: $\epsilon _{1}\sim N(150,1)$ ; 5: Extreme contamination: $\epsilon _{1},\epsilon _{2}\sim N(150,1)$ ; 6: Laplace $\cdot x_{i2,n}$ : Laplace distribution scaled by $x_{i2}$ ; and 7: $\exp (1)\cdot x_{i2,n}^{2}$ : $\exp (1)$ scaled by $x_{i2,n}^{2}$ . BIV: Best IV; GMM: GMM in Lin and Lee (Reference Lin and Lee2010); QIV: quantile IV in Su and Yang (Reference Su and Yang2011); 2SGM: Two-stage generalized M-estimator in Wagenvoort and Waldmann (Reference Wagenvoort and Waldmann2002); HMLE: Huber MLE in Tho et al. (Reference Tho, Ding, Hui, Welsh and Zou2023); HIV: Huber IV; BHIV: Best Huber IV; and HGMM: Huber GMM.

Table 2 Estimation results ( $n=400$ )

1: iid $N(0,1)$ ; 2: Mixed normal: iid $0.95N(0,1)+0.05N(0,25)$ ; 3: NC: $0.95N(0,1)+0.05t_{1}$ ; 4: Extreme contamination: $\epsilon _{1}\sim N(150,1)$ ; 5: Extreme contamination: $\epsilon _{1},\epsilon _{2}\sim N(150,1)$ ; 6: Laplace $\cdot x_{i2,n}$ : Laplace distribution scaled by $x_{i2}$ ; and 7: $\exp (1)\cdot x_{i2,n}^{2}$ : $\exp (1)$ scaled by $x_{i2,n}^{2}$ . BIV: Best IV; GMM: GMM in Lin and Lee (Reference Lin and Lee2010); QIV: quantile IV in Su and Yang (Reference Su and Yang2011); 2SGM: Two-stage generalized M-estimator in Wagenvoort and Waldmann (Reference Wagenvoort and Waldmann2002); HMLE: Huber MLE in Tho et al. (Reference Tho, Ding, Hui, Welsh and Zou2023); HIV: Huber IV; BHIV: Best Huber IV; and HGMM: Huber GMM.

Table 3 Relative efficiency of Huber estimators with respect to traditional estimators

1: iid $N(0,1)$ ; 2: Mixed normal: iid $0.95N(0,1)+0.05N(0,25)$ ; 3: NC: iid $0.95N(0,1)+0.05t_{1}$ ; 4: Extreme contamination: $\epsilon _{1}\sim N(150,1)$ ; 5: Extreme contamination: $\epsilon _{1},\epsilon _{2}\sim N(150,1)$ ; 6: Laplace $\cdot x_{i2,n}$ : Laplace distribution scaled by $x_{i2}$ ; and 7: $\exp (1)\cdot x_{i2,n}^{2}$ : $\exp (1)$ scaled by $x_{i2,n}^{2}$ . BIV: Best IV; GMM: GMM in Lin and Lee (Reference Lin and Lee2010): HIV: Huber IV; BHIV: Best Huber IV; HGMM: Huber GMM. Relative efficiency= $MSE(\hat {\theta }_{traditional})/MSE(\hat {\theta }_{Huber})$ .

Below are some interesting observations obtained from the simulation results:

(1) For all sample sizes and distributions, the biases and RMSEs of the three proposed Huber estimators are notably small. As the sample size n increases, the RMSEs of our estimators decrease at roughly the rate of $n^{-1/2}$ . For instance, under the mixed normal distribution, when $n=100$ , 200, and 400, the RMSE for $\hat {\lambda }_{HGMM}$ is 0.0661, 0.0473, and 0.0378, respectively. These results conform with the prediction of our large sample theories.

(2) When the disturbances are i.i.d. normally distributed,

$$\begin{align*}QMLE &\succsim GMM\succsim HGMM\succsim HMLE\succ2SLS\approx BIV\\ & \approx HIV\approx BHIV\succ2SGM\succ QIV, \end{align*}$$

where $\succ $ means “performing better than,” $\succsim $ means “performing slightly better than,” and $\approx $ means “performing similarly.” Although QMLE performs the best, the differences between the HGMM estimator, the GMM estimator, and the QMLE are minor. Among all robust estimators, the HGMM estimator exhibits the smallest efficiency loss under $N(0,1)$ . For all sample sizes and all parameters,

$$\begin{align*}\frac{RMSE(HGMM)-RMSE(QMLE)}{RMSE(QMLE)}<2\%, \end{align*}$$

$$\begin{align*}\frac{RMSE(HGMM)-RMSE(GMM)}{RMSE(GMM)}<2\%. \end{align*}$$

From Table 3, under normally distributed disturbances, Huber-type estimators exhibit only a slight loss in efficiency (in most cases $<6\%$ ) compared to traditional estimators. This modest efficiency trade-off is to be expected, as robustness often comes at the cost of some efficiency. The efficiency losses incurred by our estimators are minimal. The aforementioned observations can be clearly seen from Figures 1 and 2. Please note that each pair of bars sharing similar colors in the figure represents a traditional estimator and its Huber variant, respectively (for instance, 2SLS and HIV, BIV and BHIV, and GMM and HGMM estimators), while the purple bars represent 2SGM estimator and QIV estimator.

Figure 1 RMSE of $\hat {\lambda }$ of various estimators under $N(0,1)$ .

Figure 2 RMSE of $\hat {\beta }_{2}$ of various estimators under $N(0,1)$ .

(3) In extreme contamination scenarios (cases 4 and 5), the biases and RMSEs of the traditional estimators, including QMLE, 2SLS, BIV, and GMM, are extremely large, among which BIV is the worst (when $n=100$ , for $\lambda $ , the efficiency of the BHIV estimator compared to the BIV estimator is $3.2\times 10^{5}$ ). In contrast, all robust estimators exhibit significantly superior performance. The Huber MLE is slightly better than HGMM, but the difference is minimal. Both significantly outperform other robust estimators. Overall, we summarize performances of these estimators in most (although not all) cases as follows:

$$\begin{align*}HMLE &\succsim HGMM\succ BHIV\succ HIV\succ QIV\succ2SGM\succ QMLE\\ &\succ2SLS\approx GMM\succ BIV. \end{align*}$$

The conclusions for setups (2) and (3) are similar.

(4) In the presence of conditional heteroskedasticity (cases 6 and 7), the HGMM estimator and QIV estimator outperform other robust estimators, which in turn perform better than traditional estimators. We can see this more clearly from Figures 3 and 4. Results can be summarized as follows:

$$\begin{align*}HGMM &\approx QIV\succ BHIV\succ HIV\approx HMLE\succ2SGM\succ GMM\\ & \approx QMLE\succ2SLS\succ BIV. \end{align*}$$

Figure 3 RMSE of $\hat {\lambda }$ of various estimators under $Laplace\cdot x_{i2}$ .

Figure 4 RMSE of $\hat {\beta }_{2}$ of various estimators under $Laplace\cdot x_{i2}$ .

(5) In all setups (2)–(7), which involve contaminated, long-tailed, or conditionally heteroskedastic disturbances, Huber-type estimators consistently outperform their traditional counterparts:

$$\begin{align*}\begin{cases} \begin{array}{@{}cc} HGMM\succ GMM, & \\ BHIV\succ BIV, & \\ HIV\succ2SLS. & \end{array}\end{cases} \end{align*}$$

(6) All the proposed estimators are computationally efficient. For example, when the errors follow a mixed normal distribution and the sample size $n=100$ , our desktop computer (CPU: Intel(R) Xeon(R) E-2136; RAM: 32G; hard drive: 1T; and software: MATLAB) takes an average of 0.013 seconds to run the HGMM estimation and 0.25 seconds for the Huber MLE by Tho et al. (Reference Tho, Ding, Hui, Welsh and Zou2023). With $n=400$ , the HGMM estimation takes 0.047 seconds, while the Huber MLE takes 5 seconds. On average, our estimators can be computed 20 times faster than the estimator in Tho et al. (Reference Tho, Ding, Hui, Welsh and Zou2023) for $n=100$ , and 100 times faster for $n=400$ .

In conclusion, the simulation outcomes demonstrate that our estimators are robust against both contaminated or long-tailed disturbances and conditional heteroskedasticity, incur minimal efficiency loss under short-tailed disturbances, and are computationally convenient. The Huber MLE by Tho et al. (Reference Tho, Ding, Hui, Welsh and Zou2023) works well under contaminated disturbances, but is less robust to conditional heteroskedasticity, and computationally slower. QIV estimator excels under heteroskedastic long-tailed disturbances, but performs worse than the HGMM estimator and the Huber MLE in other scenarios. Overall, HGMM estimator performs the best among all these estimators in simulations.

4.3 A Test Experiment

We further carry out simulation studies to evaluate the empirical size and the power of testing $H_{0}:\beta _{20}=\beta _{30}$ vs. $H_{0}:\beta _{20}\neq {{\beta _{30}}}$ at level 5%, as suggested by an anonymous referee. The experimental setup closely follows that of Section 4.3, except that $\beta _{30}\in \{1,-1\}$ . Detailed results are reported in Tables E.7 and E.8 in the Supplementary Material. We observe that the empirical size is around 0.05 and the empirical power is close to 1 for all settings and all sample sizes. As the sample size increases, the empirical size gets closer to 0.05 and the empirical power tends toward 1. This shows that our asymptotic distribution approximates the true distribution of the estimators well.