EconBase
← Back to paper

Endogenous Heteroskedasticity in Linear Models

Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.

61,897 characters · 12 sections · 48 citation commands

Rendered from LaTeX for readability, not typeset faithfully. Citation keys are highlighted; maths is left as source; figures, tables and equation environments are summarised rather than reproduced; unrecognised commands are greyed out so nothing is silently dropped. Email addresses are removed.

Endogenous Heteroskedasticity in Linear Models

\normalem

\thispagestyle{empty}

abstract\singlespacing Linear regressions with endogeneity are widely used to estimate causal effects. This paper studies a framework that involves two common practical issues: endogeneity of the regressors and heteroskedasticity that depends on endogenous regressors, i.e., endogenous heteroskedasticity. To address the inconsistency of the two-stage least squares estimator in this scenario, and recover the causal parameters of interest, we develop a framework for practical estimation and inference based on the control function approach allowing for discrete and continuous regressors. In particular, we suggest a simple two-step estimation procedure. We establish the limiting properties of the estimator, namely, consistency and asymptotic normality. In addition, we develop practical valid inference methods by proposing an estimator for the asymptotic variance-covariance matrix, and formally establishing its consistency. Monte Carlo simulations provide evidence on the finite-sample performance of the proposed methods and evaluate different implementation strategies. We revisit an empirical application on job training to illustrate the methods.

Keywords: Two-stage least squares, instrumental variables, heteroskedasticity, control function.

JEL: C21, C26

\pagenumbering{arabic}

Introduction

This paper studies estimation of linear models with endogeneity and heteroskedasticity. At the estimation stage, endogeneity is typically addressed using instrumental variables (IV) and two-stage least squares (2SLS) procedures. Heteroskedasticity, on the other hand, is usually handled with robust inference techniques. However, when heteroskedasticity is caused by an endogenous regressor, a situation referred to as endogenous heteroskedasticity, the 2SLS estimator may become inconsistent because the necessary exogeneity condition fails. The literature on econometric models that allow for endogenous heteroskedasticity is limited (see, e.g., FlorensHeckmanMeghirVytlacil08, ChenKhan14, and AbrevayaXu23).

In this paper, we focus on this issue and develop methods for estimation and inference of parameters in linear triangular models with endogeneity and endogenous heteroskedasticity that employ a control function (CF) approach and allow for continuous and discrete regressors. Endogenous heteroskedasticity models can be viewed as a special case of nonseparable models, see, e.g., Chesher03, ChernozhukovImbensNewey07, FlorensHeckmanMeghirVytlacil08, ImbensNewey09 and Jun09. These general methods usually employ nonparametric modeling and estimation together with a CF. Differently from these papers, our interest lies in practical linear regression models in the presence of multiple covariates, such that nonparametric methods become impractical. Thus, we rely on the functional necessity of building a feasible parametric model using the CF approach.\footnote{There is a large literature on the CF approach. See, e.g., among many others, NeweyPowellVella99, KimPetrin22, and han2025setvaluedcontrolfunctions. Different alternatives have been studied in the literature. See the survey by BlundellPowell08, and ImbensNewey09 for nonparametric and semiparametric models.} We start by building on the general identification results in FlorensHeckmanMeghirVytlacil08 and tailor the conditions to accommodate continuous and discrete variables such that we are able to write the parameter of interest as an explicit function of observables, and thus derive the estimator of interest. The key restrictions for identification are a CF condition with a flexible parameterization, along with another parameterization of the skedastic function in the structural equation. Given a valid CF constructed from the first stage, and the structure of the skedastic function, the parameter of interest in the second stage can be identified using a standard ordinary least squares (OLS) regression. We highlight that our proposed approach includes the 2SLS as a special case when there is no endogenous heteroskedasticity and there is no correction for the first-stage heteroskedasticity, a condition that depends on the CF assumption.

Estimation is implemented in a two-step procedure. First, the skedastic function in the first step can be estimated using flexible scale function methods suggested in RomanoWolf17. Given the estimate of the skedastic function, the CF can be computed as the normalized errors from the first-stage regression. In the second step, one performs a simple linear regression with a CF approach. Empirically, this step uses an OLS regression of the dependent variable on the endogenous, exogenous, and CF variables augmented with their interacting terms. We establish the limiting statistical properties of the proposed two-step estimator. Mild sufficient conditions are provided for the estimator to have desired asymptotic properties, namely, consistency and asymptotic normality. The limiting theory presented here incorporates elements of both the CF and generated regressors.

In addition, we develop practical statistical inference procedures by suggesting an estimator for the asymptotic variance-covariance matrix with the generated regressors. We formally establish the consistency of the estimator of the variance-covariance matrix. The proposed methods provide foundation for general inference procedures. Testing for general linear hypotheses is easily accommodated by Wald type tests, and non-linear hypothesis are readily available by employing the Delta method.

We also evaluate the finite sample performance of the proposed methods using Monte Carlo exercises. Numerical simulations confirm that the presence of heteroskedasticity in the structural equation induces bias in the 2SLS estimator. Moreover, heteroskedasticity in the first-stage equation may amplify the bias. The proposed CF method is able to produce estimates that concentrate around the true value of the target parameter. In addition, results improve when the sample size increases.

Finally, we present an empirical application that illustrates the framework discussed in this paper. We apply the estimator to study the effect of public-sponsored training program Job Training Partnership Act (JTPA) on future earnings. We estimate the effect of interest using the CF, and for comparison, the 2SLS. Empirical results show evidence that the 2SLS-IV estimator (the typical application in other papers that used the JTPA data) is substantially smaller than the alternative CF methods.

The specific literature considering the effects of endogenous heteroskedasticity in linear models is restricted. When treatment effects are heterogeneous among observationally identical individuals, that is, in the presence of heteroskedasticity in the structural equation, causal inference for policy evaluation is known to be more difficult. FlorensHeckmanMeghirVytlacil08 use the CF approach to identify the average treatment effects and the effect of treatment on the treated in models with a continuous endogenous regressor whose impact is heterogeneous. AbrevayaXu23 show that when a binary treatment is endogenous and there is endogenous heteroskedasticity, that is, the endogenous treatment indicator also affects the scale of the outcome variable, then standard methods for estimation of average treatment effects fail. In particular, they consider the standard IV estimator with binary treatment and show that it is an inconsistent estimator for the average treatment effect. After nonparametric identification is established, estimators are provided under a linear specification for the mean and variance treatment effect, as well as the average treatment effect on the treated. In the context of local average treatment effects, ChenKhan14 discuss identification and estimation of the heteroskedasticity term under different treatment statuses. As opposed to the papers listed above, we do not require particular support restrictions for the endogenous variables nor the instrumental variables. That is, we do not require continuity as in FlorensHeckmanMeghirVytlacil08, nor discreteness as in AbrevayaXu23. Heteroskedasticity can sometimes aid in addressing endogeneity problems rather than exacerbating them, as shown in Rigobon03, KleinVella10, and Lewbell12. However, this is not the case considered in the present study.

This paper is organized as follows. Section (ref) presents the model and shows the inconsistency of 2SLS under endogenous heteroskedasticity. In Section (ref) we propose a CF approach and discuss identification. Estimation and inference procedures are detailed in Section (ref). A Monte Carlo study is provided in Section (ref). Section (ref) illustrates the methods with an empirical application. Finally, Section (ref) concludes. Mathematical proofs are relegated to the Appendix.

A linear model with endogenous heteroskedasticity

In this section, we first describe the model of interest, which contains endogenous heteroskedasticity. Second, we motivate the proposed methods by reviewing the inconsistency of the two-stage least squares in this scenario.

Model

We consider the following model:

align[align omitted — 142 chars of source]

where the main parameter of interest is the scalar \(\alpha_1\). The vector of exogenous variables \(X\) is \((p_x - 1)\)-dimensional (including a constant), and the instrumental variable (IV) \(Z\) is \(p_z\)-dimensional. The variables \(\varepsilon\) and \(V\) are unobservable and their correlation makes \(D\) potentially endogenous.

The model in equations (ref) and (ref) has several important features. Despite the linearity in $D$ in (ref), it allows for heterogeneous effects on $Y$ due to the term $g(D, X)\varepsilon$. Because of this, the parameter of interest has an interpretation as an average effect. Indeed, if we define the potential outcome for unit $i$ when $D=d$ as $Y_{di} = d\alpha_1 + X_i'\alpha_2 + g(d, X_i)\varepsilon_i$, then $E[Y_{di} - Y_{d'i}] = (d - d')\alpha_1$.\footnote{See Appendix (ref) for more about the interpretation of $\alpha_1$.} However, the exogeneity conditions \(E[\varepsilon | X, Z] = 0\) and \(E[V | X, Z] = 0\) are not sufficient to identify \(\alpha_1\) via 2SLS, since \(E[g(D, X)\varepsilon | X, Z]\) may not be zero. Since the conditional variance of the error term in \(\eqref{eq:structural0}\) is \(g(D, X)^2 Var[\varepsilon|D,X]\), we refer to it as “endogenous heteroskedasticity,” a terminology borrowed from AbrevayaXu23.

The model above is a particular case of nonadditive triangular models, where the interest lies in estimating the average linear effect, $\alpha_1$. This specification would allow us to study the role of the skedastic function $g$ in the 2SLS bias for $\alpha_1$ and specific solutions using the control function approach. A comparable model is used in AbrevayaXu23 to study the role of endogenous heteroskedasticity in treatment effects models. Furthermore, equation (ref) could be rewritten as a general model $D=m(Z,X,V)$ but we use a similar skedastic structure to match equation (ref). The parametric model in the first-stage is mainly used due to the curse of dimensionality problem, since most empirical work uses multiple control variables. In a nonparemetric case, estimation could rely on methods suggested by JiSuXiao15 and LintonXiao19.

Motivation: Inconsistency of the 2SLS estimator

The 2SLS is a popular method to address endogeneity in linear models. In this section we show that the 2SLS estimator estimates the average effect $\alpha_1$ with a bias in the presence of endogeneous heteroskedasticity. Although the result here is simple, intuitive, and similar issues were highlighted in the literature, we include this section to motivate the linear representation above.

The failure of 2SLS to identify $\alpha_1$ in the presence of endogenous heteroskedasticity was highlighted by AbrevayaXu23 in the case of a binary \(D\). When a binary IV is available, they offer a detailed analysis of the bias and propose a novel identification and estimation strategy. In a similar context with a continuous \(D\), FlorensHeckmanMeghirVytlacil08 argue that the usual exogeneity IV requirements are not sufficient for identification.

The following example illustrates the bias in the 2SLS procedure in a model described by (ref)--(ref).

exampleConsider the following simple linear case of model (ref)--(ref) with one endogenous variable, one instrument, $g(D)=1+\delta D$ and $h(Z)=1+Z\gamma$: \begin{align*} Y &=D\alpha_1 + (1+D\delta ) \varepsilon\\ D &=Z\pi_1 + (1+Z\gamma)V. \end{align*} Here, $E[D\varepsilon]\neq0$ through the possible correlation between $\varepsilon$ and $V$. We assume that $E[\varepsilon|Z]=E[V|Z]=0$. The 2SLS estimator of $\alpha_1$ is the sample counterpart of $\frac{ Cov[Z,Y]}{ Cov[Z,D]}$. Hence, we have that \begin{align*} Cov[Z,Y] &=\alpha_1Cov[Z,D] + Cov[Z,\varepsilon] + \delta Cov[Z, D\varepsilon]. \end{align*} Now, $E[\varepsilon|Z]=0$ implies $Cov[Z,\varepsilon] =0$. Therefore, we have \begin{align*} \frac{ Cov[Z,Y]}{ Cov[Z,D]} &=\alpha_1 + \delta \frac{Cov[Z, D\varepsilon]}{ Cov[Z,D]}. \end{align*} The exogeneity assumptions on $Z$ do not imply that $Cov[Z, D\varepsilon]=0$. Therefore, $\alpha_1$ is not identified by the 2SLS moment conditions whenever there is heteroskedasticity in the structural equation, i.e., $\delta\neq 0.$ Under the stronger independence assumption: $Z\perp \! \! \! \perp (V,\varepsilon)$ (typical in experimental settings when $Z$ is random), the expression simplifies to \begin{align*} \frac{ Cov[Z,Y]}{ Cov[Z,D]} &=\alpha_1 + \delta\gamma \frac{E[V\varepsilon]}{\pi_1} \end{align*} In this case, the joint heteroskedasticity together with the endogeneity (in the form of $E[V\varepsilon]\neq 0$) is inducing the bias. Thus, in this case, one could have a situation where there is heterogeneity in the structural equation with $\delta \neq 0$, but if there is no heteroskedasticity in the first-stage, $\gamma=0$, then 2SLS is able to correctly identify the parameter of interest.

\color{black}

Next, we provide a formal result deriving the asymptotic bias in the 2SLS estimator under the endogenous heteroskedasticity. This not a novel results per se, but we specify it to our particular setting.

We maintain the following assumptions.

assumption$(i)$ The sample ${(Y_{i},D_{i},X_{i}',Z_i')}_{i=1}^n$ is i.i.d. following the model (ref)--(ref). $(ii)$ The following exogeneity conditions hold: $E[\varepsilon|X,Z]=E[V|X,Z]=0$.
assumptionThere is at least one component of $\pi_1$ which is not 0, and the matrices $E\begin{bmatrix} ZZ'&ZX'\\ XZ'&XX' \end{bmatrix}$ and $E\begin{bmatrix} XD&XX'\\ ZD&ZX' \end{bmatrix}$ have full column rank.

Assumptions (ref) and (ref) are standard in the literature allowing for endogeneity of $D$, the presence of observable controls $X$, as well as for the valid instruments $Z$ to be (intrinsically) exogenous. While the instruments $Z$ are uncorrelated with $\varepsilon$ and $V$, we will see that when we consider the unobservable in (ref) to be $g(D, X)\varepsilon$, then it might no longer be the case that $E[Zg(D, X)\varepsilon]= 0$. That is, despite $Z$ being uncorrelated to the unobserved innovations, and being relevant, the heteroskedastic nature of the triangular system can invalidate the 2SLS procedure. As an alternative, we propose a control function approach in Section (ref) below.

The next result compares the probability limit of the 2SLS estimator, denoted by $\overline{\alpha}_{1}$, to the target parameter $\alpha_1$. We assume that the data is generated according to the model in (ref)--(ref), and a researcher who is interested in $\alpha_1$, instruments $D$ with $Z$ in a 2SLS regression of $Y$ on $D$ and $X$ (which includes a constant). In order to deal with the presence of the additional controls $X$, we use a population version of the Frisch-Waugh-Lovell (FWL) theorem using linear projections.\footnote{See Section (ref) in the Appendix for a brief review of linear projections.}

lemmaConsider the model in (ref)--(ref). Let $\overline{\alpha}_{1}$ be the probability limit of the 2SLS estimator of $\alpha_1$ of a 2SLS regression of $Y$ on $D$ and $X$, with $D$ instrumented by $Z$ and $X$. Under Assumptions (ref) and (ref), we have \begin{equation} \overline{\alpha}_{1} = \alpha_1 + \Sigma_{h}^{-1} \pi_1' E[\overline Zg(Z'\pi_1 + X'\pi_2+ h(X,Z)V, X)\varepsilon], \end{equation} where $\Sigma_h:=\pi_1'E[\overline Z\, \overline Z']\pi_1$, and $\overline Z' := Z'-X'E[XX']^{-1}E[XZ']$.

We provide some further remarks to understand the intuition and nature of the result in Lemma (ref).

enumerate• The asymptotic bias of the 2SLS estimator in equation (ref) can be succinctly written as $\Sigma_{h}^{-1} \pi_1' E[\overline Zg(D, X)\varepsilon]$, but this might obscure the fact that $h(\cdot)$ is also playing a role through $D$, and could, in principle, amplify the bias. • If $g(\cdot)\equiv 1$, i.e., there is no heteroskedasticity in the structural equation, then there is no asymptotic bias in the 2SLS estimator because $E[\overline Z {g(D, X)\varepsilon} ]=E[\overline Z {\varepsilon} ]=0$ by Assumption (ref). The lack of bias is regardless of $h(\cdot)$. For a general $g(\cdot)$, the form of $h(\cdot)$ may matter. • If $h(\cdot)\equiv 1$, i.e., there is no heteroskedasticity in the first-stage, then the asymptotic bias depends only on the form of $g(\cdot)$. • The strength of the instrumental variable, through the magnitude of $\pi_1$, can alleviate the bias. To see this, suppose that $d_z=1$. Then, the asymptotic bias is inversely related to $\pi_1$: \begin{align*} \overline{\alpha}_{1} = \alpha_1 + \frac{E[\overline Z {g(D, X)\varepsilon}]}{\pi_1E[\overline Z^2]} . \end{align*} Even in this simple case, the sign of the bias is not totally obvious. • If it is suspected that $b_0\leq E[\overline Zg(D, X)\varepsilon]\leq b_1$, then partial identification of $\alpha_1$ is possible. Otherwise, this can be used in a sensitivity analysis. • In this paper, we focus on the 2SLS estimator. Nevertheless, we highlight the same result is valid for a GMM estimator. The main intuition is that, under endogenous heteroskedasticity, the moment condition used for estimation would be misspecified trough the function $g(\cdot)$.

A control function approach

This section discusses sufficient conditions, based on the control function (CF) approach, that allow for identification of the parameters of interest in the model given in equations (ref)--(ref), in the presence of both endogeneity and general heteroskedasticity. In the next section, we will suggest a practical estimator and develop inference procedures. Identification in additive separable models with endogenous heteroskedasticity using a CF strategy was initially proposed FlorensHeckmanMeghirVytlacil08. Our analysis builds on their results, with the key distinction between our approach and theirs being a trade-off between making a stochastic polynomial assumption on the CF in our approach versus their approach allowing for a general nonparametric function. The advantage of our strategy is reflected in the practicality of the estimator. We will show that, by specializing the CF to be a polynomial, the identification result shows that the parameter of interest can be written as a simple explicit function of the data, which in turn can be used for estimation and inference. Moreover, our approach allows the 2SLS-IV estimator to be nested within our model, and thus being a special case.

We maintain the following assumption.

assumptionControl Function Condition: $E(\varepsilon | D,X,Z) = E(\varepsilon |V)= r(V)$.

Assumption (ref) is a standard control function condition, see, e.g., BlundellPowell08. This is the same condition as Assumption A-3 in FlorensHeckmanMeghirVytlacil08. Under this condition, the conditional expectation of the response variable in equation (ref) can be written as

align[align omitted — 156 chars of source]

It should be noted that, from the previous display, the skedastic term in the structural equation does not depend only on the control function, but also on $g$. Further assumptions are required to identify $\alpha_1$. To illustrate this, we use Example (ref), where we assume that both $D$ and $Z$ are continuous.

exampleConsider the model in Example (ref) without covariates, and assume Assumption (ref) holds. Then (ref) becomes $E[Y|D,Z] = D\alpha_1 + g(D) r(\psi(D,Z))$, where $ V:=\psi(D,Z)$ is identified from the first stage. Define $f(D,Z):=E[Y|D,Z]$. If $D$ and $Z$ are continuous, consider the following derivatives \begin{align*} \frac{\partial f}{\partial D} &=\alpha_1+g_{D}r+g r_{ V}\psi_{D} \\ \frac{\partial f}{\partial Z} &=g r_{ V}\psi_{Z}, \end{align*} where the derivatives are denoted by subscripts, and we omit the arguments in each function. The left hand-sides are identifiable. By multiplying the second equation of the above display by $\psi_{Z}'(\psi_{Z}\psi_{Z}')^{-1}$, and substituting the result into the first equation, we obtain \begin{align*} \frac{\partial f}{\partial D} &=\alpha_1+g_{D}r+\frac{\partial f}{\partial Z}\psi_{Z}'(\psi_{Z}\psi_{Z}')^{-1}\psi_{D}. \end{align*} Note that, even though $\psi_{Z}$ and $\psi_{D}$ are identified, one cannot identify $\alpha_1$ because of the additional term $g_{D}r$, that is, \begin{align} \frac{\partial f}{\partial D}-\frac{\partial f}{\partial Z}\psi_{Z}'(\psi_{Z}\psi_{Z}')^{-1}\psi_{D} &= \alpha_1+g_{D}r. \end{align} While the left-hand side of (ref) contains only identified objects, the right-hand side contains three unknowns, $\alpha_{1}$, $g_{D}$ and $r$.

The previous example shows that further assumptions are required for identification of the parameter $\alpha_{1}$. Assumption (ref) below imposes a polynomial function on the skedastic function.

assumptionPolynomial Heteroskedasticity: $g(D,X)$ is a polynomial of degree $k_g$ in $D$ and $X$. That is $g(D,X)=1+\theta_1^{d}D+X'\theta^{x}_{1}+...+\theta^{d}_{k_{g}}D^{k_{g}}+ X'^{k_{g}}\theta^{x}_{k_{g}}$.

This condition is the same as in FlorensHeckmanMeghirVytlacil08, except that here it includes additional exogenous covariates.\footnote{FlorensHeckmanMeghirVytlacil08 write Assumption (ref) directly in their model in equations (1) and (2).} In Assumption (ref), $X'^{k_{g}}$ refers to the element-by-element $k_{g}$ power of $X'$. This could be replaced by any other polynomial that may include interactions.

Finally, we impose the following polynomial structure on the CF.

assumptionPolynomial Control Function: $r(V)$ is a polynomial of degree $k_v$ in $ V$ that does not contain a constant. That is $r( V)=\sum_{j=1}^{k_v} \mu_j {V}^j$.

This assumption is a trade-off between flexible parameterization and the nonparametric model used in FlorensHeckmanMeghirVytlacil08. Finally, the identification results below also require a full rank condition that implies that the functions in Assumptions (ref) and (ref) are not linearly dependent with the full set of covariates $(D,X)$. Assumption (ref) achieves this by not including a constant term in the function $r(V)$, but different conditions can be adapted as long as the full rank is satisfied, i.e. we may allow for $r(V)$ to include a constant, but restrict $g(D, X)$ to not contain a linear term with $D$ and $X$. The intuition for these restrictions is simple: we cannot have the same covariate to appear both in the linear function and in the implied control function expanded with the skedastic function. Suppose that there is only one endogenous regressor $D$ and that both the skedastic function and the CF contain a constant as: $g(D)=1+\delta D$ and $r(V)=1+V$. Then, $g(D)r( V)=1+ V+\delta D+\delta D V$, and in this case, one is not able to identify $\alpha_{1}$ because of the term $\delta D$.

Identification of the parameter of interest requires identification of $V$, which in turn requires identification of $h$, which we can achieved up-to-scale. See for instance, RomanoWolf17 and Wooldridge12 for parametric specifications, and JiSuXiao15 and LintonXiao19 for nonparametric specifications. The next assumption imposes a normalization on the second moment of the control function, and requires the skedastic function to be positive.

assumptionEither $h\equiv 1$ or $E[V^2 | X,Z]=1$ and $h(X,Z)>0$ almost surely.

The next lemma considers the identification of $\alpha_{1}$. It shows that the parameter of interest is identified by a regression of $Y$ on transformations of $(D,X,V)$. Define the linear projection of $D$ onto $W$ at $w$ as $L_{[D|W]}(w) := w' E[WW']^{-1}E[WD]$, where $W$ denotes all the regressors in equation (ref) except $D$, i.e. this contains $X'$ and the interactions of $\{1,D^s,X^{'s}\}_{s=1,...,k_g}$ with $\{{V}^{j}\}_{j=1,...,k_v}$.

lemmaConsider the regression \begin{align} E[Y|D, X, Z] & = D\alpha_1 + X'\alpha_2 + \left( 1 + \sum_{s=1}^{k_{g}} (\theta_{s}^{d}D^{s}+X'^{s}\theta_{s}^{x}) \right) \left( \sum_{j=1}^{k_v} \mu_j{V}^j \right) \nonumber \\ & =D\alpha_1 + W'\alpha_w. \end{align} Under Assumptions (ref)--(ref), $\alpha_1$ is the coefficient of $D$ in (ref) and is given by \begin{align} \alpha_1 &= \frac{E[(D- L_{[D|W]}(W))Y]}{E[(D- L_{[D|W]}(W))^2]}, \end{align} provided $E[(D- L_{[D|W]}(W))^2]\neq 0$.

Note that if $D$ is a dummy variable, as in the empirical application, then there is no need to consider a polynomial model of it and the regression is simplified. The identification result in Lemma (ref) above shows that the parameter of interest $\alpha_{1}$ in equation (ref) can be written as a function of observable data. In the next section, we propose a simple two-step estimator for it.

Estimation and inference

This section proposes a practical estimator for the parameter of interest $\alpha_{1}$, and establishes its statistical limiting properties. We consider an estimator that consists of a two-step procedure. The overall procedure is as following. In the first-stage, one estimates the control function (CF), $ V$. In the second-stage, one performs a simple OLS regression of the dependent variable $Y$ on the regressors $(D,X')'$ and the estimated CF, $\widehat{V}$, as well as on their polynomial interactions. In addition, we provide practical inference procedures.

First stage estimator: obtaining $\hat V$

We employ OLS in the first-stage and regress the endogenous variable $D$ on the instruments, $Z$, and exogenous variables $X$, to compute the residuals

align*[align* omitted — 84 chars of source]

Next, we compute the skedastic function $h(\cdot)$. In practice, to estimate the skedastic function we employ the parametric models in RomanoWolf17 or Wooldridge12. For instance,

equation[equation omitted — 109 chars of source]

or

equation[equation omitted — 135 chars of source]

The coefficients of these models can be estimated using (possibly non-linear) OLS regressions with $\check{V}^{2}_{i}$ or $ln(\check{V}^{2}_{i})$ as the dependent variable and the absolute value of the instruments and exogenous covariates. The resulting estimator is $\widehat{h}_i= h(X_i,Z_i;\widehat{\gamma})$. Using the estimated skedastic function, $\widehat{h}$, we obtain the normalized errors as

equation*[equation* omitted — 70 chars of source]

The following assumption imposes a standard regularity condition on the skedastic function determining its smoothness and existence of moments of the conditional heteroskedasticity.

assumptionAssume that $h(X,Z)=h(X,Z;\gamma)$, a parametric function with parameter $\gamma$. (i) The map $\gamma\mapsto h(x,z;\gamma)$ is twice differentiable with gradient denoted by $ \nabla_{\gamma}h(x,z;\gamma)$; (ii) $E[\nabla_{\gamma}h(X_i,Z_i;\gamma)\nabla_{\gamma}h(X_i,Z_i;\gamma)'h(X_i,Z_i;\gamma)^2]$ has full rank; (iii) the vectors $Z_i$, $X_i$ and the components of $\nabla_{\gamma}h(X_i,Z_i;\gamma)\nabla_{\gamma}h(X_i,Z_i;\gamma)'h(X_i,Z_i;\gamma)^2$ have finite second moments; and (iv) the vectors $Z_ih(X_i,Z_i)V_i$, $X_ih(X_i,Z_i)V_i$ and $ \nabla_{\gamma}h(X_i,Z_i;\gamma) \left(V_i^2-1\right)h(X_i,Z_i;\gamma)^3$ have finite second moments.

Denote $\phi:=(\pi_1', \pi_2', \gamma')'$. The next result is the asymptotic normality of $\widehat \phi=(\hat\pi_1', \hat\pi_2', \hat\gamma')'$. The precise form of the influence function $\psi_{\phi,i}$ is detailed in the appendix.

propositionUnder Assumptions (ref), (ref), (ref), and (ref), \begin{align*} \sqrt n(\widehat \phi-\phi)=\frac{1}{\sqrt n}\sum_{i=1}^n\psi_{\phi,i}+o_p(1)\overset{d}{\to}\mathcal N(0,\Omega_\phi), \end{align*} where $\Omega_\phi=E[\psi_{\phi}\psi_{\phi}']$.

Second stage estimator: augmented control function estimator

In the second step, the goal is to estimate $\alpha_1$ given in equation (ref). Recall that from the identification results discussed in the previous section, the population parameter is given by

align*[align* omitted — 77 chars of source]

where $W$ is defined as all the regressors on the right hand side of conditional average equation in (ref), but $D$.

Denote the preliminary estimator of $ V$ denoted by $\widehat{V}$. The estimator is given by

align*[align* omitted — 169 chars of source]

where $\widehat{L}[D|\widehat{W}_{i}] = \widehat{W}_{i}'(\widehat{W}'\widehat{W})^{-1}\widehat{W}' D.$ Here $\widehat{W}_{i}$ denotes the vector $W_{i}$ with $ V_i$ replaced by $\widehat{V}_i$.

Given that the standard CF model (i.e. adding $\check{V}$, the first-stage residuals, to the regression equation) is in fact the 2SLS estimator, this determines that the latter is a special case of our augmented CF procedure. That is, 2SLS-IV estimator is nested in our model. The presence of heteroskedasticity in the structural model and its subsequent effect can be evaluated by adding the interaction term $\check{V} D$ and checking its statistical significance. If the interaction is not statistically significant, this provides some evidence that the 2SLS might not be biased. This has a similar intuition to a test for endogeneity in linear regression models when checking the statistical significance of $\check{V}$ in the CF-2SLS model Wooldridge10. As such, our proposed method is robust to the presence of homoskedasticity and heteroskedasticity in the structural equation. However, the inclusion of additional terms and the previous estimation of the CF in the reduced form OLS estimation might impact the efficiency of the augmented CF procedure.

For notational brevity, we write equation (ref) as

equation[equation omitted — 71 chars of source]

where $W(\phi)$ contains all the regressors except $D$ as given in Lemma (ref), and by construction, $E[U|D,W]=0$. Furthermore, $\alpha=(\alpha_1,\alpha_w')'$. The next result is the asymptotic normality of $\widehat \phi=(\hat\pi_1', \hat\pi_2', \hat\gamma')'$.

The next theorem contains the limiting distribution of the proposed two-step CF estimator of $\alpha=(\alpha_1,\alpha_w')'$, denoted $\widehat\alpha=(\widehat\alpha_1,\widehat\alpha_w')'$. This provides details for the asymptotics for the estimator of interest, i.e., $\alpha_1$, which is an application of OLS with generated regressors. We consider the following assumption:

assumptionLet $W(\phi)$ be all the regressors in (ref) (except $D$) as given in Lemma (ref), and define $R_i(\phi):=(D_i,W_i(\phi)')'$. Assume that uniformly in $\phi$ (i) $E[R_i(\phi)R_i'(\phi)]$ is non-singular; (ii) a (uniform) law of large numbers holds for $\frac{1}{n}\sum_{i=1}^n R_i(\phi)R_i'(\phi)$, $\frac{1}{n}\sum_{i=1}^n(U_iJ_{\phi,i}- R_i(\phi)\alpha'J_{\phi,i})$, where $J_\phi$ is the Jacobian of $\phi\mapsto R(\phi)$, $J_\phi=\nabla_\phi R(\phi)$; and (iii) $E[\|RR'U^2\|]<\infty$.

Assumption (ref) $(i)$ ensures that the CF additional terms to the structural equation are not perfectly collinear to the main regressors, and Assumptions (ref) $(ii)$ and $(iii)$ are required for the asymptotic representations using the influence functions.\footnote{Sufficient conditions for a uniform law of large numbers can be found in Lemma 2.4 in NeweyMcFadden1994. Essentially, it depends on a dominating function for $R(\phi)$ or, alternatively, $h(x,z;\gamma)$.}

theoremUnder Assumptions (ref)-(ref), \begin{align*} \sqrt n (\widehat \alpha-\alpha) &= \frac{1}{\sqrt n}\sum_{i=1}^n\psi_{\alpha,i} + o_p(1)\overset{d}{\to}\mathcal N(0, \Omega_\alpha) \end{align*} where $\psi_{\alpha,i}:=E[RR']^{-1}R_i U_i + E[RR']^{-1}E[UJ_\phi- R(\phi)\alpha'J_\phi]\psi_{\phi,i} $, and $\Omega_\alpha=E[\psi_{\alpha}\psi_{\alpha}'].$

Note that $\psi_{\phi,i}$ is contributing to $\psi_{\alpha,i}$, the influence function of $\hat\alpha$. From the above result, we obtain $\sqrt n(\widehat \alpha_1 - \alpha_1)\stackrel{d}{\rightarrow} \mathcal N (0,\Omega_{\alpha,(1,1)}),$ where $\Omega_{\alpha,(1,1)}$ is the $(1,1)$ element of $\Omega_{\alpha}$.

Inference

Now we turn our attention to inference in the model under endogenous heteroskedasticity, and suggest a practical estimator for the asymptotic variance-covariance matrix. We also formally establish its consistency.

Estimation of $\Omega_{\alpha,(1,1)}$ can be carried out by resorting to the sample counterpart of $\Omega_{\alpha}$ as following:

align*[align* omitted — 106 chars of source]

with

align*[align* omitted — 339 chars of source]

where

align*[align* omitted — 732 chars of source]

Next, we state a result establishing the consistency of the above estimator.

lemmaIf Assumptions (ref) and (ref) hold, $E[\|R\|^4]<\infty$, $E[\|J_{\phi}\|^4]<\infty$, $E[U^4]<\infty$, $E[\|Z\|^4]<\infty$, $E[\|X\|^4]<\infty$,$E[|h(X,Z;\gamma)V|^4]<\infty$, $c_i(\phi):=\nabla_{\gamma}h(X_i,Z_i;\gamma) \left( V_i^2-1\right)h(X_i,Z_i;\gamma)^3$ is differentiable in $\phi$ with derivative $\nabla_\phi c_i(\phi)$, $E\big[\|\nabla_\phi c_i(\phi)\|^2\big] < \infty$, and satisfies a uniform law of large numbers, then $\widehat \Omega_{\alpha}\overset{p}{\to}\Omega_{\alpha}.$

Given the result in Lemma (ref) and Theorem (ref), inference is standard. A $(1-\theta)\%$ confidence interval for $\alpha_1$ can be constructed as $\widehat \alpha_1\pm Z_{\theta/2}\sqrt {\widehat \Omega_{\alpha,(1,1)}}$, where $\widehat \Omega_{\alpha,(1,1)}$ is the $(1,1)$ element of $\widehat \Omega_{\alpha}$, and $Z_{\theta/2}$ is the corresponding critical value ($\theta/2$ area to the left) from the standard normal distribution.

Moreover, given the asymptotic normality in Theorem (ref), general hypotheses on the parameter $\alpha_{1}$, other more generally $\alpha$, can be accommodated by simple Wald-type tests. The Wald process and associated limiting theory provide a natural foundation for the hypothesis $R\alpha =r$ when $r$ is known, with a Chi-square limiting distribution. Non-linear testing can be easily performed by employing the Delta method.

Monte Carlo Simulations

This section provides numerical simulations to assess the finite sample performance of the proposed methods. We use the following version of the model in equations (ref) and (ref),

align[align omitted — 204 chars of source]

This allows for endogeneity and different types of heteroskedasticity in both the first- and second-stages. The parameter $\lambda$ allows for $D$ being exogenous (i.e. $\lambda=0$) or endogenous (i.e. $\lambda\neq0$) in the structural equation (ref).

We compute simulation results for the cases of $U\sim N(0,1)$, $V\sim N(0,1)$, $Z$ following a folded normal distribution, that is $Z\sim |N(0,1)|$. Additionally, they are all independent of each other. We set $X\equiv 1$ so that all models have a constant term. Moreover, we use sample sizes of $n\in\{250,500,1000\}$. The number of replications is set to $2000$.

We consider different types of heteroskedastic models given by

align[align omitted — 144 chars of source]

Parameters $(\delta_1,\delta_2)$ allow for heteroskedasticity of the endogenous variable in the structural equation (i.e. $\delta_1\neq0$ or $\delta_2\neq0$) and nonlinear effects (i.e. $\delta_2\neq0$), while the parameter $\delta_{3}$ allows for heteroskedasticity of the exogenous variable, and in this case to have a constant term. Since the model includes a constant this ensures that the variance is not zero under the absence of endogenous heteroskedasticity. The parameter $\gamma_1$ also allows for $D$ having a conditional homoskedasticity (i.e. $\gamma_1=0$) or heteroskedasticity (i.e. $\gamma_1\neq0$) in the first-stage equation (ref). The parameters are $\alpha_1=\alpha_2=\pi_1=\pi_2=\delta_3=\gamma_2=1$, $\gamma_1\in\{0,1\}$, $\delta_1\in\{0,1\}$ and $\delta_2\in\{0,0.2\}$, thus covering a wide range of scenarios.

We present results for different estimators: OLS, standard 2SLS (note that this is equivalent to a standard CF model where $\check V$ is included as an additional regressor in the structural equation), and the proposed CF with one interaction term CF1 $(\widehat{V} ,\widehat{V} D)$, as well as with two interaction terms CF2 $(\widehat{V} ,\widehat{V} D,\widehat{V} D^2)$.

For the CF models we consider one type of skedastic function given in RomanoWolf17, as given in equation (ref), $h(X,Z;\gamma)= (|Z|\gamma_{1}+|X|\gamma_{2})^{1/2}=(\gamma_0+|Z|\gamma_{1})^{1/2}$, such that $\hat{h}=h(X,Z;\hat\gamma)$. Let $\check{V}$ be the residuals of the first-stage regression of $D$ on $Z$ and a constant, then $h^2$ is estimated as the regression of $\check{V}^2$ on $|Z|$ and a constant. Note that this is not the correctly specified skedastic function given in equation (ref).

Tables (ref) and (ref) report the simulation results for the case of $\lambda=1$, i.e. endogeneity, for $\gamma_1=0$ and $\gamma_1=1$, respectively. For completeness Tables (ref) and (ref) report the simulation results for the case of $\lambda=0$, i.e. no endogeneity, also for $\gamma_1=0$ and $\gamma_1=1$, respectively. In all cases we report empirical bias and variance of the estimators OLS, 2SLS and CF estimators, and for the latter we compute the average of the estimated variance and 95% empirical coverage using the Gaussian confidence intervals.

The results in Tables (ref) and (ref) clearly indicate that OLS is biased and that 2SLS is also biased in the presence of endogenous heteroskedasticity (i.e. $\delta_1\neq0$ and/or $\delta_2\neq 0$). Moreover, note that the parameter $\gamma_1$, i.e. heteroskedasticity in the first-stage, aggravates the 2SLS bias. However, the proposed CF estimators provide unbiased results. Note that this depends on the functional form of the heteroskedasticity-inducing function $g$, but not on correctly specifying $h$. When $\delta_2=0$, both CF1 and CF2 works, but when $\delta_2=0.2$, only CF2 completely eliminates the bias. Overall, this indicates that the CF approach eliminates the bias, but nonlinear terms may be needed for different functional forms of the heteroskedasticity inducing functions $g$.

The last two columns of the tables provide the average estimated variance and the 95% coverage of the CF estimators. In all cases, the estimated variance gets close to the simulated variance and it is very similar when $n=1000$. Moreover, the empirical coverage is close to the 95% whenever we have an approximately unbiased estimator.

Tables (ref) and (ref) indicate that when there is no endogeneity and no correction is needed from OLS. the CF estimator is unbiased and behaves in a similar fashion to the 2SLS counterpart.

table[table omitted — 2,046 chars of source]
table[table omitted — 1,958 chars of source]
table[table omitted — 1,966 chars of source]
table[table omitted — 1,967 chars of source]

Empirical application

We apply the estimator to the study of the effect of public-sponsored training programs. As argued in LaLonde95, public programs of training and employment are designed to improve participant's productive skills, which in turn would increase their earning potential and decrease dependency on social welfare benefits. We use data from the Job Training Partnership Act (JTPA) program, that has been extensively studied in the literature. For example, see Bloometal97 for a description, and AbadieAngristImbens02 for applications. The JTPA was a large publicly-funded training program that began funding in October 1983 and continued until late 1990's. We focus on the Title II subprogram, which was offered only to individuals with “barriers to employment” (long-term use of welfare, being a high-school drop-out, 15 or more recent weeks of unemployment, limited English proficiency, physical or mental disability, reading proficiency below 7th grade level or an arrest record). Individuals in the randomly assigned JTPA treatment group were offered training, while those in the control group were excluded for a period of 18 months. Our interest lies in measuring the effect of actual training on of participants' future earnings.

We use the database in AbadieAngristImbens02 that contains information about adult male and female JTPA participants and non-participants. Let $Z$ denote the indicator variable for those receiving a JTPA offer. Of those offered, 60% completed the training; of those in the control group completion rate was less than 2%. The dependent variable $Y$ is the logarithm of 30 month accumulated earnings (we exclude individuals without earnings), $Z$ is a dummy variable for the JTPA offer, $D$ is the endogenous dummy variable corresponding to the JTPA training, $X$ is a set of exogenous covariates containing individual characteristics. The variables used as detailed in the database are: sex hsorged black hispanic married wkless13 afdc age2225 age2629 age3035 age3644 age4554. Note that there are then many covariates that will not allow us to implement a nonparametric strategy, and this serves for motivation for our parametric implementation. The parameter of interest is the effect of JTPA training on earnings.

table[table omitted — 2,073 chars of source]

Results for different models appear in Table (ref). Column $earnings$ uses earnings as dependent variable, while $earnings>0$ uses only observations with non-zero earnings, same as the $log(earnings)$ column. Besides OLS and 2SLS, we also consider two control function methods. The first model (CF0) does not consider a stochastic function in the first-stage and uses only the first-stage residuals $\check V$. The second model (CF1) applies a first-stage skedastic function using a regression of $\check{V}^2$ on the instrument and a constant ($\check{V}^2=\gamma_0+Z\gamma_{1}+error$ as in the Monte Carlo section) with $\widehat V=\check V/h(X,Z;\hat\gamma)$ and $h(X,Z;\gamma)=(\gamma_0+Z\gamma_{1})^{1/2}$. Both CF models add the estimated $V$ from the first-stage (either $\check V$ for CF0 or $\hat V$ for CF1) and the interaction with $D$ (either $\check V\times D$ for CF0 or $\hat V\times D$ for CF1). All models include other covariates $X$ but coefficients are not reported.

The empirical results in Table (ref) show that the OLS point estimates are larger than the 2SLS estimates for all models in levels and logarithms. However, the CF models have $\alpha_1$ estimates closer to OLS but larger than the IV-2SLS estimators. CF0 estimates have very large standard errors for $\alpha_1$ when compared with all other models, and the estimates of $\alpha_{\check{V}D}$ are not statistically significant. Moreover, CF1 have standard errors that are approximately half the CF0 level, and of similar magnitude of 2SLS ones.

The empirical set-up is similar to that in Example (ref) for the case of the IV being independent of the error components ($Z$, JTPA offer, is random). For convenience we copy the main result here:

align*[align* omitted — 119 chars of source]

Note that the mentioned example suggests that the 2SLS estimate is biased if there is both first-stage ($\gamma$) and second-stage endogenous heteroskedasticity ($\delta$). For the former, the CF1 point estimates of the skedastic first-stage model reveal that $\gamma_1>0$ (i.e. $\gamma$ in the Example) is statistically significant. For the latter, CF1 also shows that $\alpha_{\hat{V}D}<0$ (i.e. $\delta$ in the Example) is also significant, thus indicating endogenous heteroskedasticity. If we assume that $E[V\varepsilon]>0$, i.e. individuals with higher unobservable components in the earnings equation are more likely to select into training, then the 2SLS bias is negative (provided that $\delta<0$, $\gamma>0$, $\pi_1>0$). As such, 2SLS might be downward biased to estimate $\alpha_1$. The CF1 estimates are in line with this argument.

Overall, the CF results show empirical evidence that the job training is effective in increasing earning, but estimates are more modest than those from simple OLS regressions, but larger than 2SLS estimates.

Conclusion

This paper shows that the estimation of linear models with endogeneous heteroskedasticity using the popular 2SLS method may not be estimating the parameter of interest. As an alternative we provide a simple to implement augmented control function estimator. The CF approach provides an intuitive modeling strategy to handle endogeneity and lies in the center of general approaches to this particular set-up. We show that this can be adapted to the least-squares linear estimation methods to correct for the 2SLS bias. In fact, the proposed solution is to augment the standard CF model (which corresponds to 2SLS) with interactions of the structural equation variables with the residuals from the first-stage.