A multi-step kernel–based regression estimator that adapts to error distributions of unknown form

Abstract

For linear regression models, we propose and study a multi-step kernel density-based estimator that is adaptive to unknown error distributions. We establish asymptotic normality and almost sure convergence. An efficient EM algorithm is provided to implement the proposed estimator. We also compare its finite sample performance with five other adaptive estimators in an extensive Monte Carlo study of eight error distributions. Our method generally attains high mean-square-error efficiency. An empirical example illustrates the gain in efficiency of the new adaptive method when making statistical inference about the slope parameters in three linear regressions.

Keywords:

• Adaptive estimation
• EM algorithm
• kernel density estimate
• least squares estimate
• linear regression

MATHEMATICS SUBJECT CLASSIFICATION:

• 62F35
• 62J05

1. Introduction

In parametric linear regression analysis one often imposes the model assumptions that the errors are independent and normally distributed. The normality assumption is convenient, as it is well known that the maximum likelihood estimator (MLE) of the unknown parameter vector simplifies to the least squares estimator (LSE). Naturally, an invalid assumption on the error distribution F comes at a cost; the MLE is in general neither consistent nor asymptotically efficient under model misspecification. Moreover, in practice, it can lead to inaccurate or invalid statistical inference; see Sec. 5. This has motivated the search for alternative, (semi)parametric, estimators that retain asymptotic efficiency when F is unknown.

One approach is adaptive estimation, which “adapts” to an unknown, or incorrectly specified, distribution F by maximizing an estimated likelihood function based on an initial estimate of the error distribution; see Bickel (1982), Linton and Xiao (2007), Yuan and De Gooijer (2007), and the references therein. The adaptive idea has been studied for (non)linear regression models using non- and semiparametric methods to estimate F or its probability density function (pdf) f.

There are various alternative adaptive estimation methods for non- and semiparametric regression problems with errors of unknown distributional form. For instance, the empirical likelihood method of Owen (2001) has been used to obtain adaptive confidence limits and likelihood ratio test statistics for regression parameters; see also Owen (1988, 1990, 1991), Qin and Lawless (1994), Kitamura (1997), Kitamura (2007), among many others. Another example is the multivariate adaptive regression splines (MARS) (Friedman 1991) which is a global adaptive nonparametric method for fitting nonlinear regression models. PolyMARS (Kooperberg, Bose, and Stone 1997) is an extension of MARS that allows for multiple polychotomous regression. Time series MARS, or TSMARS, can be used for nonlinear time series analysis and forecasting; see, e.g., De Gooijer (2017). More recently, Wang and Yao (2012) proposed a minimum average variance estimation method, a dimension reduction technique, which can be adaptive to different error distributions. Also, Chen, Wang, and Yao (2015) developed an adaptive estimation method for varying coefficient models.

Recently, Yao and Zhao (2013) proposed an adaptive kernel density-based estimator for classical linear regression models, called KDRE. In particular, with an estimate of the “true” parameter vector, f is modeled by a kernel density-based estimator of the regression residuals. In the second step, parameter estimates are obtained by maximizing a local, kernel-based, log-likelihood function using the first-step estimated density function as the true one. Through a simulation study, Yao and Zhao (2013) show that the resulting KDRE is asymptotically equivalent to the oracle estimator in which the true error pdf is known.

Now, it is well known that each LSE-based residual is the sum of two components: one is the true error, the other is a linear function of the entire vector of errors. Since, in finite samples, the second term will tend to be normally distributed (as long as the errors have finite variance), the residuals for small samples will appear more normal than would the unobserved values of the errors. This tendency is called supernormality; see Bassett and Koenker (1982), Bloomfield (1974), and White and MacDonald (1980). Hence, for the two-step KDRE method, it is likely that the finite-sample properties of the proposed estimator strongly depend on the empirical distribution of the residuals, more closely resembling normality than would be the case by using the pdf of the errors itself, if this were in fact possible. This makes the KDRE method nonoptimal.

In this paper, we remedy this deficiency of the KDRE method by further iteration. That is, we maximize over different kernel-based likelihood functions. These different likelihood functions follow iteratively from parameter estimates that result from maximizing the (previous) likelihood function. The algorithm iterates until the parameter estimates (and, hence, the estimated likelihood function) reaches a fixed point. The new estimator is called multi-step KDRE (M-KDRE). In finite samples, one may expect that M-KDRE yields better estimation results; see Robinson (1988). Actually, from an extensive Monte Carlo study where we compare the finite-sample performance of M-KDRE with estimates based on five other, parametric and (semi)nonparametric adaptive estimators (AEs), we find that iterating over the likelihood functions can indeed strongly increase finite-sample performance. In fact, we find that M-KDRE outperforms all other considered AEs. In an empirical example based on Andrabi, Das, and Khwaja (2017), we show that LSE estimates may be misleading. M-KDRE outperforms the other considered estimators in terms of out-of-sample prediction performance. Furthermore, M-KDRE provides strong evidence that the treatment effect as described in Andrabi, Das, and Khwaja (2017) does not exist.

Theoretically, we establish strong (almost sure) convergence to the true parameter vector, under relative weak conditions. We also show that the M-KDRE method is adaptive, i.e., asymptotically normal and efficient. Furthermore, its computation is made convenient by proposing an EM type algorithm.

The rest of the paper is organized as follows. Section 2 introduces the new adaptive M-KDRE method and contains our theoretical results. An efficient EM type algorithm to implement the proposed estimator is also given in this section. In Section 3, we describe and explain five alternative adaptive estimation methods. Section 4 contains results of a simulation-based study of the finite sample properties of the M-KDRE method and compares it with those of the adaptive estimation methods discussed in Section 3. In Section 5, we present the empirical application of our method to the educational data set as used in Andrabi, Das, and Khwaja (2017). Section 6 gives a summary and some concluding remarks. Proofs are presented in Appendix A.

2. Multi-Step kernel Density-Based regression estimation

2.1. Model and method

Consider the general linear regression problem with observations (1) $$y_i={\mathbf{x}}_i^{\prime }{\mathit{\beta }}_0+{\epsilon }_i,\hspace{1em}(i=1,\dots ,n)$$(1) where y_i is a univariate response variable, ${\mathbf{x}}_i=(x_{i,1},\dots ,x_{i,p})\prime$ is a p-dimensional (p < n) vector of covariates, and ${\mathit{\beta }}_0\in \mathcal{B}\subseteq {\mathbb{R}}^p$ is an unknown parameter vector including an intercept. Here $(y_i,{\mathbf{x}}_i,{\epsilon }_i)$ are independent and identically distributed (i.i.d.) realizations from a common random source $(y,\mathbf{x},\mathit{\epsilon }).$ Moreover, the ${\epsilon }_i$‘s are assumed to have some common unknown pdf $f(\epsilon ),$ and $E[{\epsilon }_i|{\mathbf{x}}_i,{\mathit{\beta }}_0]=0$ and $E[\left|\epsilon \right||{\mathbf{x}}_i,{\mathit{\beta }}_0]<\infty$ $(i=1,\dots ,n).$ Model (1(1) $$y_i={\mathbf{x}}_i^{\prime }{\mathit{\beta }}_0+{\epsilon }_i,\hspace{1em}(i=1,\dots ,n)$$(1) ) is semiparametric with $\mathit{\beta }$ and $f(·)$ its parametric and non-parametric part, respectively.

Let ${\widehat{\mathit{\beta }}}_{\text{LSE}}$ be the LSE of ${\mathit{\beta }}_0$ in (1(1) $$y_i={\mathbf{x}}_i^{\prime }{\mathit{\beta }}_0+{\epsilon }_i,\hspace{1em}(i=1,\dots ,n)$$(1) ), which is a natural estimator to start the M-KDRE method. Also, let ${\widehat{\mathit{\beta }}}^{(u)}$ denote an estimator of ${\mathit{\beta }}_0$ at iteration step $u=0,1,\dots .$ Then, with the conditions introduced above, the proposed M-KDRE can be obtained as follows.

1. Initial step: At u = 0, start with ${\widehat{\mathit{\beta }}}^{(0)}={\widehat{\mathit{\beta }}}_{\text{LSE}}.$ Compute the residuals ${\widehat{\epsilon }}_i^{(0)}=y_i-{\mathbf{x}}_i^{\prime }{\widehat{\mathit{\beta }}}^{(0)}$ $(i=1,\dots ,n).$

2. Compute the Rosenblatt-Parzen kernel-based estimator ${\widehat{f}}_n^{(u)}(·)$ of $f(·).$ That is (2) $${\widehat{f}}_n^{(u)}(x)={(n{h}_n)}^{-1}\sum _{i=1}^{n}K(\frac{x-{\widehat{\epsilon }}_i^{(u-1)}}{h_n})$$(2) where h_n > 0 is the bandwidth.

3. Let $c(\mathit{\beta })=n^{-1}{\sum }_{i=1}^n(y_i-{\mathbf{x}}_i^{\prime }\mathit{\beta }).$ Then, using (2(2) $${\widehat{f}}_n^{(u)}(x)={(n{h}_n)}^{-1}\sum _{i=1}^{n}K(\frac{x-{\widehat{\epsilon }}_i^{(u-1)}}{h_n})$$(2) ), compute (3) $${\widehat{\mathit{\beta }}}^{(u)}=\arg {\max }_{\mathit{\beta }\in \mathcal{B}}{\widehat{Q}}_u(\mathit{\beta })\hspace{1em}\mathrm{s}.\mathrm{t}. \mathrm{c}({\mathit{\beta }}^{(\mathrm{u})})=0$$(3) where ${\widehat{Q}}_u(\mathit{\beta })$ is the local log-likelihood function (4) $${\widehat{Q}}_u(\mathit{\beta })=\sum _{i=1}^n\hspace{0.17em}\text{ln}\hspace{0.17em}{\widehat{f}}_n^{(u)}(y_i-{\mathbf{x}}_i^{\prime }\mathit{\beta })$$(4)

4. Repeat steps (ii)–(iii) until convergence at iteration step $(u+1)$ $(u=1,2,\dots ).$

It is worth mentioning that Yao and Zhao (2013) only iterate over the empirical likelihood function of the first initial, LSE, estimate while in the above case we allow for stochastic fluctuations in ${\widehat{f}}_n^{(u)}(·).$ As we show, this generalization will increase mean-square error parameter efficiency.

2.2. Asymptotic properties

In this section we give the asymptotic properties of the uth M-KDRE, denoted hereafter by ${\widehat{\mathit{\beta }}}^{(u)}.$ For convenience, when we emphasize the dependence on $\mathit{\beta },$ we use $f(y_i|\mathit{\beta })$ to denote $f(y_i-{\mathbf{x}}_i^{\prime }\mathit{\beta }).$ Technical details, lemmas, and proofs are given in Appendix A.

Theorem 2.1.

(Almost sure (a.s.) convergence) Under the assumptions of Lemma A.1 and

1. $f(y_i|\widehat{\mathit{\beta }})$ is non-parametrically identifiable,

2. $\mathit{\beta }\in \mathcal{B}$ where $\mathcal{B}\subseteq {\mathbb{R}}^p$ is compact,

3. $f(y_i|\mathit{\beta })$ is continuous for each $\mathit{\beta }\in \mathcal{B},$

4. $E[\underset{\mathit{\beta }\in \mathcal{B}}{\sup }|\hspace{0.17em}\text{ln}\hspace{0.17em}f(y_i|\mathit{\beta })|]<\infty ,$

   then

(5) $${\widehat{\mathit{\beta }}}^{(u)}\stackrel{a.s.}{\to }{\mathit{\beta }}_0$$(5)

The following notation is used throughout the next part of the paper. Let $L_n(y_i|\mathit{\beta })=\hspace{0.17em}\text{ln}\hspace{0.17em}f(y_i|\mathit{\beta })=\hspace{0.17em}\text{ln}\hspace{0.17em}f(y_i-{\mathbf{x}}_i^{\prime }\mathit{\beta }).$ Then ${\mathbf{d}}_{\beta }(\mathit{\beta })=\partial L_n(y_i|\mathit{\beta })/\partial \mathit{\beta }=-f\prime (y_i|\mathit{\beta })f^{-1}(y_i|\mathit{\beta }){\mathbf{x}}_i,{\mathbf{d}}_{\beta \beta }(\mathit{\beta })=\partial L_n^2(y_i|\mathit{\beta })/\partial \mathit{\beta }\partial {\mathit{\beta }}^{\prime }\mathrm{}=\mathrm{}(f(y_i|\mathit{\beta })f″(y_i|\mathit{\beta })-f{\prime }^2(y_i|\mathit{\beta })/f^2(y_i|\mathit{\beta })){\mathbf{x}}_i{\mathbf{x}}_i^{\prime },$ where $f\prime (u)\mathrm{}=\mathrm{}\partial f(u)/\partial u,$ $f″(u)=\partial f\prime (u)/\partial u,$ and $f{\prime }^2(u)=f\prime (u)f\prime (u).$ Given these notations, the Fisher information matrix (FIM) for the unconstrained linear regression problem, evaluated at ${\mathit{\beta }}_0,$ can be defined as (6) $${\mathcal{I}}_{\beta \beta }=E({\mathbf{d}}_{\beta }({\mathit{\beta }}_0){\mathbf{d}}_{\beta }({\mathit{\beta }}_0)\prime )=-E({\mathbf{d}}_{\beta \beta }({\mathit{\beta }}_0))$$(6) where the second equality holds under mild regularity conditions. If ${\mathcal{I}}_{\beta \beta }$ is nonsingular, then ${\mathcal{I}}_{\beta \beta }^{-1}$ is the unconstrained Cramér-Rao bound (CRB) for the mean-square error (MSE) covariance matrix of any unbiased estimator of ${\mathit{\beta }}_0.$

Observe that incorporation of the one-dimensional linear constraint $c({\mathit{\beta }}^{(u)})=0$ in step (iii) of the M-KDRE method leads to a p-dimensional parameter vector that has only p − 1 independent components. As a consequence, the FIM is singular and the CRB may not be an informative lower bound on the MSE matrix of the resulting estimator. So the asymptotic distribution of ${\widehat{\mathit{\beta }}}^{(u)}$ degenerates. For deterministic linear parameter constraints, Stoica and Ng (1998) formulated a constrained CRB (CCRB) that explicitly incorporates the active constraint information with the original FIM, singular or nonsingular. Their general setting is for q (q < p) continuously differentiable constraints $g(\mathit{\beta })=0.$ Assuming $\mathit{\beta }$ is regular in the active set of linear constraints, the q × p gradient matrix $\mathbf{G}=\partial g(\mathit{\beta })/\partial {\mathit{\beta }}^{\prime }$ has full row rank q, with G independent of $\mathit{\beta }.$ Hence, there exists a matrix $\mathbf{U}\in {\mathbb{R}}^{p\times (p-q)}$ whose p-dimensional columns form an orthonormal null space of the range space of the row vectors in G, i.e., such that (7) $$\mathbf{GU}=0 \text{and} \mathbf{U}\prime \mathbf{U}={\mathbf{I}}_{p-q}$$(7) where ${\mathbf{I}}_{p-q}$ denotes the identity matrix of size p – q. For nonlinear deterministic constraints, G and U are functions of $\mathit{\beta };$ see, e.g., Moore, Sadler, and Kozick (2008).

Recently, Ren et al. (2015) extended the deterministic CCRB of Stoica and Ng (1998) to a hybrid parameter vector with both nonrandom and random parameter constraints. In the case of the M-KDRE method the constraint is not deterministic, depending on the random variables ${\mathbf{x}}_i.$ Then the matrix U in (7(7) $$\mathbf{GU}=0 \text{and} \mathbf{U}\prime \mathbf{U}={\mathbf{I}}_{p-q}$$(7) ) depends on the estimate $\widehat{\mathit{\beta }}$ of ${\mathit{\beta }}_0,$ say $\mathbf{U}(\widehat{\mathit{\beta }}).$ The resulting CCRB, as a special case of the hybrid CCRB of Ren et al. (2015), can be stated as follows.

Theorem 2.2.

Let $\widehat{\mathit{\beta }}$ be an unbiased estimate of ${\mathit{\beta }}_0$ satisfying the active functional constraints $g(\mathit{\beta })=\mathbf{0}$, and let $\mathbf{U}=\mathbf{U}(\widehat{\mathit{\beta }})$ be defined in (7)(7) $$\mathbf{GU}=0 \text{and} \mathbf{U}\prime \mathbf{U}={\mathbf{I}}_{p-q}$$(7) . Then, under certain regularity conditions, if $\mathbf{U}\prime {\mathcal{I}}_{\beta \beta }\mathbf{U}$ is nonsingular, $$E((\widehat{\mathit{\beta }}-{\mathit{\beta }}_0)(\widehat{\mathit{\beta }}-{\mathit{\beta }}_0)\prime )\ge \mathbf{U}{(\mathbf{U}\prime {\mathcal{I}}_{\beta \beta }\mathbf{U})}^{-1}\mathbf{U}\prime$$ where the equality is achieved if and only if $\widehat{\mathit{\beta }}-{\mathit{\beta }}_0=\mathbf{U}{(\mathbf{U}\prime {\mathcal{I}}_{\beta \beta }\mathbf{U})}^{-1}\mathbf{U}\prime {\mathbf{d}}_{\beta }(\mathit{\beta })$, in the mean- square sense.

Remark 1.

Note that rather than requiring a nonsingular FIM ${\mathcal{I}}_{\beta \beta },$ the alternative condition is that $\mathbf{U}\prime {\mathcal{I}}_{\beta \beta }\mathbf{U}$ is nonsingular. Thus, the unconstrained FIM may be singular, or, equivalently, the unconstrained model unidentifiable, but the constrained model must be identifiable, at least locally. Ren et al. (2015) show that the difference between the standard CRB-based covariance matrix and the CCRB-based covariance matrix is a positive semi-definite matrix. This result is expected since the presence of parameter constraints can be considered as additional information to improve the performance of the estimator under study.

Theorem 2.3.

(Normality and efficiency) If model (1)(1) $$y_i={\mathbf{x}}_i^{\prime }{\mathit{\beta }}_0+{\epsilon }_i,\hspace{1em}(i=1,\dots ,n)$$(1) holds, ${\left\{{\epsilon }_i\right\}}_{i=1}^n$ are i.i.d. with unknown density f(x) where f is a continuous function symmetric around zero with bounded continuous derivatives that satisfy

1. $\int xf(x)\mathrm{d}x=0,$

2. $E[{\left(\frac{\partial \hspace{0.17em}\text{ln}\hspace{0.17em}f(x)}{\partial x}\right)}^2+|\frac{{\partial }^2\hspace{0.17em}\text{ln}\hspace{0.17em}f(x)}{\partial x^2}|+|\frac{{\partial }^3\hspace{0.17em}\text{ln}\hspace{0.17em}f(x)}{\partial x^3}\left|\right]<\infty ,$

   ${\left\{{\mathbf{x}}_i\right\}}_{i=1}^{\infty }$ satisfies,

3. $\exists \mathrm{ }\mathrm{ }0<M<\infty$ such that $\Vert \mathbf{x}\Vert <M,$

   $K(·)$ is a symmetric and four times continuously differentiable function such that

4. $\exists \mathrm{ }0<\rho <\infty$ such that K(x) = 0 $\forall x:\Vert \mathbf{x}\Vert \ge \rho$

   holds, and

5. when $n\to \infty ,n{h}_n^4\to \infty$ and $n{h}_n^8\to 0,$

then ${\widehat{\mathit{\beta }}}^{(u)}$ $(u=1,2,\dots )$ is asymptotically normal and efficient. That is, as $n\to \infty ,$ (8) $$\sqrt{n}({\widehat{\mathit{\beta }}}^{(u)}-{\mathit{\beta }}_0)\stackrel{d}{\to }N(0,\mathbf{U}{(\mathbf{U}\prime {\mathcal{I}}_{\beta \beta }\mathbf{U})}^{-1}\mathbf{U}\prime )$$(8)

Remark 2.

All conditions are practical and easy to satisfy. Condition (ii) is used to guarantee the adaptiveness of ${\widehat{\mathit{\beta }}}^{(u)}.$

2.3. EM algorithm

In this section, we propose an EM type algorithm by noticing that (4(4) $${\widehat{Q}}_u(\mathit{\beta })=\sum _{i=1}^n\hspace{0.17em}\text{ln}\hspace{0.17em}{\widehat{f}}_n^{(u)}(y_i-{\mathbf{x}}_i^{\prime }\mathit{\beta })$$(4) ) has a mixture log-likelihood form with an imposed constraint. Specifically, given an initial parameter estimate ${\widehat{\mathit{\beta }}}^{(0)}$ and the set of initial estimates ${\left\{{\widehat{\epsilon }}_i^{(0)}\right\}}_{i=1}^n,$ the $(k+1)$ th iteration of the EM algorithm to maximize (4(4) $${\widehat{Q}}_u(\mathit{\beta })=\sum _{i=1}^n\hspace{0.17em}\text{ln}\hspace{0.17em}{\widehat{f}}_n^{(u)}(y_i-{\mathbf{x}}_i^{\prime }\mathit{\beta })$$(4) ) (the uth likelihood function) is as follows.

E-step: Calculate the classification probabilities, (9) $$p_{ij,(k+1)}^{(u)}=\frac{\hspace{0.17em}\text{exp}\hspace{0.17em}\{-\frac{1}{2h_n^2}{(y_i-{\mathbf{x}}_i^{\prime }{\widehat{\mathit{\beta }}}_{(k)}^{(u)}-{\widehat{\epsilon }}_j^{(u-1)})}^2\}}{\sum _{\mathcal{\ell }=1}^n\hspace{0.17em}\text{exp}\hspace{0.17em}\{-\frac{1}{2h_n^2}{(y_i-{\mathbf{x}}_i^{\prime }{\widehat{\mathit{\beta }}}_{(k)}^{(u)}-{\widehat{\epsilon }}_{\mathcal{\ell }}^{(u-1)})}^2\}}$$(9)

M-step: Update ${\widehat{\mathit{\beta }}}_{(k)}^{(u)}$ with (10) $$\begin{array}{c}{\widehat{\mathit{\beta }}}_{(k+1)}^{(u)}={\widehat{\mathit{\beta }}}_{\text{LSE}}-{\left(\sum _{i=1}^n{\mathbf{x}}_i{\mathbf{x}}_i^{\prime }\right)}^{-1}\sum _{i=1}^n\left({\mathbf{x}}_i\sum _{j=1}^n{p}_{ij,(k+1)}^{(u)}{\widehat{\epsilon }}_j^{(u-1)}\right)\\ \hspace{1em}+\frac{1}{n}{\left(\sum _{i=1}^n{\mathbf{x}}_i{\mathbf{x}}_i^{\prime }\right)}^{-1}\sum _{i=1}^n{\mathbf{x}}_i\left(\sum _{i=1}^n\sum _{j=1}^n{p}_{ij,(k+1)}^{(u)}{\widehat{\epsilon }}_j^{(u-1)}\right)\end{array}$$(10) where (10(10) $$\begin{array}{c}{\widehat{\mathit{\beta }}}_{(k+1)}^{(u)}={\widehat{\mathit{\beta }}}_{\text{LSE}}-{\left(\sum _{i=1}^n{\mathbf{x}}_i{\mathbf{x}}_i^{\prime }\right)}^{-1}\sum _{i=1}^n\left({\mathbf{x}}_i\sum _{j=1}^n{p}_{ij,(k+1)}^{(u)}{\widehat{\epsilon }}_j^{(u-1)}\right)\\ \hspace{1em}+\frac{1}{n}{\left(\sum _{i=1}^n{\mathbf{x}}_i{\mathbf{x}}_i^{\prime }\right)}^{-1}\sum _{i=1}^n{\mathbf{x}}_i\left(\sum _{i=1}^n\sum _{j=1}^n{p}_{ij,(k+1)}^{(u)}{\widehat{\epsilon }}_j^{(u-1)}\right)\end{array}$$(10) ) follows from using a Gaussian kernel for density estimation. The choice of the kernel is not critical. Any symmetric kernel can be used for our method. However, the Gaussian second order kernel provides an explicit solution of the EM algorithm.

Theorem 2.4.

Under the linearity constraint $c(\mathit{\beta })=0$, each iteration of the above E- and M-steps will monotonically increase the local log-likelihood (4(4) $${\widehat{Q}}_u(\mathit{\beta })=\sum _{i=1}^n\hspace{0.17em}\text{ln}\hspace{0.17em}{\widehat{f}}_n^{(u)}(y_i-{\mathbf{x}}_i^{\prime }\mathit{\beta })$$(4) ), i.e., ${\widehat{Q}}_n({\widehat{\mathit{\beta }}}_{(k+1))}^{(u)})\ge {\widehat{Q}}_n({\widehat{\mathit{\beta }}}_{(k)}^{(u)})$, for all k.

Remark 3.

For the EM type algorithm, we use a full-kernel method rather than a leave-one-out method as in Yao and Zhao (2013). The approach has the following advantage. If a certain residual ${\widehat{\epsilon }}_j^{(u-1)}$ is extremely large, $p_{ij,(k+1)}^{(u)}$ will be close to zero for all $i\ne j$ and close to one for i = j. This implies that the effect of the residual is limited to the observation for which the following iteration of $\mathit{\beta }$ is likely to lead to a residual that is similar in magnitude. Hence, the effect of a large residual on the maximization of (4(4) $${\widehat{Q}}_u(\mathit{\beta })=\sum _{i=1}^n\hspace{0.17em}\text{ln}\hspace{0.17em}{\widehat{f}}_n^{(u)}(y_i-{\mathbf{x}}_i^{\prime }\mathit{\beta })$$(4) ) is small. In the leave-one-out method, the effect of the residual may be considerably larger as $p_{ij,(k+1)}^{(u)}$ is likely to have a substantial value for several observations.

Remark 4.

The EM type algorithm is considered converged when $\max |{\widehat{\mathit{\beta }}}_{(k)}^{(u)}-{\widehat{\mathit{\beta }}}_{(k+1)}^{(u)}|$ is smaller than a threshold value, where $\max \left|\mathbf{A}\right|$ denotes the largest (absolute) element in A. In the uth step of the M-KDRE method, the EM algorithm is initialized by the estimate at the $(u-1)$ th step. That is, ${\widehat{\mathit{\beta }}}_{(0)}^{(u)}={\widehat{\mathit{\beta }}}^{(u-1)}.$

3. Some alternative adaptive estimation methods

3.1. SBS method

Stone (1975), Bickel (1982), and Schick (1993) (henceforth SBS) introduce a two-step AE. Let $\tilde{\mathit{\beta }}$ be a certain $\sqrt{n}$-consistent estimator of ${\mathit{\beta }}_0.$ Then, an infeasible two-step estimator can be defined as $$\widehat{\mathit{\beta }}=\tilde{\mathit{\beta }}+n^{-1}{\mathcal{I}}_{\beta \beta }^{-1}(\tilde{\mathit{\beta }}){\mathbf{d}}_{\beta }(\tilde{\mathit{\beta }})$$ where ${\mathcal{I}}_{\beta \beta }(\tilde{\mathit{\beta }})$ is Fisher’s information matrix evaluated at $\tilde{\mathit{\beta }}$ and ${\mathbf{d}}_{\beta }(\tilde{\mathit{\beta }})$ is a corresponding $p\times 1$ score vector. The infeasibility of $\widehat{\mathit{\beta }}$ follows from the fact that f is unknown, and hence ${\mathcal{I}}_{\beta \beta }$ and ${\mathbf{d}}_{\beta }$ are unknown. The approach of SBS is to replace ${\mathbf{d}}_{\beta }(\tilde{\mathit{\beta }})$ by ${\widehat{\mathbf{d}}}_{\beta }(\tilde{\mathit{\beta }})=-{\sum }_{i=1}^n\widehat{f}{\prime }_n(y_i|\tilde{\mathit{\beta }}){\widehat{f}}_n^{-1}(y_i|\tilde{\mathit{\beta }}){\mathbf{x}}_i$ where ${\widehat{f}}_n(x)$ is defined in a similar way as in (2(2) $${\widehat{f}}_n^{(u)}(x)={(n{h}_n)}^{-1}\sum _{i=1}^{n}K(\frac{x-{\widehat{\epsilon }}_i^{(u-1)}}{h_n})$$(2) ) and $\widehat{f}{\prime }_n(x)$ is its derivative with respect to x. Similarly, ${\mathcal{I}}_{\beta \beta }^{-1}(\tilde{\mathit{\beta }})$ is replaced with $n^2{[{\sum }_{i=1}^n{\mathbf{x}}_i{\mathbf{x}}_i^{\prime }\sum ({\widehat{f}}_n^{\prime }(y_i|\tilde{\mathit{\beta }}){\widehat{f}}_n^{-1}(y_i|\tilde{\mathit{\beta }}))}^2{]}^{-1}.$

The conditions under which the two-step AE approach can be shown to be asymptotically efficient have been researched extensively. Most importantly, the kernel estimator of the score function must be (i) i.i.d., and (ii) independent of ${\mathbf{x}}_i.$ These conditions are restrictive and not easy to verify in practice; see, e.g., Yuan and De Gooijer (2007). Bickel (1982) solved the i.i.d. problem by splitting the sample in two; one sub-sample to estimate the score and another to solve for $\mathit{\beta }.$ However, Manski (1984) finds that the estimator works much better when the sample is not split, i.e., if the estimated score and $\tilde{\mathit{\beta }}$ are both computed using the entire sample. If (i) and (ii) are satisfied, a sufficient condition for adaptiveness is that (iii) $E[{(f\prime (y_i|\mathit{\beta })f^{-1}(y_i|\mathit{\beta }){\widehat{f}}_n^{\prime }(y_i|\tilde{\mathit{\beta }}){\widehat{f}}_n^{-1}(y_i|\tilde{\mathit{\beta }}))}^2]\to 0,$ as $n\to \infty .$

Since ${\widehat{f}}_n$ is present in the denominator of ${\widehat{\mathbf{d}}}_{\beta },$ unstable estimates may follow for near-zero values of ${\widehat{f}}_n(·).$ Hence, Bickel (1982) suggests to trim the estimator of the kernel score as follows $$\frac{{\widehat{f}}_n^{\prime }(y_i|\tilde{\mathit{\beta }})}{{\widehat{f}}_n(y_i|\tilde{\mathit{\beta }})}=\{\begin{array}{cc}\frac{{\widehat{f}}_n^{\prime }(y_i|\tilde{\mathit{\beta }})}{{\widehat{f}}_n(y_i|\tilde{\mathit{\beta }})},\hfill & \text{if}\mathrm{ }\left|y_i-{\mathbf{x}}_i^{\prime }\mathit{\beta }\right|\le t_1,\mathrm{ }{\widehat{f}}_n(y_i|\tilde{\mathit{\beta }})>t_2,\mathrm{ }\text{and}\mathrm{ }\frac{{\widehat{f}}_n^{\prime }(y_i|\tilde{\mathit{\beta }})}{{\widehat{f}}_n(y_i|\tilde{\mathit{\beta }})}<t_3,\hfill \\ 0,\hfill & \text{otherwise}.\hfill \end{array}$$

This trimming mechanism ensures that near-zero values do not have unreasonably large influence on the estimate. If $t_1\to \infty ,$ $t_2\to 0,t_3\to \infty ,h_n\to 0,t_1/n{h}_n^3\to 0,$ and $h_n{t}_1\to 0$ as $n\to \infty$ then condition (iii) is satisfied. Hence, adaptiveness is established under the proper trimming parameters and conditions (i) and (ii).

Naturally, the growth rates of the trimming parameters are of little use to a practitioner, and as such the choice for the trimming parameter is a practical disadvantage. Hsieh and Manski (1987) reduce the problem to selecting a one-dimensional parameter t by suggesting the following relation between the trimming parameters: $t_1=t,t_2=\hspace{0.17em}\text{exp}\hspace{0.17em}(-t^2/2),$ $t_3=t.$ These authors vary t between 3, 4, and 8. For a sample size of 50, t = 8 works best in almost all cases under study.

3.2. LGMM and LGMMS methods

Newey (1988) describes a two-step AE that avoids kernel estimation. His approach is based on moment conditions that can be derived from certain assumptions on the error distribution. Two situations are analyzed. First, the case where the error terms are i.i.d. and independent of ${\mathbf{x}}_i.$ This model implies the moment condition that any function of the errors are uncorrelated with any function of the regressors. Second, the case where the distribution of ${\epsilon }_i$ is symmetric (S) around zero conditional on ${\mathbf{x}}_i.$ The assumption that the errors are symmetrically distributed around zero yields the moment conditions that any odd function of the errors are uncorrelated with any function of the regressors. Hence, in both situations we can exploit moment restrictions. In particular, in the first case, we refer to the linearized general method of moment estimator as LGMM. In the second case, we use the short-hand notation LGMMS. For LGMM, natural moment conditions arise from the fact that $E[{\mathbf{x}}_i({\epsilon }_i^j-E({\epsilon }_i^j))]=\mathbf{0}$ for $j=1,2,\dots .$ However, Newey (1988) finds that these high-order “raw” moments, $m_j({\epsilon }_i)={\epsilon }_i^j,$ are sensitive to a fat-tailed error distribution.

Estimates that are more robust against fat tails can be obtained using the “transformed” powers with $m_j({\epsilon }_i)={({\epsilon }_i/(1+\left|{\epsilon }_i\right|))}^j$ or the “weighted” powers with $m_j({\epsilon }_i)=\hspace{0.17em}\text{exp}\hspace{0.17em}(-{\epsilon }_i/2){\epsilon }_i^j.$ Similarly, for LGMMS we may use $E[{\mathbf{x}}_i{\epsilon }_i^{2j-1}]=\mathbf{0}$ $(j=1,2,\dots ).$ As for LGMM, performance may be improved if we use the odd powers of the transformed method instead. However, for technical reasons the weighted powers can not be used for LGMMS; Newey (1988). In general, both for LGMM and LGMMS, we use the moment condition $E[{\mathbf{x}}_i(m_j({\epsilon }_i)-{\mu }_j)]=\mathbf{0}$ $(j=1,2,\dots )$ where ${\mu }_j=E[m_j({\epsilon }_i)].$

To define the LGMM(S) estimator, we introduce the following notation for some fixed value $J=J(n)$ of j, (11) $${\mathit{\zeta }}_i=(m_1({\epsilon }_i)-{\mu }_1,\dots ,m_J({\epsilon }_i)-{\mu }_J)\prime , \mathbf{w}=E[(m_{1,\epsilon }({\epsilon }_i),\dots ,m_{J,\epsilon }({\epsilon }_i))\prime ] {\mathbf{V}}_{\zeta \zeta }=\text{Cov}({\mathit{\zeta }}_i)$$(11) where $m_{j,\epsilon }(·)=\partial m_j(·)/\partial \epsilon .$ Let ${\left\{{\widehat{\epsilon }}_i\right\}}_{i=1}^n$ denote the residuals corresponding to the initial estimate $\widehat{\mathit{\beta }},$ then the quantities in (11(11) $${\mathit{\zeta }}_i=(m_1({\epsilon }_i)-{\mu }_1,\dots ,m_J({\epsilon }_i)-{\mu }_J)\prime , \mathbf{w}=E[(m_{1,\epsilon }({\epsilon }_i),\dots ,m_{J,\epsilon }({\epsilon }_i))\prime ] {\mathbf{V}}_{\zeta \zeta }=\text{Cov}({\mathit{\zeta }}_i)$$(11) ) can be consistently estimated by their corresponding sample statistics, i.e., $${\widehat{\mathit{\zeta }}}_i=(m_1({\widehat{\epsilon }}_i)-{\widehat{\mu }}_1,\dots ,m_J(\widehat{{\epsilon }_i})-{\widehat{\mu }}_J)\prime ,\hspace{1em}\widehat{\mathbf{w}}=(n^{-1}\sum _i{m}_{1,\epsilon }({\widehat{\epsilon }}_i),\dots ,n^{-1}\sum _i{m}_{J,\epsilon }({\widehat{\epsilon }}_i))\prime$$ and $${\widehat{\mathbf{V}}}_{\zeta \zeta }=n^{-1}\sum _i{\widehat{\mathit{\zeta }}}_i{\widehat{\mathit{\zeta }}}_i^{\prime }$$

The LGMM(S) estimator is given by (12) $$\begin{array}{c}{\widehat{\mathit{\beta }}}_{\text{LGMM}(\mathrm{S})}=\widehat{\mathit{\beta }}+{\left[\left(\widehat{\mathbf{w}}\prime \otimes \mathbf{X}\prime \mathbf{X}\right)\left({\widehat{\mathbf{V}}}_{\zeta \zeta }^{-1}\otimes {\left[\mathbf{X}\prime \mathbf{X}\right]}^{-1}\right)\left(\widehat{\mathbf{w}}\otimes \mathbf{X}\prime \mathbf{X}\right)\right]}^{-1}\\ \hspace{1em}\times \left(\widehat{\mathbf{w}}\prime \otimes \mathbf{X}\prime \mathbf{X}\right)\left({\widehat{\mathbf{V}}}_{\zeta \zeta }^{-1}\otimes {\left[\mathbf{X}\prime \mathbf{X}\right]}^{-1}\right)\left({\mathbf{I}}_J\otimes \mathbf{X}\prime \right)\text{vec}(\widehat{\mathit{\zeta }}),\end{array}$$(12) where $\widehat{\mathit{\zeta }}$ is the n × J matrix $({\widehat{\mathit{\zeta }}}_1^{\prime },\dots ,{\widehat{\mathit{\zeta }}}_n^{\prime })\prime ,$ ${\mathbf{I}}_J$ is an J × J identity matrix, and X an n × p matrix with its first column an $n\times 1$ vector of ones. Under certain assumptions, Newey (1988) proves asymptotic normality of the LGMM and LGMMS estimators. In particular, it should hold that $J\to \infty$ and $J\hspace{0.17em}\text{ln}\hspace{0.17em}J/\hspace{0.17em}\text{ln}\hspace{0.17em}n\to 0,$ as $n\to \infty .$ Only for LGMMS, asymptotic efficiency is obtained, but not for LGMM. By means of simulation, Newey (1988) finds for LGMM that J = 3 performs best for sample sizes between n = 50 and n = 200. However, the MSE efficiency of the estimator as a function of J flattens out as n increases. Also, the transformed method is in general preferred over the weighted method.

3.3. KDRE method

The KDRE method of Yao and Zhao (2013) can be viewed as the unconstrained two-step version of M-KDRE. That is, it follows from unconstrained maximization of the kernel likelihood function that is estimated on the basis of the residuals corresponding to an initial estimate. Under conditions (i)–(v) of Theorem 2.3 these authors prove that the KDRE method is adaptive. For technical reasons, this property is proven for a trimmed version. The untrimmed maximizer of the kernel-based likelihood is the solution to $n^{-1}{\sum }_{i=1}^n\{{\widehat{f}}_n^{\prime }(y_i|\mathit{\beta })/{\widehat{f}}_n(y_i|\mathit{\beta })\}{\mathbf{x}}_i^{\prime }=0.$ The trimmed version is then defined as the solution to $$\frac{1}{n}\sum _{i=1}^n\frac{{\widehat{f}}_n^{\prime }(y_i|\mathit{\beta })}{{\widehat{f}}_n(y_i|\mathit{\beta })}{\mathbf{x}}_i^{\prime }{G}_b\left({\widehat{f}}_n\left(y_i|\mathit{\beta }\right)\right)=0.$$

Here, (13) $$G_b(x)=\{\begin{array}{cc}0,\hfill & \text{if}\mathrm{ }x<b,\hfill \\ {\int }_b^x{g}_b(z)\mathrm{d}z\hfill & \text{if}\mathrm{ }b\le x\le 2b,\hfill \\ 1,\hfill & \text{if}\mathrm{ }x>2b,\hfill \end{array}$$(13) where $g_b(·)$ is a four times continuously differentiable function with support on $[b, 2b],$ and $b\to 0$ if $n\to \infty .$ This trimming function is introduced by Linton and Xiao (2007) and they suggest the use of the beta function. For the purpose of KDRE, the trimming parameter is only used to simplify the proof and is not a part of the actual implementation of the method.

3.4. YDG method

Yuan and De Gooijer (2007) (hereafter YDG) propose another estimator of ${\mathit{\beta }}_0$ based on estimating the error density by means of a kernel. The method is a one-step approach and as such does not require an initial estimate. The proposed estimator is given by (14) $$\widehat{\mathit{\beta }}=\arg \underset{\mathit{\beta }\in \mathcal{B}}{\max }\sum _{i=1}^n\hspace{0.17em}\text{ln}\hspace{0.17em}\frac{1}{(n-1)h_n}\sum _{j\ne i}^{n}K(\frac{r(y_i-{\mathbf{x}}_i^{\prime }\mathit{\beta })-r(y_j-{\mathbf{x}}_j^{\prime }\mathit{\beta })}{h_n})$$(14) where $r(z)=10\times e^z/(1+e^z).$ The nonlinear function $r(·)$ is introduced to avoid cancelation of the intercept coefficient in ${\mathbf{x}}_j^{\prime }\mathit{\beta }-{\mathbf{x}}_i^{\prime }\mathit{\beta }.$ However, as Yao and Zhao (2013) note, this comes with an efficiency loss; r(z) = z is efficient in the sense that even though the intercept is canceled out, the slope coefficients are adaptively estimated. They suggest the use of the following estimator (15) $${\widehat{\mathit{\beta }}}_{\text{YDG}}^{*}=\arg \underset{\mathit{\beta }\in \mathcal{B}}{\max }\sum _{i=1}^n\hspace{0.17em}\text{ln}\hspace{0.17em}\frac{1}{(n-1)h_n}\sum _{j\ne i}^{n}K(\frac{(y_i-{\mathbf{x}}_i^{*\prime }{\mathit{\beta }}^{*})-(y_j-{\mathbf{x}}_j^{*\prime }{\mathit{\beta }}^{*})}{h_n})$$(15) where ${\mathbf{x}}_i^{*}={(1, {\mathbf{x}}_i^{\prime })}^{\prime },$ and ${\widehat{\mathit{\beta }}}_{\text{YDG}}={({\widehat{\alpha }}_{\text{YDG}}, {\widehat{\mathit{\beta }}}_{\text{YDG}}^{*\prime })}^{\prime }$ and ${\widehat{\alpha }}_{\text{YDG}}=n^{-1}\sum _i(y_i-{\mathbf{x}}_i^{*\prime }{\widehat{\mathit{\beta }}}_{\text{YDG}}^{*}).$ The intercept estimate ${\widehat{\alpha }}_{\text{YDG}},$ however, is not in general asymptotically efficient.

4. Simulation study

4.1. Setup

In order to assess the finite sample practical performance of all reviewed AEs, we conduct a Monte Carlo study. We generate i.i.d. data ${\{(x_i,y_i)\}}_{i=1}^n$ from the regression model (16) $$y_i={\mathbf{x}}_i^{\prime }\mathit{\beta }+{\epsilon }_i,\hspace{1em}\mathit{\beta }={(1, -1, 2, -0.5, 3, 1, -1, 2, -0.5, 3)}^{\prime }$$(16) where $\mathit{\beta }$ is a $p\times 1$ parameter vector containing an intercept and the parameters corresponding to p − 1 explanatory variables. Here p = 10, but we also consider the case p = 2 and p = 5 consisting of the first two and first five coefficients of $\mathit{\beta },$ respectively. The sample size is set at $n=50, 100, 500,$ and $1,\mathrm{}000.$ All simulation results are based on 500 replications. The explanatory variables in ${\mathbf{x}}_i$ are independent realizations of an N(0, 1) distribution. The errors ${\epsilon }_i$ are i.i.d., and we consider the following eight error distributions:

1. standard normal;

2. variance-contaminated normal, the mixture $0.9N(0,1/9)+0.1N(0,9);$

3. t-distribution with two degrees of freedom;

4. bimodal symmetric mixture of two normals, $0.5N(-3,1)+0.5N(3,1);$

5. Unif$\left[-\sqrt{3};\sqrt{3}\right];$

6. Gamma(2,2);

7. skewed mixture of normals, $0.3N(-1.4,1)+0.7N(0.6,0.16);$ and

8. log-normal, being the distribution of $\text{exp}\hspace{0.17em}(z)$ for $z\sim N(0,1).$

The distributions are centered and scaled to have mean zero and unit variance, when necessary and possible. The t(2)-distribution is left unscaled as its variance is infinite.

If required, we use ${\widehat{\mathit{\beta }}}_{\text{LSE}}$ as an initial parameter estimate. In addition, we adopt the standard normal kernel density $K(·).$ For the SBS method, we set the trimming parameter at t = 8. Following Newey (1988), we compute the LGMM and LGMMS estimators for the transformed moments with J = 3. Implementing the kernel density-based estimators requires a method for choosing the bandwidth h_n. There is a vast literature on this topic, ranging from simple to involved methods. But none of the proposed methods has overall performance. In an extensive simulation study of model (1(1) $$y_i={\mathbf{x}}_i^{\prime }{\mathit{\beta }}_0+{\epsilon }_i,\hspace{1em}(i=1,\dots ,n)$$(1) ) with n = 100 and p = 2, Reichardt (2017) concludes that for M-KDRE $h_{n,1}=1.06{\widehat{\sigma }}_n{n}^{-1/5}$ is preferable for symmetric error distributions in terms of root mean squared error (RMSE) of the estimators. Here ${\widehat{\sigma }}_n$ is the standard deviation of the data. For skewed distributions, he recommends $h_{n,2}=0.9A{n}^{-1/5}$ where $A=\min ({\widehat{\sigma }}_n,R/1.34)$ with R the inter-quartile range of the data. The KDRE and YDG estimators perform best under $h_{n,1}.$ The SBS estimator generally shows the smallest RMSE for $h_{n,2}.$ Hence, throughout the simulations, we use the above bandwidths.

4.2. Results

Averaged over all replications, Table 1 provides summary information on the computed RMSEs of the slope and intercept coefficients for all n and p, and across all error distributions. Clearly, M-KDRE shows the best overall performance of all estimators for the slope coefficient, with 64 lowest RMSE values out of a total of 84, i.e., 7 estimation methods (M-KDRE, KDRE, YDG, SBS, LGMM, LGMMS, LSE), 3 values for p, and 4 values for n. On the other hand, the SBS estimator has only one lowest RMSE value for $n=1,\mathrm{}000$ and p = 2. The other estimators have low RMSE values lying in between the above two values, with the LSE results as a benchmark. Finally, from the last column of Table 1, we see that M-KDRE performs equally well with YDG, LGGM and LGMMS in estimating the intercept term, and M-KDRE markedly outperforms KDRE.

Table 1. Number of times the RMSE attains its lowest value for seven regression estimators, for each n and p and across all eight error distribution functions, for the slope coefficient and, in parentheses, for the intercept.

p	2				5				10
n	50	100	500	1,000	50	100	500	1,000	50	100	500	1,000	Total
M-KDRE	5${}_{(3)}$	6${}_{(3)}$	5${}_{(3)}$	6${}_{(2)}$	$4_{(1)}$	5${}_{(2)}$	6${}_{(2)}$	6${}_{(4)}$	4${}_{(3)}$	4${}_{(3)}$	6${}_{(3)}$	7${}_{(3)}$	64${}_{(32)}$
KDRE	0${}_{(0)}$	0${}_{(0)}$	3${}_{(1)}$	3${}_{(0)}$	0${}_{(0)}$	0${}_{(0)}$	2${}_{(1)}$	3${}_{(1)}$	0${}_{(1)}$	0${}_{(0)}$	1${}_{(0)}$	4${}_{(1)}$	16${}_{(5)}$
YDG	1${}_{(2)}$	0${}_{(3)}$	2${}_{(5)}$	2${}_{(4)}$	3${}_{(1)}$	1${}_{(1)}$	3${}_{(4)}$	3${}_{(5)}$	1${}_{(0)}$	3${}_{(1)}$	2${}_{(5)}$	3${}_{(5)}$	24${}_{(36)}$
SBS	0${}_{(0)}$	0${}_{(0)}$	0${}_{(0)}$	1${}_{(0)}$	0${}_{(0)}$	0${}_{(0)}$	0${}_{(1)}$	0${}_{(1)}$	0${}_{(0)}$	0${}_{(0)}$	0${}_{(0)}$	0${}_{(1)}$	1${}_{(3)}$
LGMM	1${}_{(1)}$	2${}_{(4)}$	4${}_{(3)}$	3${}_{(3)}$	1${}_{(5)}$	0${}_{(4)}$	2${}_{(1)}$	4${}_{(4)}$	1${}_{(1)}$	1${}_{(2)}$	1${}_{(3)}$	4${}_{(3)}$	24${}_{(34)}$
LGMMS	0${}_{(4)}$	1${}_{(4)}$	3${}_{(4)}$	2${}_{(3)}$	1${}_{(3)}$	1${}_{(3)}$	0${}_{(4)}$	3${}_{(4)}$	1${}_{(3)}$	1${}_{(3)}$	1${}_{(3)}$	3${}_{(4)}$	17${}_{(42)}$
LSE	1${}_{(1)}$	1${}_{(3)}$	1${}_{(3)}$	1${}_{(3)}$	1${}_{(1)}$	1${}_{(1)}$	1${}_{(2)}$	1${}_{(3)}$	1${}_{(1)}$	1${}_{(1)}$	1${}_{(2)}$	1${}_{(3)}$	12${}_{(24)}$
Total	8	10	18	18	10	8	14	20	8	10	12	22	158
	(11(11) $${\mathit{\zeta }}_i=(m_1({\epsilon }_i)-{\mu }_1,\dots ,m_J({\epsilon }_i)-{\mu }_J)\prime , \mathbf{w}=E[(m_{1,\epsilon }({\epsilon }_i),\dots ,m_{J,\epsilon }({\epsilon }_i))\prime ] {\mathbf{V}}_{\zeta \zeta }=\text{Cov}({\mathit{\zeta }}_i)$$(11) )	(17(17) $$Y_{i,j}={\alpha }_j+{\beta }_j{\text{RC}}_i+{\gamma }_j{\text{Y}}_{i,j}^{*}+{\delta }_j{\text{X}}_{i,j}+{\epsilon }_i,\hspace{1em}(i=1,\dots ,n;j=1, 2, 3)$$(17) )	(19)	(15(15) $${\widehat{\mathit{\beta }}}_{\text{YDG}}^{*}=\arg \underset{\mathit{\beta }\in \mathcal{B}}{\max }\sum _{i=1}^n\hspace{0.17em}\text{ln}\hspace{0.17em}\frac{1}{(n-1)h_n}\sum _{j\ne i}^{n}K(\frac{(y_i-{\mathbf{x}}_i^{*\prime }{\mathit{\beta }}^{*})-(y_j-{\mathbf{x}}_j^{*\prime }{\mathit{\beta }}^{*})}{h_n})$$(15) )	(11(11) $${\mathit{\zeta }}_i=(m_1({\epsilon }_i)-{\mu }_1,\dots ,m_J({\epsilon }_i)-{\mu }_J)\prime , \mathbf{w}=E[(m_{1,\epsilon }({\epsilon }_i),\dots ,m_{J,\epsilon }({\epsilon }_i))\prime ] {\mathbf{V}}_{\zeta \zeta }=\text{Cov}({\mathit{\zeta }}_i)$$(11) )	(11(11) $${\mathit{\zeta }}_i=(m_1({\epsilon }_i)-{\mu }_1,\dots ,m_J({\epsilon }_i)-{\mu }_J)\prime , \mathbf{w}=E[(m_{1,\epsilon }({\epsilon }_i),\dots ,m_{J,\epsilon }({\epsilon }_i))\prime ] {\mathbf{V}}_{\zeta \zeta }=\text{Cov}({\mathit{\zeta }}_i)$$(11) )	(15(15) $${\widehat{\mathit{\beta }}}_{\text{YDG}}^{*}=\arg \underset{\mathit{\beta }\in \mathcal{B}}{\max }\sum _{i=1}^n\hspace{0.17em}\text{ln}\hspace{0.17em}\frac{1}{(n-1)h_n}\sum _{j\ne i}^{n}K(\frac{(y_i-{\mathbf{x}}_i^{*\prime }{\mathit{\beta }}^{*})-(y_j-{\mathbf{x}}_j^{*\prime }{\mathit{\beta }}^{*})}{h_n})$$(15) )	(22)	(9(9) $$p_{ij,(k+1)}^{(u)}=\frac{\hspace{0.17em}\text{exp}\hspace{0.17em}\{-\frac{1}{2h_n^2}{(y_i-{\mathbf{x}}_i^{\prime }{\widehat{\mathit{\beta }}}_{(k)}^{(u)}-{\widehat{\epsilon }}_j^{(u-1)})}^2\}}{\sum _{\mathcal{\ell }=1}^n\hspace{0.17em}\text{exp}\hspace{0.17em}\{-\frac{1}{2h_n^2}{(y_i-{\mathbf{x}}_i^{\prime }{\widehat{\mathit{\beta }}}_{(k)}^{(u)}-{\widehat{\epsilon }}_{\mathcal{\ell }}^{(u-1)})}^2\}}$$(9) )	(10(10) $$\begin{array}{c}{\widehat{\mathit{\beta }}}_{(k+1)}^{(u)}={\widehat{\mathit{\beta }}}_{\text{LSE}}-{\left(\sum _{i=1}^n{\mathbf{x}}_i{\mathbf{x}}_i^{\prime }\right)}^{-1}\sum _{i=1}^n\left({\mathbf{x}}_i\sum _{j=1}^n{p}_{ij,(k+1)}^{(u)}{\widehat{\epsilon }}_j^{(u-1)}\right)\\ \hspace{1em}+\frac{1}{n}{\left(\sum _{i=1}^n{\mathbf{x}}_i{\mathbf{x}}_i^{\prime }\right)}^{-1}\sum _{i=1}^n{\mathbf{x}}_i\left(\sum _{i=1}^n\sum _{j=1}^n{p}_{ij,(k+1)}^{(u)}{\widehat{\epsilon }}_j^{(u-1)}\right)\end{array}$$(10) )	(16(16) $$y_i={\mathbf{x}}_i^{\prime }\mathit{\beta }+{\epsilon }_i,\hspace{1em}\mathit{\beta }={(1, -1, 2, -0.5, 3, 1, -1, 2, -0.5, 3)}^{\prime }$$(16) )	(20)	(176)

Note: For the purpose of this table, the RMSE is measured up to three decimal figures such that more than one estimator may have the lowest RMSE.

Reichardt (2017) reports RMSE values for each of the eight error distributions. For the sake of space we omit details. However, for the slope coefficient the simulation results can be summarized as follows.

• In terms of RMSE, the M-KDRE method performs very well for the log-normal error distribution (h). That is, the RMSE of the second most efficient estimators (KDRE and LGMM) is approximately 40% larger, even for $n=1,\mathrm{}000.$ Under error distributions (b) and (c), M-KDRE is also most efficient, but here the efficiency is gained mostly for n = 50 and n = 100. Furthermore, M-KDRE has a superior performance in small samples of the t(2) error distribution. It is also close to best for error distributions (d)–(g).

• The YDG method performs well for error distributions (d)–(g), but fails quite dramatically for error distributions with fat tails, i.e. (b), (c), and (h).

• Efficiency of the SBS estimator is in general low relative to alternatives, but performance is especially weak under error distributions (c) and (d).

• Overall, LGMM is a reasonable estimator, but its efficiency is lost under error distributions (e) and (f). This efficiency loss persists even for $n=1,\mathrm{}000.$

• LGMMS is by construction inefficient when the error distribution is skewed: (f)–(h). More interestingly, the LGMMS-estimate of the slope coefficient is no improvement over LGMM under symmetric error distributions. The inefficiency of LGMMS with respect to LGMM is likely to be due to the fact that LGMMS uses moment restrictions on odd powers of the disturbances only and, hence, for a particular value of J, uses higher order moments that may lead to less efficient estimation.

Table 2 shows summary results for the bias for both slope and intercept estimators for all n and p, and across all error distributions. For the slope coefficient, M-KDRE has the best performance in terms of the lowest bias values. Again, from Reichardt (2017) we learn that the intercept bias of the different estimators is usually of similar magnitude in the symmetric cases. Under the asymmetric error distributions, the bias of the intercept is much larger for KDRE, SBS and LGMMS than for the other estimators.

Table 2. Number of times the bias of seven regression estimators attains it lowest value for each n and p, and across all eight error distribution functions for the slope coefficient and, in parentheses, for the intercept.

p	2				5				10
n	50	100	500	1,000	50	100	500	1,000	50	100	500	1,000	Total
M-KDRE	0${}_{(4)}$	3${}_{(2)}$	7${}_{(5)}$	8${}_{(6)}$	4${}_{(1)}$	6${}_{(6)}$	5${}_{(6)}$	8${}_{(3)}$	5${}_{(1)}$	7${}_{(3)}$	7${}_{(5)}$	6${}_{(5)}$	66${}_{(47)}$
KDRE	0${}_{(1)}$	1${}_{(2)}$	7${}_{(3)}$	7${}_{(3)}$	3${}_{(0)}$	4${}_{(3)}$	6${}_{(3)}$	8${}_{(3)}$	0${}_{(0)}$	4${}_{(1)}$	6${}_{(2)}$	7${}_{(2)}$	53${}_{(23)}$
YDG	3${}_{(2)}$	3${}_{(4)}$	4${}_{(3)}$	7${}_{(6)}$	0${}_{(2)}$	2${}_{(4)}$	5${}_{(5)}$	5${}_{(4)}$	2${}_{(5)}$	2${}_{(1)}$	6${}_{(8)}$	5${}_{(6)}$	44${}_{(50)}$
SBS	3${}_{(2)}$	1${}_{(1)}$	5${}_{(3)}$	4${}_{(3)}$	3${}_{(0)}$	1${}_{(2)}$	7${}_{(5)}$	7${}_{(3)}$	3${}_{(1)}$	2${}_{(0)}$	6${}_{(3)}$	4${}_{(1)}$	46${}_{(24)}$
LGMM	2${}_{(4)}$	3${}_{(2)}$	5${}_{(4)}$	5${}_{(6)}$	1${}_{(1)}$	4${}_{(2)}$	6${}_{(5)}$	8${}_{(3)}$	3${}_{(1)}$	3${}_{(2)}$	5${}_{(6)}$	4${}_{(4)}$	49${}_{(40)}$
LGMMS	3${}_{(3)}$	3${}_{(3)}$	5${}_{(4)}$	6${}_{(5)}$	2${}_{(4)}$	3${}_{(1)}$	5${}_{(4)}$	6${}_{(4)}$	1${}_{(2)}$	2${}_{(4)}$	6${}_{(1)}$	5${}_{(4)}$	47${}_{(39)}$
LSE	2${}_{(5)}$	1${}_{(2)}$	2${}_{(5)}$	4${}_{(6)}$	0${}_{(3)}$	0${}_{(5)}$	2${}_{(4)}$	2${}_{(3)}$	1${}_{(1)}$	1${}_{(3)}$	2${}_{(5)}$	1${}_{(2)}$	18${}_{(45)}$
Total	13	15	35	41	13	20	36	44	15	21	38	32	323
	(21)	(16(16) $$y_i={\mathbf{x}}_i^{\prime }\mathit{\beta }+{\epsilon }_i,\hspace{1em}\mathit{\beta }={(1, -1, 2, -0.5, 3, 1, -1, 2, -0.5, 3)}^{\prime }$$(16) )	(27)	(35)	(11(11) $${\mathit{\zeta }}_i=(m_1({\epsilon }_i)-{\mu }_1,\dots ,m_J({\epsilon }_i)-{\mu }_J)\prime , \mathbf{w}=E[(m_{1,\epsilon }({\epsilon }_i),\dots ,m_{J,\epsilon }({\epsilon }_i))\prime ] {\mathbf{V}}_{\zeta \zeta }=\text{Cov}({\mathit{\zeta }}_i)$$(11) )	(23)	(32)	(23)	(11(11) $${\mathit{\zeta }}_i=(m_1({\epsilon }_i)-{\mu }_1,\dots ,m_J({\epsilon }_i)-{\mu }_J)\prime , \mathbf{w}=E[(m_{1,\epsilon }({\epsilon }_i),\dots ,m_{J,\epsilon }({\epsilon }_i))\prime ] {\mathbf{V}}_{\zeta \zeta }=\text{Cov}({\mathit{\zeta }}_i)$$(11) )	(14(14) $$\widehat{\mathit{\beta }}=\arg \underset{\mathit{\beta }\in \mathcal{B}}{\max }\sum _{i=1}^n\hspace{0.17em}\text{ln}\hspace{0.17em}\frac{1}{(n-1)h_n}\sum _{j\ne i}^{n}K(\frac{r(y_i-{\mathbf{x}}_i^{\prime }\mathit{\beta })-r(y_j-{\mathbf{x}}_j^{\prime }\mathit{\beta })}{h_n})$$(14) )	(30)	(25)	(268)

5. Empirical application

Andrabi, Das, and Khwaja (2017) study the impact of providing information in the form of school report cards on educational outcomes such as school fees, test scores, and enrollment in markets with multiple public and private providers. The report cards given to both households and schools in n randomly sampled villages across three districts, in the Punjab province of Pakistan, include information on the performance of the child, the average score of different schools in the village, and the average village score in mathematics, English, and Urdu. The following three linear regression models are of interest: (17) $$Y_{i,j}={\alpha }_j+{\beta }_j{\text{RC}}_i+{\gamma }_j{\text{Y}}_{i,j}^{*}+{\delta }_j{\text{X}}_{i,j}+{\epsilon }_i,\hspace{1em}(i=1,\dots ,n;j=1, 2, 3)$$(17) where $Y_{i,1},Y_{i,2},$ and $Y_{i,3}$ are average fees, test scores, and enrollment rate in the post-intervention year of village i, respectively. $Y_{i,j}^{*}$ denotes the baseline measurement of the same variables. ${\text{RC}}_i$ is the treatment dummy assignment to village i, which makes β_j the variable of interest, an estimate of the impact of the report card assignment. $X_{i,j}$ is a vector of village-level baseline controls. All models in the paper are estimated using LSE.

Table 3, column 1, shows the LSEs of β_j $(j=1, 2, 3)$ and their corresponding standard errors (in parentheses) as shown in, respectively, Table 3(1) panel C, Table 3(4) panel C, and Table 4(1) panel C of Andrabi, Das, and Khwaja (2017). The Shapiro-Wilk test for normal data indicates that the LSE residuals are far from normally distributed, with p-values 0.000 $(j=1,n=104),$ 0.002 $(j=2,n=112),$ and 0.000 $(j=3,n=112).$ Indeed, in all cases, diagnostic statistics show that the residuals have fatter tails than one would expect based on normality. Based on the LSE results, Andrabi, Das, and Khwaja (2017) report the following main findings. First, private schools decreased their annualized fees $(Y_{i,1})$ by an average of 187 rupees, about 17% of their baseline fees, in response to the report card intervention. Second, test scores $(Y_{i,2})$ increased by 0.11 standard deviation. Third, primary enrollment $(Y_{i,3})$ increased by 3.2 percent points or 4.5% in treatment villages.

Table 3. Effect of report cards on school fees, test scores, and enrollment as given by parameter estimates of β_j $(j=1, 2, 3)$ using seven estimation methods. For columns 2–7, standard errors (in parentheses) are based on 500 bootstrap replicates.

	LSE	M-KDRE	KDRE	YDG	SBS	LGMM	LGMMS
${\widehat{\beta }}_1$	−187.0	−4.418	−85.55	−154.4	−47.74	−172.3	−182.9
	(65.91)	(38.76)	(63.42)	(59.90)	(73.34)	(64.54)	(67.65)
${\widehat{\beta }}_2$	0.114	0.084	0.092	0.094	0.088	0.088	0.099
	(0.046)	(0.050)	(0.048)	(0.060)	(0.054)	(0.045)	(0.050)
${\widehat{\beta }}_3$	0.032	0.028	0.030	0.018	0.030	0.029	0.027
	(0.014)	(0.014)	(0.012)	(0.019)	(0.013)	(0.012)	(0.012)

Table 4. Median absolute prediction error (MAPE) of six AEs relative to LSE.

Model (j)	M-KDRE	KDRE	YDG	SBS	LGMM	LGMMS
1	0.720	0.858	0.906	0.769	0.992	1.001
2	1.002	1.016	0.998	1.005	0.997	1.020
3	0.972	0.975	1.029	0.963	0.965	0.971

Table 3, columns 2–7, shows estimates of the six AEs for models $j=1, 2,$ and 3. We see that these estimators pull the estimated treatment effect toward zero for all models. The results for model 1 are especially striking. The M-KDRE of β is more than 40 times smaller than the LSE, in absolute value. Also, for model 1, the AEs differ substantially. In that respect, it is interesting to investigate the prediction performance of the respective methods.

Table 4 shows the ratio of the median absolute prediction error (MAPE) of an estimator relative to the LSE. The training set is a random sample of the data of size $\lceil 0.8n\rceil .$ We see that M-RKDRE has the lowest MAPE for model 1. In addition, observe that the prediction performance is generally better for AEs with a low estimate of β₁ such as KDRE and SBS. This suggests that the effect of the report cards on school fees, if it exists at all, is much lower than reported. For models 2 and 3, there is less difference between the estimates of the AEs. Also, the estimates are adjusted less strongly with respect to LSE.

Table 5 reports two bootstrapped 95% confidence intervals for β_j $(j=1, 2, 3)$ as estimated by M-KDRE. The confidence interval termed “Normal” is based on an asymptotic normality assumption, and the column called “Percentile” is based on the 95% inter quartile range obtained from the empirical distribution of 500 bootstrap replicates. Both intervals show that the estimated treatment effect for model 1 is not significantly different from zero.

Table 5. Bootstrapped 95% confidence intervals of the effect of report cards for the M-KDRE method.

Model (j)	Normal	Percentile
1	[−80.38, 71.55]	[−88.11, 59.25]
2	[−0.015, 0.183]	[−0.011, 0.187]
3	[0.001, 0.055]	[0.002, 0.055]

In summary, the above results demonstrate the practical relevance of the AEs in general and that of the proposed M-KDRE method in particular. None of the other AEs adjusted the treatment effect on school fees as far toward zero as M-KDRE, while prediction performance suggests that this method should be preferred over other methods for this particular linear regression model and sample size. Thus, there is no support for the first finding of Andrabi, Das, and Khwaja (2017) at any reasonable significance level. Further, Table 3 shows that the AEs find that the effect of report cards on test scores is not significantly different from zero at the 5% nominal level. Lastly, the effect of report cards on the enrollment rate seems, even though marginally significant for M-KDRE, also questionable.

6. Summary and concluding remarks

In this paper, we proposed an adaptive multi-step kernel density-based regression estimator for linear regression models. We have established the theoretical properties of our estimation method, including asymptotic normality and almost sure convergence. In an extensive simulation study, we have shown that the finite sample performance of M-KDRE is second to none of five alternative AEs. For several error distributions, it is up to twice as efficient in terms of RMSE than the second best estimator. Further, for any other error distribution, it is either the most efficient or very close to the most efficient estimator. All other AEs show a loss of efficiency for certain specific error distributions. Our empirical application provides a good illustration of many of these issues. In particular, using the M-KDRE method and its corresponding bootstrap standard errors, we found fairly compelling evidence that the treatment effects that Andrabi, Das, and Khwaja (2017) find are not significantly different from zero.

The results raise several questions for further research. For instance, one may wish to estimate nonlinear regressions via the M-KDRE method. In that case the EM algorithm, at least in its present form, needs to be adjusted. Another issue concerns the fact that the multi-step method makes use of ${\widehat{\beta }}_{\text{LSE}}$ in the initial step. This choice was primarily based on computational convenience. Perhaps, efficiency may be further enhanced by a more prudent choice of the initial estimator. It may also be of interest to assess the robustness of M-KDRE to a violation of the independence assumption. In particular, adaptive estimation is not in general possible when the vector of covariates and the error process are not mutually independent. We leave these questions for future research.

Acknowledgments

The authors would like to thank two anonymous referees for their valuable comments and suggestions.

Appendix A: Proofs of results

Lemma A.1.

(Zhang 1990, Theorem 5) If model (1) holds, ${\left\{{\epsilon }_i\right\}}_{i=1}^n$ are i.i.d. with unknown density f(x) where $f(·)$ is a uniformly continuous function that satisfies (i) $\int xf(x)\mathrm{d}x=0$, (ii) $0<\int x^{2}f(x)\mathrm{d}x<\infty$, the set of covariates ${\left\{{\mathbf{x}}_i\right\}}_{i=1}^{\infty }$ satisfies (iii) $\exists \mathrm{ }0<M<\infty$ such that $\Vert {\mathbf{x}}_i\Vert <M$ $\forall i=1,\dots ,n$, (iv) ${\mathbf{S}}_n\to \mathbf{Q}=E(\mathbf{x}{\mathbf{x}}^{\prime })$ where ${\mathbf{S}}_n=n^{-1}{\sum }_{i=1}^n{\mathbf{x}}_i{\mathbf{x}}_i^{\prime }$, and the following assumptions on the kernel function $K(·)$ hold: (v) K(x) is uniformly bounded and $\exists \mathrm{ }0<\rho <\infty$ such that K(x) = 0 $\forall x:\Vert \mathbf{x}\Vert \ge \rho$ (vi) K(x) is Riemann integrable on $[-\rho , \rho ]$ (vii) when $n\to \infty ,0<h_n\to 0$ and $\sqrt{n}{h}_n/\hspace{0.17em}\text{ln}\hspace{0.17em}n\to \infty$, then (A.1) $$\underset{x\in \mathbb{R}}{\sup } \Vert {\widehat{f}}_{n,\text{LSE}}(x)-f(x)\Vert \stackrel{a.s.}{\to }0$$(A.1) where ${\widehat{f}}_{n,\text{LSE}}$ is the kernel density-based estimator of the $\text{LSE}$ residuals ${\widehat{\epsilon }}_{i,\text{LSE}}=y_i-{\mathbf{x}}_i^{\prime }{\widehat{\mathit{\beta }}}_{\text{LSE}}$ $(i=1,\dots ,n).$

Lemma A.2.

Suppose that the assumptions of Lemma A.1 hold. Then, if any estimator ${\mathit{\beta }}^{*}$ satisfies $Pr(\underset{n\to \infty }{\lim }\max |{\mathit{\beta }}^{*}-{\mathit{\beta }}_0|\le \max |{\mathbf{S}}_n^{-1}| \Vert \hspace{0.17em}\text{ln}\hspace{0.17em}\max |{\mathbf{S}}_n^{-1}| \Vert )=1$ where $\max \left|\mathbf{A}\right|={\max }_{i,j}\left|a_{ij}\right|$ with a_ij the elements of a matrix A, (A.2) $$\underset{x\in \mathbb{R}}{\sup }\Vert {\widehat{f}}_n^{*}(x)-f(x)\Vert \stackrel{a.s.}{\to }0$$(A.2) where ${\widehat{f}}_n^{*}$ is the kernel density-based estimator of the residuals ${\widehat{\epsilon }}_i=y_i-{\mathbf{x}}_i^{\prime }{\mathit{\beta }}^{*}$ $(i=1,\dots ,n).$

Proof.

This result follows immediately from Theorem 5 and Lemma 4 in Zhang (1990) in conjunction with Theorem 4 and (29) in Chai, Li, and Tian (1991). □

Lemma A.3.

If there is a function $Q_0(\mathit{\beta })$ such that (i) $Q_0(\mathit{\beta })$ is uniquely maximized at ${\mathit{\beta }}_0$, (ii) $\mathcal{B}$ is compact, (iii) $Q_0(\mathit{\beta })$ is continuous, (iv) $\underset{\mathit{\beta }\in \mathcal{B}}{\sup }\Vert {\widehat{Q}}_n(\mathit{\beta })-Q_0(\mathit{\beta })\Vert \stackrel{a.s.}{\to }0$, then for $u=1,2,\dots ,$ ${\widehat{\mathit{\beta }}}^{(u)}\stackrel{a.s.}{\to }{\mathit{\beta }}_0$, where $\widehat{\mathit{\beta }}$ maximizes the objective function ${\widehat{Q}}_n(\mathit{\beta })$ subject to $\mathit{\beta }\in \mathcal{B}$. The weak convergence result, i.e., $\widehat{\mathit{\beta }}\stackrel{p}{\to }{\mathit{\beta }}_0$ can be obtained by replacing condition (iv) by $\underset{\mathit{\beta }\in \mathcal{B}}{\sup }\Vert {\widehat{Q}}_n(\mathit{\beta })-Q_0(\mathit{\beta })\Vert \stackrel{p}{\to }0.$

Proof:

The proof is similar to the proof of Theorem 2.1 of Newey and McFadden (1994). □

Lemma A.4.

If $f_n:\mathcal{B}\to \mathbb{R}$ is a continuous function, $\mathcal{B}$ is compact, and $f_n\stackrel{a.s.}{\to }f$, then (A.3) $$\underset{n\to \infty }{\lim }{\int }_{\mathcal{B}}{f}_n\mathrm{d}u={\int }_{\mathcal{B}}f\mathrm{d}u$$(A.3)

Proof.

Since $\mathcal{B}$ is compact and f_n is continuous, the image $f_n(\mathcal{B})$ is a compact subset of $\mathbb{R}$ and hence, closed and bounded. Then, the result follows from the bounded convergence theorem; see, e.g., Wade (1974). □

Proof of Theorem 2.1.

Following (Newey and McFadden 1994, Thm. 2.5), we verify the conditions in Lemma A.3. Note that conditions (i)–(iii) are on the density $f(·)$ of $\mathit{\epsilon },$ and on the parameter space $\mathcal{B}.$ These conditions hold under the usual regularity conditions of MLE. Condition (iv) of Lemma A.3 implies that we have to prove that $$\underset{\mathit{\beta }\in \mathcal{B}}{\sup }\Vert \frac{1}{n}\sum _{i=1}^n\hspace{0.17em}\text{ln}\hspace{0.17em}{\widehat{f}}_n^{(1)}(y_i|\mathit{\beta })-E[\text{ln}\hspace{0.17em}f(y_i|\mathit{\beta })]\Vert \stackrel{a.s.}{\to }0.$$

Since ${\widehat{f}}_n^{(1)}={\widehat{f}}_{n,\text{LSE}},$ we have by Lemma A.1 $$\underset{x\in \mathbb{R}}{\sup }\Vert {\widehat{f}}_n^{(1)}(x)-f(x)\Vert \stackrel{a.s.}{\to }0.$$

Now note that $$\underset{\mathit{\beta }\in \mathcal{B}}{\sup }\Vert {\widehat{f}}_n^{(1)}(y_i|\mathit{\beta })-f(y_i|\mathit{\beta })\Vert \le \underset{\mathit{\beta }\in {\mathbb{R}}^p}{\sup }\Vert {\widehat{f}}_n^{(1)}(y_i|\mathit{\beta })-f(y_i|\mathit{\beta })\Vert \le \underset{x\in \mathbb{R}}{\sup } \Vert {\widehat{f}}_n^{(1)}(x)-f(x)\Vert$$ implying that (A.4) $$\underset{\mathit{\beta }\in \mathcal{B}}{\sup } \Vert {\widehat{f}}_n^{(1)}(y_i|\mathit{\beta })-f(y_i|\mathit{\beta })\Vert \stackrel{a.s.}{\to }0.$$(A.4)

Condition (iii) implies that $\underset{\mathit{\beta }\in \mathcal{B}}{\inf }f(y_i|\mathit{\beta })>0.$ Thus, $\exists \epsilon >0$ such that $\underset{\mathit{\beta }\in \mathcal{B}}{\inf }f(y_i|\mathit{\beta })>\epsilon .$ Also, by (A.4) for any $\epsilon >0,$ $$\text{Pr}(\underset{n\to \infty }{\lim }\underset{\mathit{\beta }\in \mathcal{B}}{\sup }\Vert {\widehat{f}}_n^{(1)}(y_i|\mathit{\beta })-f(y_i|\mathit{\beta })\Vert <\epsilon )=1\Rightarrow \text{Pr}(\underset{n\to \infty }{\lim }\underset{\mathit{\beta }\in \mathcal{B}}{\inf }{\widehat{f}}_n^{(1)}(y_i|\mathit{\beta })>0)=1.$$

This, together with condition (ii), ensures that for n large enough both $\text{ln}\hspace{0.17em}f(y_i|\mathit{\beta })$ and $\text{ln}\hspace{0.17em}{\widehat{f}}_n^{(1)}(y_i|\mathit{\beta })$ are uniformly continuous with probability one such that by the uniform continuous mapping theorem $$\underset{\mathit{\beta }\in \mathcal{B}}{\sup }\Vert \hspace{0.17em}\text{ln}\hspace{0.17em}{\widehat{f}}_n^{(1)}(y_i|\mathit{\beta })-\hspace{0.17em}\text{ln}\hspace{0.17em}f(y_i|\mathit{\beta })\Vert \stackrel{a.s.}{\to }0.$$

Note that by conditions (ii) and (iii), $\text{ln}\hspace{0.17em}{\widehat{f}}_n^{(1)}(y_i|\mathit{\beta })$ is bounded and we may invoke the uniform law of large numbers such that, (A.5) $$\underset{\mathit{\beta }\in \mathcal{B}}{\sup }\Vert \frac{1}{n}\sum _{i=1}^n\hspace{0.17em}\text{ln}\hspace{0.17em}{\widehat{f}}_n^{(1)}(y_i|\mathit{\beta })-E[\text{ln}\hspace{0.17em}{\widehat{f}}_n^{(1)}(y_i|\mathit{\beta })]\Vert \stackrel{a.s.}{\to }0.$$(A.5)

Also, by Lemma A.4, (A.6) $$\underset{n\to \infty }{\lim }E[\text{ln}\hspace{0.17em}{\widehat{f}}_n^{(1)}(y_i|\mathit{\beta })]=E[\text{ln}\hspace{0.17em}f(y_i|\mathit{\beta })]$$(A.6)

Now define the following variables $$\begin{array}{c}A\equiv \Vert \frac{1}{n}\sum _{i=1}^n\hspace{0.17em}\text{ln}\hspace{0.17em}{\widehat{f}}_n^{(1)}(y_i|\mathit{\beta })-E[\text{ln}\hspace{0.17em}f(y_i|\mathit{\beta })]\Vert ,\hspace{1em}{A}_1\equiv \Vert \frac{1}{n}\sum _{i=1}^n\hspace{0.17em}\text{ln}\hspace{0.17em}{\widehat{f}}_n^{(1)}(y_i|\mathit{\beta })-E[\text{ln}\hspace{0.17em}{\widehat{f}}_n^{(1)}(y_i|\mathit{\beta })]\Vert ,\\ A_2\equiv \Vert E[\text{ln}\hspace{0.17em}{\widehat{f}}_n^{(1)}(y_i|\mathit{\beta })]-E[\text{ln}\hspace{0.17em}f(y_i|\mathit{\beta })]\Vert .\end{array}$$

Then, by the triangle inequality, we have $A\le A_1+A_2,$ and by (A.5(A.5) $$\underset{\mathit{\beta }\in \mathcal{B}}{\sup }\Vert \frac{1}{n}\sum _{i=1}^n\hspace{0.17em}\text{ln}\hspace{0.17em}{\widehat{f}}_n^{(1)}(y_i|\mathit{\beta })-E[\text{ln}\hspace{0.17em}{\widehat{f}}_n^{(1)}(y_i|\mathit{\beta })]\Vert \stackrel{a.s.}{\to }0.$$(A.5) ) and (A.6(A.6) $$\underset{n\to \infty }{\lim }E[\text{ln}\hspace{0.17em}{\widehat{f}}_n^{(1)}(y_i|\mathit{\beta })]=E[\text{ln}\hspace{0.17em}f(y_i|\mathit{\beta })]$$(A.6) ), $\underset{\mathit{\beta }\in \mathcal{B}}{\sup }{A}_1\stackrel{a.s.}{\to }0$ and $\underset{n\to \infty }{\lim }\underset{\mathit{\beta }\in \mathcal{B}}{\sup }{A}_2=0.$ From condition (iv) of Lemma A.3 it follows that (A.7) $$\underset{\mathit{\beta }\in \mathcal{B}}{\sup }\Vert \frac{1}{n}\sum _{i=1}^n\hspace{0.17em}\text{ln}\hspace{0.17em}{\widehat{f}}_n^{(1)}(y_i|\mathit{\beta })-E[\text{ln}\hspace{0.17em}f(y_i|\mathit{\beta })]\Vert \stackrel{a.s.}{\to }0.$$(A.7)

Thus, by Lemma A.3, (A.8) $${\widehat{\mathit{\beta }}}^{(1)}\stackrel{a.s.}{\to }{\mathit{\beta }}_0$$(A.8)

For the sake of completeness, we prove also that the constraint $c(\mathit{\beta })=0$ does not affect this result. Let $\mathcal{C}\subseteq \mathcal{B}$ be the subset for which $c(\mathit{\beta })=0.$ That is, $\mathcal{C}=\{\mathit{\beta }\in \mathcal{B}:c(\mathit{\beta })=0\}.$ First, note that $\mathcal{C}$ is the level set of the continuous function $c(\mathit{\beta })$ such that $\mathcal{C}$ is closed. Also, $\mathcal{C}$ is bounded since $\mathcal{C}\subseteq \mathcal{B}$ and $\mathcal{B}$ is bounded. Hence, $\mathcal{C}$ is compact such that ${\widehat{\mathit{\beta }}}^{(1)}=\arg \underset{\mathit{\beta }\in \mathcal{C}}{\sup }{\sum }_{i=1}^n\hspace{0.17em}\text{ln}\hspace{0.17em}{\widehat{f}}_n^{(1)}(y_i|\mathit{\beta }).$ Denote $\tilde{\mathit{\beta }}=\arg \underset{\mathit{\beta }\in \mathcal{B}}{\sup }{\sum }_{i=1}^n\hspace{0.17em}\text{ln}\hspace{0.17em}{\widehat{f}}_n^{(1)}(y_i|\mathit{\beta })$ as the global maximizer of the objective function over $\mathcal{B}.$ (Newey and McFadden 1994, p. 2122) show that for (A.8) to hold, it suffices to prove that (A.9) $$\frac{1}{n}\sum _{i=1}^n\hspace{0.17em}\text{ln}\hspace{0.17em}{\widehat{f}}_n^{(1)}\left(y_i\right|{\widehat{\mathit{\beta }}}^{(1)})-\frac{1}{n}\sum _{i=1}^n\hspace{0.17em}\text{ln}\hspace{0.17em}{\widehat{f}}_n^{(1)}(y_i|\tilde{\mathit{\beta }})\stackrel{a.s.}{\to }0.$$(A.9)

For that purpose, define $$\begin{array}{c}B\equiv \Vert \underset{\mathit{\beta }\in \mathcal{C}}{\sup }\frac{1}{n}\sum _{i=1}^n\hspace{0.17em}\text{ln}\hspace{0.17em}{\widehat{f}}_n^{(1)}(y_i|\mathit{\beta })-\underset{\mathit{\beta }\in \mathcal{B}}{\sup }\frac{1}{n}\sum _{i=1}^n\hspace{0.17em}\text{ln}\hspace{0.17em}{\widehat{f}}_n^{(1)}(y_i|\mathit{\beta })\Vert ,\\ B_1\equiv \Vert \underset{\mathit{\beta }\in \mathcal{C}}{\sup }\frac{1}{n}\sum _{i=1}^n\hspace{0.17em}\text{ln}\hspace{0.17em}{\widehat{f}}_n^{(1)}(y_i|\mathit{\beta })-\underset{\mathit{\beta }\in \mathcal{C}}{\sup }E[\text{ln}\hspace{0.17em}f(y_i|\mathit{\beta })]\Vert ,\\ B_2\equiv \Vert \underset{\mathit{\beta }\in \mathcal{C}}{\sup }E[\text{ln}\hspace{0.17em}f(y_i|\mathit{\beta })]-\underset{\mathit{\beta }\in \mathcal{B}}{\sup }E[\text{ln}\hspace{0.17em}f(y_i|\mathit{\beta })]\Vert ,B_3\equiv \mathrm{}\Vert \underset{\mathit{\beta }\in \mathcal{B}}{\sup }E[\text{ln}\hspace{0.17em}f(y_i|\mathit{\beta })]-\underset{\mathit{\beta }\in \mathcal{B}}{\sup }\frac{1}{n}\sum _{i=1}^n\hspace{0.17em}\text{ln}\hspace{0.17em}{\widehat{f}}_n^{(1)}(y_i|\mathit{\beta })\Vert .\end{array}$$

Again by the triangle inequality, $B\le B_1+B_2+B_3.$ From (A.7), it is easy to show that $B_1\stackrel{a.s.}{\to }0$ and $B_3\stackrel{a.s.}{\to }0.$ To show that $B_2\stackrel{a.s.}{\to }0,$ first observe that by conditions (i) and (ii) of Lemma A.1 and the strong law of large numbers, (A.10) $$\text{Pr}\left(\underset{n\to \infty }{\lim }\right[\frac{1}{n}\sum _{i=1}^n(y_i-{\mathbf{x}}_i^{\prime }{\mathit{\beta }}_0)]=0)=1.$$(A.10)

This implies $$\begin{array}{c}\text{Pr}(\underset{n\to \infty }{\lim }{\mathit{\beta }}_0\in \mathcal{C})=1\Rightarrow \text{Pr}(\underset{n\to \infty }{\lim }[\arg \underset{\mathit{\beta }\in \mathcal{B}}{\sup }E[\text{ln}\hspace{0.17em}f(y_i|\mathit{\beta })]]\in \mathcal{C})=1\\ \Rightarrow \text{Pr}(\underset{n\to \infty }{\lim }{B}_2=0)=1,\end{array}$$ and the last result implies, by definition of almost sure convergence, that $B_2\stackrel{a.s.}{\to }0.$ Hence, $B\stackrel{a.s.}{\to }0$ and the constraint does not affect the result.

Lastly, to show that the algorithm asymptotically converges to ${\mathit{\beta }}_0,$ remark that (A.8) implies by Lemma A.2 that $\underset{x\in \mathbb{R}}{\sup }\Vert {\widehat{f}}_n^{(2)}(x)-f(x)\Vert \stackrel{a.s.}{\to }0$ where ${\widehat{f}}_n^{(2)}(x)$ is the kernel density-based estimator of the residuals corresponding to ${\widehat{\mathit{\beta }}}^{(1)}.$ Thus, by identical reasoning, we obtain ${\widehat{\mathit{\beta }}}^{(u)}\stackrel{a.s.}{\to }{\mathit{\beta }}_0$ for $u=1,2,\dots .$ □

Remark 5.

Conditions (i)–(iv) of Theorem 2.1 are the regularity conditions that are necessary for the convergence of MLE under the true density. Thus, the only additional conditions imposed are those in Lemma A.1 of which condition (i) of zero mean goes without loss of generality in the context of linear regression as we can always adjust the intercept parameter in $\mathit{\beta }$ if the center of $f(·)$ is not zero. Condition (ii) of Lemma A.1 may be restrictive in some cases as it rules out, for instance, the t(v)-distribution with $1<v\le 2.$ However, in Section 4.2 we observed that M-KDRE performs well for t(2(2) $${\widehat{f}}_n^{(u)}(x)={(n{h}_n)}^{-1}\sum _{i=1}^{n}K(\frac{x-{\widehat{\epsilon }}_i^{(u-1)}}{h_n})$$(2) ). In fact, its performance is best of all considered estimators under that error distribution. Hence, the practical use of M-KDRE does not seem to be restricted to distributions with finite variance. Conditions (iii) and (iv) of Lemma A.1 are easy to verify in practice, and conditions (v)–(vii) are technical requirements on the kernel and bandwidth. Note that (v) is not satisfied by the Gaussian kernel since that kernel does not have bounded support. In practice, however, the Gaussian kernel entails a significant computational advantage.

Proof of Theorem 2.3

(sketch): For the case of q (q < p) with linear random or deterministic equality constraints, the proof of consistency and asymptotic distribution of $\widehat{\mathit{\beta }}$ can be based on results in Crowder (1984) and Osborne (2000). In particular, given these results it follows that ${\widehat{\mathit{\beta }}}^{(1)}$ is asymptotically normal and efficient. The only condition on the initial estimator ${\widehat{\mathit{\beta }}}^{(0)}$ is that $({\widehat{\mathit{\beta }}}^{(0)}-{\widehat{\mathit{\beta }}}_0)=O_p(n^{-1/2}).$ For ${\widehat{\mathit{\beta }}}^{(1)}$ this follows from the proof of (Yao and Zhao 2013, Thm. 2.1). Hence, all subsequent estimates ${\widehat{\mathit{\beta }}}^{(u)}$ $(u=2,3,\dots )$ also satisfy (8(8) $$\sqrt{n}({\widehat{\mathit{\beta }}}^{(u)}-{\mathit{\beta }}_0)\stackrel{d}{\to }N(0,\mathbf{U}{(\mathbf{U}\prime {\mathcal{I}}_{\beta \beta }\mathbf{U})}^{-1}\mathbf{U}\prime )$$(8) ). □

Proof of Theorem 2.4.

Under the Gaussian kernel, the linear constraint $c(\mathit{\beta })=0,$ and a full-kernel method, the M-step in (10(10) $$\begin{array}{c}{\widehat{\mathit{\beta }}}_{(k+1)}^{(u)}={\widehat{\mathit{\beta }}}_{\text{LSE}}-{\left(\sum _{i=1}^n{\mathbf{x}}_i{\mathbf{x}}_i^{\prime }\right)}^{-1}\sum _{i=1}^n\left({\mathbf{x}}_i\sum _{j=1}^n{p}_{ij,(k+1)}^{(u)}{\widehat{\epsilon }}_j^{(u-1)}\right)\\ \hspace{1em}+\frac{1}{n}{\left(\sum _{i=1}^n{\mathbf{x}}_i{\mathbf{x}}_i^{\prime }\right)}^{-1}\sum _{i=1}^n{\mathbf{x}}_i\left(\sum _{i=1}^n\sum _{j=1}^n{p}_{ij,(k+1)}^{(u)}{\widehat{\epsilon }}_j^{(u-1)}\right)\end{array}$$(10) ) becomes (A.11) $${\widehat{\mathit{\beta }}}_{(k+1)}^{(u)}=\arg {\min }_{\mathit{\beta }}\sum _{i=1}^n\sum _{j=1}^n(p_{ij,(k+1)}^{(u)}{(y_i-{\mathbf{x}}_i^{\prime }\mathit{\beta }-{\widehat{\epsilon }}_j^{(u-1)})}^2)\mathrm{ }\mathrm{s}.\mathrm{t}.\mathrm{ }\sum _{\mathrm{i}=1}^{\mathrm{n}}({\mathrm{y}}_{\mathrm{i}}-{\mathbf{x}}_{\mathrm{i}}^{\prime }\mathit{\beta })=0$$(A.11)

This can be solved by Lagrangian optimization. Define the Lagrangian $\mathcal{L}$ as (A.12) $$\mathcal{L}(\mathit{\beta },\lambda )=\sum _{i=1}^n\sum _{j=1}^n(p_{ij,(k+1)}^{(u)}{(y_i-{\mathbf{x}}_i^{\prime }\mathit{\beta }-{\widehat{\epsilon }}_j^{(u-1)})}^2)-\lambda \sum _{i=1}^n(y_i-{\mathbf{x}}_i^{\prime }\mathit{\beta })$$(A.12) with first-order conditions (A.13) $$\frac{\partial \mathcal{L}}{\partial \mathit{\beta }}=-2\sum _{i=1}^n\sum _{j=1}^n(p_{ij,(k+1)}^{(u)}{\mathbf{x}}_i(y_i-{\mathbf{x}}_i^{\prime }\mathit{\beta }-{\widehat{\epsilon }}_j^{(u-1)}))-\lambda \sum _{i=1}^n{\mathbf{x}}_i=0$$(A.13) (A.14) $$\frac{\partial \mathcal{L}}{\partial \lambda }=\sum _{i=1}^n(y_i-{\mathbf{x}}_i^{\prime }\mathit{\beta })=0.$$(A.14)

By setting $x_{i,1}\equiv 1$ $(i=1,\dots ,n),$ the first element of the first-order condition in (A.13(A.13) $$\frac{\partial \mathcal{L}}{\partial \mathit{\beta }}=-2\sum _{i=1}^n\sum _{j=1}^n(p_{ij,(k+1)}^{(u)}{\mathbf{x}}_i(y_i-{\mathbf{x}}_i^{\prime }\mathit{\beta }-{\widehat{\epsilon }}_j^{(u-1)}))-\lambda \sum _{i=1}^n{\mathbf{x}}_i=0$$(A.13) ) implies (A.15) $$\begin{array}{cc}\lambda =& -\frac{2}{n}\sum _{i=1}^n\sum _{j=1}^n\left(p_{ij,(k+1)}^{(u)}(y_i-{\mathbf{x}}_i^{\prime }\mathit{\beta }-{\widehat{\epsilon }}_j^{(u-1)})\right)\\ =& -\frac{2}{n}\sum _{i=1}^n(y_i-{\mathbf{x}}_i^{\prime }\mathit{\beta })\sum _{j=1}^n{p}_{ij,(k+1)}^{(u)}+\frac{2}{n}\sum _{i=1}^n\sum _{j=1}^n{p}_{ij,(k+1)}^{(m)}{\widehat{\epsilon }}_j^{(u-1)}\\ =& \frac{2}{n}\sum _{i=1}^n\sum _{j=1}^n{p}_{ij,(k+1)}^{(u)}{\widehat{\epsilon }}_j^{(u-1)},\end{array}$$(A.15) where the last equality follows from (A.14(A.14) $$\frac{\partial \mathcal{L}}{\partial \lambda }=\sum _{i=1}^n(y_i-{\mathbf{x}}_i^{\prime }\mathit{\beta })=0.$$(A.14) ). Then, by plugging λ in (A.13(A.13) $$\frac{\partial \mathcal{L}}{\partial \mathit{\beta }}=-2\sum _{i=1}^n\sum _{j=1}^n(p_{ij,(k+1)}^{(u)}{\mathbf{x}}_i(y_i-{\mathbf{x}}_i^{\prime }\mathit{\beta }-{\widehat{\epsilon }}_j^{(u-1)}))-\lambda \sum _{i=1}^n{\mathbf{x}}_i=0$$(A.13) ), rearranging terms and using ${\sum }_{j=1}^n{p}_{ij,(k+1)}^{(u)}=1,$ we obtain (A.16) $$\begin{array}{cc}{\widehat{\mathit{\beta }}}_{(k+1)}^{(u)}& ={\left(\sum _{i=1}^n{\mathbf{x}}_i{\mathbf{x}}_i^{\prime }\right)}^{-1}\sum _{i=1}^n{\mathbf{x}}_i{y}_i-{\left(\sum _{i=1}^n{\mathbf{x}}_i{\mathbf{x}}_i^{\prime }\right)}^{-1}\sum _{i=1}^n({\mathbf{x}}_i\sum _{j=1}^n{p}_{ij,(k+1)}^{(u)}{\widehat{\epsilon }}_j^{(u-1)})\\ & \hspace{1em}+\frac{1}{n}{\left(\sum _{i=1}^n{\mathbf{x}}_i{\mathbf{x}}_i^{\prime }\right)}^{-1}\sum _{i=1}^n{\mathbf{x}}_i(\sum _{i=1}^n\sum _{j=1}^n{p}_{ij,(k+1)}^{(u)}{\widehat{\epsilon }}_j^{(u-1)}).\end{array}$$(A.16)

Recognize that the first term is equal to ${\widehat{\mathit{\beta }}}_{\text{LSE}}.$ Then, the fact that (9(9) $$p_{ij,(k+1)}^{(u)}=\frac{\hspace{0.17em}\text{exp}\hspace{0.17em}\{-\frac{1}{2h_n^2}{(y_i-{\mathbf{x}}_i^{\prime }{\widehat{\mathit{\beta }}}_{(k)}^{(u)}-{\widehat{\epsilon }}_j^{(u-1)})}^2\}}{\sum _{\mathcal{\ell }=1}^n\hspace{0.17em}\text{exp}\hspace{0.17em}\{-\frac{1}{2h_n^2}{(y_i-{\mathbf{x}}_i^{\prime }{\widehat{\mathit{\beta }}}_{(k)}^{(u)}-{\widehat{\epsilon }}_{\mathcal{\ell }}^{(u-1)})}^2\}}$$(9) ) and (10(10) $$\begin{array}{c}{\widehat{\mathit{\beta }}}_{(k+1)}^{(u)}={\widehat{\mathit{\beta }}}_{\text{LSE}}-{\left(\sum _{i=1}^n{\mathbf{x}}_i{\mathbf{x}}_i^{\prime }\right)}^{-1}\sum _{i=1}^n\left({\mathbf{x}}_i\sum _{j=1}^n{p}_{ij,(k+1)}^{(u)}{\widehat{\epsilon }}_j^{(u-1)}\right)\\ \hspace{1em}+\frac{1}{n}{\left(\sum _{i=1}^n{\mathbf{x}}_i{\mathbf{x}}_i^{\prime }\right)}^{-1}\sum _{i=1}^n{\mathbf{x}}_i\left(\sum _{i=1}^n\sum _{j=1}^n{p}_{ij,(k+1)}^{(u)}{\widehat{\epsilon }}_j^{(u-1)}\right)\end{array}$$(10) ) are the E- and M-step, respectively, of an EM type algorithm follows trivially from the proof of Theorem 2.2 in Yao and Zhao (2013). □