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.
107,792 characters · 16 sections · 60 citation commands
Inference for Two-Stage Extremum EstimatorsFor comments and suggestions, we are grateful to Arnaud Dufays, Ulrich Hounyo, Mathieu Marcoux, Antoine Djogbenou, Frank Windmeijer, Xiaohong Chen, Jean-Marie Dufour, Jad Beyhum, Prosper Dovonon, Désiré Kédagni, Pamela Giustinelli and Florian Pelgrin. We also thank the participants of the EDHEX Business School econometric seminar, the CIREQ econometric seminar, the 58th Annual Meetings of the CEA, and the 2024 conference of IAAE. Replication codes for the results from this research are available at https://github.com/ahoundetoungan/InferenceTSE.
{4pt} {4pt}
Two-stage estimation approaches are widely used to address challenges such as endogeneity, selection bias, non-identification, missing data, and high dimensionality hirano2003efficient, jofre2003estimation, newey2003instrumental, freyberger2022identification, chernozhukov2022locally, ichimura2022influence, boucher2020estimating. These methods involve estimating a function or parameter in the first stage (FS), followed by incorporating the estimator into a second-stage (SS) model to estimate another parameter, denoted as $\boldsymbol{\theta}_0$. The estimator of $\boldsymbol{\theta}_0$ is referred to as a two-stage (TS) or plug-in estimator. For inference, researchers often require the asymptotic distribution of $\sqrt{n}(\boldsymbol{\hat{\theta}}_n - \boldsymbol{\theta}_0)$, where $\boldsymbol{\hat{\theta}}_n$ is the TS estimator and $n$ is the sample size. However, deriving this distribution can be complex due to the sampling error introduced in the FS.
Three important issues may arise regarding the asymptotic distribution of $\sqrt{n}(\boldsymbol{\hat \theta}_n - \boldsymbol{\theta}_0)$:
While the literature has often focused on the first issue newey1994large, ackerberg2012practical, chen2015sieve, the latter two have not received much attention, despite their potential to arise in many contexts. Standard inference methods typically impose conditions for the asymptotic distribution of $\sqrt{n}(\boldsymbol{\hat \theta}_n - \boldsymbol{\theta}_0)$ to be normal with a mean of zero. If these conditions are not met, such methods may become inappropriate and lead to incorrect conclusions in hypothesis tests.
In this paper, we propose a simulation-based inference approach for TS estimators, focusing on extremum estimators in the SS. We accommodate a broad range of FS estimators, including extremum estimators, high-dimensional estimators, and other types of estimators (e.g., Bayesian estimators). The novelty of our approach lies in its ability to estimate the asymptotic cumulative distribution function (CDF) of $\sqrt{n}(\boldsymbol{\hat \theta}_n - \boldsymbol{\theta}_0)$, even when the distribution is non-normal and the asymptotic mean is different from zero. By estimating this mean, we demean $\sqrt{n}(\boldsymbol{\hat \theta}_n - \boldsymbol{\theta}_0)$ and introduce a debiased estimator that demonstrates improved finite sample properties. Our debiased estimator is straightforward and easily applicable to complex models. Unlike resampling methods, it avoids the need for repeated calculations of the TS estimator.
Regularization bias can arise in various situations, particularly when the FS model is high-dimensional or when the number of observations in the FS is small relative to $n$. The slow convergence rate of the FS estimator can significantly distort the SS estimator chernozhukov2017double, belloni2014high, belloni2017program, cattaneo2019two. Moreover, when dealing with complex models, researchers may compute a Bayesian estimator, such as a posterior mean, which is then used in the SS to obtain a standard extremum estimator breza2020using, lubold2023identifying. As noted by zellner1984bayesian, Bayesian estimators may not be asymptotically normally distributed, potentially leading to a non-normal asymptotic distribution in the SS. Non-normal asymptotic distributions can also occur in the FS with certain frequentist estimators that are robust to outliers, such as the trimmed mean and U-statistics stigler1973asymptotic, ma2011asymptotic. Our approach covers a broad range of models, including these contexts where standard methods are no longer applicable.
Our approach consists first of examining the conditional asymptotic distribution of \break$\sqrt{n}(\boldsymbol{\hat{\theta}}_n - \boldsymbol{\theta}_0)$, given any realization of the FS estimator. Since we are disregarding the FS sampling error at this stage, the conditional distribution resembles that of a single-step estimator. Consequently, under certain assumptions, it will be normal with some mean and variance that depend on the FS estimator. Next, we characterize the unconditional asymptotic CDF of $\sqrt{n}(\boldsymbol{\hat{\theta}}_n - \boldsymbol{\theta}_0)$ by using both the asymptotic normality at the SS, conditional on the FS estimator, and the asymptotic distribution of the FS estimator. We demonstrate that this CDF may not follow a normal distribution and may have a mean that is not necessarily zero.
We then combine simulations from an estimator of the distribution of the FS estimator with the asymptotic normality in the second stage, conditional on the FS estimator, to simulate the asymptotic distribution of $\sqrt{n}(\boldsymbol{\hat{\theta}}_n - \boldsymbol{\theta}_0)$. This simulation allows us to construct a sample for $\sqrt{n}(\boldsymbol{\hat{\theta}}_n - \boldsymbol{\theta}_0)$, which we use to approximate its CDF and obtain confidence intervals (CIs) for $\boldsymbol{\theta}_0$. Additionally, the average of the constructed sample estimates the asymptotic mean of $\sqrt{n}(\boldsymbol{\hat{\theta}}_n - \boldsymbol{\theta}_0)$, which helps reduce the bias of $\boldsymbol{\hat{\theta}}_n$.
Simulating from the estimated distribution of the FS estimator is a key requirement of our method and is feasible for a broad range of models. When conducting inference in the FS using certain frequentist approaches, the asymptotic distribution is known (generally normal), which allows us to draw from a normal distribution. In cases where a Bayesian estimator is used in the FS, the posterior distribution can serve as an estimator, and samples can be obtained via Gibbs sampling or the Metropolis-Hastings algorithm casella1992explaining, chib1995understanding.\footnote{An estimator of the distribution of the FS estimator can also be derived using resampling methods when standard approaches are not easily applicable. In general, our approach avoids the need for multiple computations of the estimators, except when resampling is required to obtain the asymptotic distribution in the FS.}
We evaluate the performance of our method through an extensive simulation study, covering various models, including an instrumental variable (IV) model with many weak instruments cattaneo2019two, a Poisson model with unobserved variables, and a multivariate GARCH model gonccalves2022bootstrapping, with a number of returns increasing in $n$. Our simulation results confirm that our method performs well on these models.
Furthermore, to demonstrate the effectiveness of our method, we revisit the empirical analysis by fortin2015peer on peer effects in adolescent fast-food consumption. The estimation of endogenous peer effects often relies on IV methods, where instruments are constructed based on the characteristics of friends' friends bramoulle2009identification. We address the issue of biased estimates that can arise when instruments are weak. By using the characteristics of both close and distant friends, we expand the IV set to include many weak instruments. We correct the finite sample bias of the IV estimator and provide valid inference. Our findings suggest that a one-point increase in the average friend's fast-food consumption frequency leads to a 0.23 increase in one's fast-food consumption frequency.
This paper contributes to the extensive and growing literature on inference methods for sequential estimators. Most existing studies impose conditions for the asymptotic distribution of $\sqrt{n}(\boldsymbol{\hat \theta}_n - \boldsymbol{\theta}_0)$ to be normal with a zero mean newey1984method, andrews1994asymptotics, hotz1993conditional, murphy2002estimation, boucher2020estimating, houndetoungan2023ident. Examples include cases where both the first- and second-stage estimators are finite-dimensional extremum estimators converging at the same rate, and scenarios where the SS estimator is asymptotically invariant to infinitesimal variations in the FS estimator. We contribute to this literature by proposing a new method that relies on more flexible conditions. Additionally, we do not impose a specific class of estimators in the FS, making our approach more general. Importantly, our method yields results comparable to those obtained from classical asymptotic inference methods when applicable.
Even though the asymptotic distribution of the plug-in estimator is normal, accurately calculating its variance can be complex. Various methods have been proposed in the literature to estimate the asymptotic variance newey1984method, newey1994asymptotic, ackerberg2012practical, chen2015sieve, and resampling methods have also been employed for this purpose gonccalves2005bootstrap, kline2012score, armstrong2014fast, honore2017poor, gonccalves2022bootstrapping. We discuss scenarios where our simulation method can be used to estimate the asymptotic variance from the estimated CDF of $\sqrt{n}(\boldsymbol{\hat \theta}_n - \boldsymbol{\theta}_0)$. Moreover, our method not only estimates the asymptotic variance but also the asymptotic mean, thereby reducing the bias of $\boldsymbol{\hat \theta}_n$.
Our framework is also related to the literature on extremum estimators in the presence of high-dimensional nuisance parameters ai2007estimation, belloni2014high, belloni2017program, farrell2015robust, chernozhukov2015valid, mikusheva2022inference. Although $\boldsymbol{\hat \theta}_n$ can still be consistent in this context, the limiting distribution of $\sqrt{n}(\boldsymbol{\hat \theta}_n - \boldsymbol{\theta}_0)$ may not have a zero mean. Resampling methods are often used to approximate the regularization bias of $\boldsymbol{\hat \theta}_n$, but the choice of method depends on the studied model. For instance, in IV approaches with many weak instruments, the standard bootstrap method fails to infer $\boldsymbol{\theta}_0$, while the Jackknife method is appropriate cattaneo2019two. The standard bootstrap also fails when the FS involves variable selection methods such as Lasso chatterjee2011bootstrapping. We contribute to this literature because our bias reduction technique can be applied to a broad range of models. As long as simulating from the distribution of the FS estimator is possible, the asymptotic mean of $\sqrt{n}(\boldsymbol{\hat \theta}_n - \boldsymbol{\theta}_0)$ can be estimated to reduce the bias of $\boldsymbol{\hat \theta}_n$.
The regularization bias of extremum estimators in the presence of high-dimensional nuisance parameters can also be addressed using the double or debiased machine learning (DML) approach chernozhukov2017double, chernozhukov2018double, chernozhukov2022locally. This method combines Neyman-orthogonal moments (or scores) with cross-fitting to produce an estimator whose bias converges to zero faster than $1/\sqrt{n}$. chernozhukov2017double present this approach in the context of several models, including partially linear regression models and interactive regression models. However, obtaining orthogonal moments can be challenging in some cases, particularly for nonlinear models or those with weakly dependent data. Our approach differs from the DML method in that it does not require modifying moment functions or obtaining orthogonal scores. Instead, we retain the standard estimation procedure and apply bias reduction post-estimation. We illustrate the bias-reducing performance of our method using a multivariate GARCH model, where obtaining an orthogonal score can be challenging.
\paragraph{Plan of the Paper} The remainder of the paper is organized as follows. In Section (ref), we present our framework. Section (ref) provides an overview of our approach using a leading example. In Section (ref), we present our main results. Section (ref) provides a simulation study to assess the finite sample performance of our approach. In Section (ref), we present an empirical analysis with peer effects. Section (ref) concludes the paper.
\paragraph{Notation} The symbols $\mathbb{E}$ and $\mathbb{V}$ denote expectation and variance, respectively. $\lVert . \rVert$ is the $\ell_2$-norm. $\partial_{x}$ is the derivative with respect to $x$. The symbol $\mathbbm{1}\{\cdot\}$ is the indicator function. If $\boldsymbol a = (a_1,~ \dots, ~ a_d)^{\prime}$, $\boldsymbol b=(b_1,~ \dots, ~ b_d)^{\prime}$ are vectors in $\mathbb{R}^d$, then $\boldsymbol a \preceq \boldsymbol b$ means that $a_k \leq b_k$ for all $k$. $\boldsymbol{I}_{d}$ is the $d$-dimensional identity matrix. We use $\lim$ and $\operatorname{plim}$ to denote the standard limit and the limit in probability as $n$ grows to infinity, respectively. For a positive definite matrix $\mathbf{M}$, we use $\mathbf{M}^{1/2}$ to denote its Cholesky decomposition and $\mathbf{M}^{-1/2}$ to denote the Cholesky decomposition of its inverse.
This section introduces the class of plug-in estimators that are studied in this paper. For expositional ease, we consider the case of M-estimators in the second stage (SS). However, our findings can be extended to any extremum estimator since inference methods for extremum estimators are similar, regardless of whether it is an M-estimator, GMM estimator, or MD estimator amemiya1985advanced. Moreover, due to this similarity, we interchangeably use the terms M-estimator and extremum estimator in the paper.
In the SS, we assume that the practitioner maximizes an objective function given by
where $\mathbf{y}_n = (y_{1},~\dots,~y_{n})^{\prime}$, $\mathbf{X}_n = (\boldsymbol{x}_{1},~\dots,~\boldsymbol{x}_{n})^{\prime}$, $\mathbf{\hat B}_n = (\boldsymbol{\hat\beta}_{n,1},~\dots,~\boldsymbol{\hat\beta}_{n,n})^{\prime}$, and $q$ is a known function. In Equation (ref), $y_{i}$ and $\boldsymbol{x}_{i}$ are observed variables for the $i$-th unit in the sample (e.g., $y_{i}$ is a dependent variable and $\boldsymbol{x}_{i}$ are explanatory variables). $\boldsymbol{\hat\beta}_{n,1},~\dots,~\boldsymbol{\hat\beta}_{n,n}$ are estimators from some first stage (FS) regression (e.g., prediction of some variable in a preliminary regression). The estimator $\boldsymbol{\hat\beta}_{n,i}$ may be a scalar or finite-dimensional vector. The subscript $i$ may also refer to time in time-series models. We will refer to $\mathbf{\hat B}_n$ as the FS estimator. Importantly, we do not require $\boldsymbol{\hat\beta}_{n,i}$ to originate from an extremum estimation, or to have a particular asymptotic distribution (like the normal distribution). However, we assume that $\boldsymbol{\hat\beta}_{n,i}$ uniformly converges in probability to some $\boldsymbol{\beta}_{0,i}$, the true value of the parameter that it is designed to estimate (see Assumption (ref) below).
Let $\boldsymbol{\hat{\theta}}_n$ be the estimator that maximizes the objective function (ref). $\boldsymbol{\hat{\theta}}_n$ is called two-stage (TS) or plug-in estimator. We denote by $\boldsymbol{\theta}_0$ the true value of the parameter $\boldsymbol{\theta}$; i.e., the value taken by $\boldsymbol{\theta}$ in the data-generating process (DGP).
Special cases within our framework arise when $\boldsymbol{\beta}_{0,i} = f(\boldsymbol{z}_i, \boldsymbol{\gamma}_0)$ for some function $f$, where $\boldsymbol{z}_i$ is a control variable that may overlap components of $\boldsymbol{x}_i$, and $\boldsymbol{\gamma}_0$ is a parameter. In this case, we have $\boldsymbol{\hat\beta}_{n,i} = f(\boldsymbol{z}_i, \boldsymbol{\hat\gamma}_n)$, where $\boldsymbol{\hat\gamma}_n$ is an estimator of $\boldsymbol{\gamma}_0$. An example of this situation is the instrumental variable (IV) approach with $\boldsymbol{z}_i$ being the instrument and $\boldsymbol{\hat\beta}_{n,i}$ is the predicted value of the endogenous variable to be plugged into the second stage cattaneo2019two. The function $f$ may also be constant; that is, $\boldsymbol{\beta}_{0,i} = \boldsymbol{\beta}_{0}$ for any $i$, where $\boldsymbol{\beta}_{0}$ is a finite-dimensional vector to be estimated in the first stage murphy2002estimation. In the case of semiparametric or nonparametric specification, we can have $\boldsymbol{\beta}_{0,i} = f_n(\boldsymbol{z}_i, \boldsymbol{\gamma}_{0,n})$, where the specification of the function $f_n$ and the dimension of the parameter $\boldsymbol{\gamma}_{0,n}$ depends on the sample size $n$. Examples of estimators in this situation are power series, splines, and Fourier series approximations belloni2015some.
Let $\mathcal{Y}\subseteq \mathbb{R}^{K_{y}}$, $\mathcal{X}\subseteq \mathbb{R}^{K_{x}}$, and $\mathcal{B} \subseteq \mathbb{R}^{K_{\beta}}$ be the supports of $y$, $\boldsymbol{x}$, and $\boldsymbol{\beta}_{0,i}$, respectively, where $K_{y}$, $K_{x}$, and $K_{\beta}$ are the corresponding dimensions. Let also $\boldsymbol{\Theta} \subset \mathbb{R}^{K_{\theta}}$ be the space of $\boldsymbol{\theta}_0$, where $K_{\theta}$ is the dimension of $\boldsymbol{\theta}_0$. We introduce the following assumptions.
Assumption (ref) is a common requirement when dealing with FS estimators that may be infinite-dimensional chen2003estimation, ichimura2010characterization. The condition will hold in many applications. For the case where $\boldsymbol{\beta}_{0,i} = f(\boldsymbol{z}_i, \boldsymbol{\gamma}_0)$, Assumption (ref) requires $\boldsymbol{\hat\gamma}_n$ to be a consistent estimator and $f(\boldsymbol{z}_i, \boldsymbol{\gamma})$ to be continuously differentiable in $\boldsymbol{\gamma}$, with bounded derivative uniformly in $i$.\footnote{The result follows from the mean value theorem: $\boldsymbol{\hat{\beta}}_{n,i} - \boldsymbol{\beta}_{0,i} = \partial_{\boldsymbol{\gamma}^{\prime}}f(\boldsymbol{z}_i, \boldsymbol{\gamma}^+)(\boldsymbol{\hat\gamma}_n - \boldsymbol{\gamma}_0)$, for some $\boldsymbol{\gamma}^+$ that lies between $\boldsymbol{\hat\gamma}_n$ and $\boldsymbol{\gamma}_0$.} For nonparametric sieve estimators, lower-level conditions can be imposed to satisfy Assumption (ref) belloni2015some. Certain of these conditions are discussed by cattaneo2019two in their online appendix. Assumption (ref) also holds in the case where $\boldsymbol{\beta}_{0,i}$ represents fixed effects from another model dzemski2019empirical, yan2019statistical. Assumption (ref) sets regularity conditions on the objective function's behavior. These conditions are generally imposed for classical M-estimators and do not involve the FS estimator amemiya1985advanced.
Despite the risk of regularization bias, $\boldsymbol{\hat{\theta}}_n$ generally converges to $\boldsymbol{\theta}_0$ in probability, even when the FS estimator is high-dimensional chen2003estimation, cattaneo2019two. The regularization bias specifically affects $\sqrt{n}(\boldsymbol{\hat \theta}_n - \boldsymbol{\theta}_0)$, which may not have a zero mean asymptotically due to the $\sqrt{n}$ factor. We acknowledge the consistency of $\boldsymbol{\hat{\theta}}_n$ as a high-level assumption.
The proof of this consistency is context-dependent, and the required conditions may vary. In Online Appendix (OA) (ref), we discuss primitive conditions for Assumption (ref) by adapting Theorem 4.1.1 of amemiya1985advanced to our framework.
Before presenting the theoretical framework underlying our approach and the formal results, this section provides an overview through an illustrative example involving a latent variable model. In this example, $\sqrt{n}(\boldsymbol{\hat \theta}_n - \boldsymbol{\theta}_0)$ asymptotically follows a normal distribution with a mean of zero; this result can be established using standard methods newey1994large. This section demonstrates how our method can be effectively applied to a simple case before extending it to more complex scenarios.
We present our main results in this section. Technical details of proofs can be found in Appendix (ref). The first-order condition of the maximization of (ref) is $\frac{1}{n}\sum_{i = 1}^n\partial_{\boldsymbol{\theta}}q(\boldsymbol{\hat \theta}_n, ~y_{i}, ~\boldsymbol{x}_{i}, ~\boldsymbol{\hat{\beta}}_{n,i}) = 0$. By applying the mean value theorem to $\frac{1}{n}\sum_{i = 1}^n\partial_{\boldsymbol{\theta}}q(\boldsymbol{\theta}, ~y_{i}, ~\boldsymbol{x}_{i}, ~\boldsymbol{\hat{\beta}}_{n,i})$ with respect to (w.r.t) $\boldsymbol{\theta}$, we obtain:
where $\boldsymbol{\dot q}_{n,i}(y_i,~\boldsymbol{\hat{\beta}}_{n,i}) = \partial_{\boldsymbol{\theta}} q(\boldsymbol{\theta}_0, ~y_{i}, ~\boldsymbol{x}_{i}, ~\boldsymbol{\hat{\beta}}_{n,i})$ and $\mathbf{A}_n = -\frac{1}{n}\sum_{i = 1}^n \partial_{\boldsymbol{\theta}}\partial_{\boldsymbol{\theta}^\prime}q(\boldsymbol{\theta}^+_n, ~y_{i}, ~\boldsymbol{x}_{i}, ~\boldsymbol{\hat{\beta}}_{n,i})$, for some $\boldsymbol{\theta}^+_n$ that lies between $\boldsymbol{\hat \theta}_n$ and $\boldsymbol{\theta}_0$. As $\operatorname{plim} \boldsymbol{\hat \theta}_n = \boldsymbol{\theta}_0$, we also have $\operatorname{plim} \boldsymbol{\theta}^+_n = \boldsymbol{\theta}_0$. In large samples, $\mathbf{A}_n$ is assumed to be nonsingular (see Assumption (ref) below). Let $\boldsymbol{\dot q}_{n}(\mathbf{y}_n,~\mathbf{\hat{B}}_n) = \frac{1}{\sqrt{n}} \sum_{i = 1}^n \boldsymbol{\dot q}_{n,i}(y_i,~\boldsymbol{\hat{\beta}}_{n,i})$. We will refer to $\boldsymbol{\dot q}_{n}(\mathbf{y}_n,~\mathbf{\hat{B}}_n)$ as the influence function (IF).
In the case of a single-step estimator, the central limit theorem (CLT) implies (under regularity conditions) that the IF is asymptotically normally distributed with zero mean amemiya1985advanced. A crucial condition that is required by the CLT is that the dependence among the variables $\boldsymbol{\dot q}_{n,i}(y_i,~\boldsymbol{\hat{\beta}}_{n,i})$'s is "weak." Roughly speaking, if we define a certain order between the subscripts $i$'s (e.g., if $i$ is time), the correlation between $\boldsymbol{\dot q}_{n,i}(y_i,~\boldsymbol{\hat{\beta}}_{n,i})$ and $\boldsymbol{\dot q}_{n,j}(y_j,~\boldsymbol{\hat{\beta}}_{n,j})$ must vanish at a certain rate as $\lvert i - j\rvert$ grows to infinity withers1981central, romano2000more, ekstrom2014general. For TS estimators, the variables $\boldsymbol{\dot q}_{n,i}(y_i,~\boldsymbol{\hat{\beta}}_{n,i})$'s are dependent on each other because they all depend on the same FS estimator. Consequently, the weak dependence condition does not hold in general, even though $\boldsymbol{\hat{\beta}}_{n,i}$ uniformly converges in probability to $\boldsymbol{\beta}_{0, i}$. Without imposing additional conditions, there is no general CLT that guarantees asymptotic normality in this case.
Our approach does not involve applying the CLT directly to the IF. Instead, we assume that the conditional distribution of the IF, given $\mathbf{\hat{B}}_n$, is asymptotically normal. This condition is less restrictive in many cases because, conditional on $\mathbf{\hat{B}}_n$, the sampling error from the first stage is ignored, and the IF can be viewed as that of a single-step extremum estimator.
We introduce the following regularity assumptions.
Condition ((ref)) of Assumption (ref) imposes weak requirements on the existence of the conditional and unconditional moments of the IF. Condition ((ref)) will hold in general because $\mathbf{V}_n$ can be expressed as an average, whereas $\mathcal{E}_n$ is a sum of $n$ random variables scaled by $1/\sqrt{n}$. By the uniform Law of Large Numbers (LLN), if $\mathbf{V}_n$ is smooth in $\mathbf{\hat{B}}_n$, then it will converge in probability to a constant because the FS estimator is consistent (see Example (ref)). Additionally, in many cases, it is possible to express $\mathcal{E}_n$ as a function of $\mathbf{C}_n(\boldsymbol{\hat\gamma}_n - \boldsymbol{b}_n)$, for some sequences $\mathbf{C}_n$ and $\boldsymbol{b}_n$ and some estimator $\boldsymbol{\hat\gamma}_n$, such that $\mathbf{C}_n(\boldsymbol{\hat\gamma}_n - \boldsymbol{b}_n)$ has a limiting distribution.\footnote{In Example (ref), $\textstyle\mathcal{E}_n$ can be approximated using a first-order Taylor expansion around $\boldsymbol{\gamma}_0$ as $\mathcal{E}_n \approx - (\frac{2\theta_0}{n}\sum_{i = 1}^n\boldsymbol{z}_i^{\prime}\boldsymbol{\gamma}_0\boldsymbol{z}_i^{\prime})\sqrt{n}(\boldsymbol{\hat\gamma}_n - \boldsymbol{\gamma}_0)$. Consequently, $\mathcal{E}_0$ is normally distributed with a zero mean.} Importantly, Condition ((ref)) is flexible enough to encompass many models. Specifically, we allow for $\mathcal{E}_0$ not to be normally distributed and its mean may not be zero.
Assumption (ref) ensures the consistency of the Hessian matrix. Under some weak regularity conditions, for instance, if $\partial_{\boldsymbol{\theta}}\partial_{\boldsymbol{\theta}^\prime}q(\boldsymbol{\theta}_0, ~y_{i}, ~\boldsymbol{x}_{i}, ~\boldsymbol{\beta}_{0,i})$ is ergodic stationary across $i$ with a finite variance, the LLN implies that $\operatorname{plim}\frac{1}{n}\sum_{i = 1}^n \partial_{\boldsymbol{\theta}}\partial_{\boldsymbol{\theta}^\prime}q(\boldsymbol{\theta}_0, ~y_{i}, ~\boldsymbol{x}_{i}, ~\boldsymbol{\beta}_{0,i}) = \mathbf{A}_0$. Assumption (ref) extends this convergence to the case where $\boldsymbol{\theta}_0$ and $\boldsymbol{\beta}_{0,i}$ are replaced with consistent estimators. It adapts Conditions (B) in Theorem 4.1.3 of amemiya1985advanced to TS estimation approaches. We discuss primitive conditions for Assumption (ref) in OA (ref). These conditions require the Hessian at $\boldsymbol{\theta}_n^+$ to be smooth in $\boldsymbol{\theta}_n^+$ and $\boldsymbol{\hat{\beta}}_{n,i}$.
In the rest of this section, we first present some theoretical results on the asymptotic distribution of $\Delta_n$. Subsequently, we discuss the finite sample approximation of this distribution and introduce our debiased plug-in estimator.
We first examine the conditional distribution of the IF, given $\mathbf{\hat{B}}_n$. Indeed, treating $\mathbf{\hat{B}}_n$ as a predetermined sequence in $n$ allows us to approach the problem as in the case of single-step M-estimators. However, the IF may not have a zero mean, even asymptotically, thereby leading to a limiting distribution of $\Delta_n$ that is not centered at zero chernozhukov2018double. Therefore, we define the standardized IF as $\boldsymbol{u}_{n}(\mathbf{y}_n,~\mathbf{\hat{B}}_n) := \mathbf{V}_n^{-1/2}(\boldsymbol{\dot q}_{n}(\mathbf{y}_n,~\mathbf{\hat{B}}_n) - \mathcal{E}_n)$, which has a zero mean and variance $\boldsymbol{I}_{K_{\theta}}$. We impose the following assumption.
Conditional on $\mathbf{\hat{B}}_n$, the variable $\boldsymbol{u}_{n}(\mathbf{y}_n,~\mathbf{\hat{B}}_n)$ can be viewed as the standardized IF of a single-step estimator. Consequently, a conditional CLT may imply Assumption (ref) Rubshtein1996ACL.\footnote{A similar interpretation of Assumption (ref) by kato2011note is that ${\sup_{g\in LB}}\lvert \mathbb{E}[g(\boldsymbol{u}_{n}(\mathbf{y}_n,~\mathbf{\hat{B}}_n))|\mathbf{\hat{B}}_n] - {\displaystyle\varint}_{\mathbb{R}} g(t)d\Phi(\boldsymbol{t})\rvert = o_p(1)$, where $LB$ is the set of all functions on $\mathbb{R}$ with Lipschitz norm bounded by one. See also fligner1979use, van2000asymptotic who used conditional asymptotic normality.} For example, if $y_i$ is independent across $i$, conditional on $\mathbf{\hat{B}}_n$, then the variables $\boldsymbol{\dot q}_{n,i}(y_i,~\boldsymbol{\hat \beta}_{n,i})$'s would also be independent across $i$, conditional on $\mathbf{\hat{B}}_n$. Thus, the Lyapunov CLT or Lindeberg CLT, conditional on $\mathbf{\hat{B}}_n$, implies Assumption (ref) (under similar conditions to that of single-step estimators). When $y_i$'s are dependent, we may use a more general CLT for dependent processes if the dependence is weak conditional on $\mathbf{\hat{B}}_n$.
Assumption (ref) allows us to separate the model error in the SS (conditional on the FS sampling error) from the sampling error in the FS. Importantly, since $\boldsymbol{u}_{n}(\mathbf{y}_n,~\mathbf{\hat{B}}_n)$ is standardized, its first two conditional moments, given $\mathbf{\hat{B}}_n$, are independent of $\mathbf{\hat{B}}_n$. Consequently, even unconditionally, $\boldsymbol{u}_{n}(\mathbf{y}_n, \mathbf{\hat{B}}_n)$ asymptotically follows a standard normal distribution. We claim and show this result in Lemma (ref) in Appendix (ref). However, this result does not extend to the non-standardized IF, as its conditional moments depend on $\mathbf{\hat{B}}_n$.
The following theorem establishes the asymptotic distribution of $\Delta_n$.
The proof of Theorem (ref) is presented in Appendix (ref). Since $\Delta_n = \mathbf{A}_n^{-1}\mathbf{V}_n^{1/2}\boldsymbol{u}_{n}(\mathbf{y}_n,~\mathbf{\hat{B}}_n) + \mathbf{A}_n^{-1}\mathcal{E}_n$ and $\boldsymbol{u}_{n}(\mathbf{y}_n,~\mathbf{\hat{B}}_n)$ is asymptotically normally distributed, we show that $\boldsymbol{u}_{n}(\mathbf{y}_n,~\mathbf{\hat{B}}_n)$ can be substituted with $\boldsymbol{\zeta}$. It is noteworthy that this substitution is not trivial because $\boldsymbol{u}_{n}(\mathbf{y}_n,~\mathbf{\hat{B}}_n)$ and $\mathcal{E}_n$ may be correlated. However, we show that they are asymptotically independent; that is, both the asymptotic conditional distribution of $\boldsymbol{u}_{n}(\mathbf{y}_n,~\mathbf{\hat{B}}_n)$, given $\mathbf{\hat{B}}_n$ and the asymptotic unconditional distribution are the same (see Lemma (ref)). In the expression of $\boldsymbol{\psi}_n$, the sampling error from the FS is captured by the term $\mathbf{A}_0^{-1}\mathcal{E}_n$, whereas $\mathbf{A}_0^{-1}\mathbf{V}_0^{1/2}\boldsymbol{\zeta}$ captures variability in the SS.\footnote{In Appendix (ref), we extend Theorem (ref) to the uniform convergence. We show that $\sup_{\boldsymbol{t}\in\mathbb{R}^{K_{\theta}}}\lvert \mathbb{P}(\mathbf{V}_n^{-1/2}\mathbf{A}_0\Delta_n \preceq \boldsymbol{t}) - G(\boldsymbol{t})\rvert = o_p(1)$, where $G$ is the limiting distribution function of $\boldsymbol{\zeta} + \mathbf{V}_n^{-1/2}\mathcal{E}_n$.}
Similarly, we can also decompose the asymptotic variance of $\Delta_n$. Let $\boldsymbol{\Sigma}_n = \mathbb{V}(\boldsymbol{\dot q}_{n}(\mathbf{y}_n,~\mathbf{\hat{B}}_n))$ be the variance of the IF and $\boldsymbol{\Sigma}_0 = \lim \boldsymbol{\Sigma}_n$ be the limit of this variance. By Slutsky's theorem, the asymptotic variance of $\Delta_n $ is $\mathbb{V}(\textstyle\Delta_0):=\mathbf{A}_0^{-1}\boldsymbol{\Sigma}_0\mathbf{A}_0^{-1}$. This expression is similar to the asymptotic variance formula for single-step M-estimators. Yet, a notable difference here is that the sampling error from the FS estimator is incorporated into $\boldsymbol{\Sigma}_0$. By the law of iterated variances, we have $\boldsymbol{\Sigma}_n = \mathbb{E}(\mathbf{V}_n) + \mathbb{V}(\mathcal{E}_n)$. The first term on the right-hand side (RHS) is the asymptotic variance of the IF conditional on the FS sampling error. The second term is the variance due to the FS estimation. Using this equation, we establish the following result.
Furthermore, Theorem (ref) implies that $\Delta_n$ is asymptotically normally distributed if $\mathcal{E}_0$ follows a normal distribution. However, $\Delta_n$ may exhibit a regularization bias that depends on the expectation of $\mathcal{E}_0$. This leads to the following result.
The regularization bias of $\Delta_n$ is given by $\mathbf{A}_0^{-1}\mathbb{E}(\mathcal{E}_0)$. Corollary (ref) shares similarities with Theorem 1 of cattaneo2019two. In the context of IV approaches with many instruments, they show that $\mathbf{V}_n^{-1/2}\mathbf{A}_0(\Delta_n - \mathbf{A}_0^{-1}\mathbb{E}(\mathcal{E}_n))$ is asymptotically normally distributed, where $\mathbf{A}_0^{-1}\mathbb{E}(\mathcal{E}_n)$ represents the bias of $\Delta_n$. Corollary (ref) generalizes this result to a broad class of models.
In this section, we discuss how to simulate the asymptotic distribution of $\boldsymbol{\psi}_n = \mathbf{A}_0^{-1}\mathbf{V}_0^{1/2}\boldsymbol{\zeta} + \mathbf{A}_0^{-1}\mathcal{E}_n$ in finite samples. Since $\boldsymbol{\psi}_n$ and $\Delta_n$ share the same asymptotic distribution (Theorem (ref)), we can use the simulated distribution to infer $\boldsymbol{\theta}_0$.
We construct the sample $\mathcal{S}_{\kappa}=\{\boldsymbol{\hat \psi}_{n,s}: s = 1, ~\dots, ~ \kappa\}$ for some integer $\kappa \geq 1$, where $\boldsymbol{\hat \psi}_{n,1}$, \dots, $\boldsymbol{\hat \psi}_{n,\kappa}$ are independent variables with the same asymptotic distribution as $\boldsymbol{\psi}_n$. To obtain $\boldsymbol{\hat \psi}_{n,s}$, we replace the unknown nonstochastic quantities $\mathbf{A}_0$ and $\mathbf{V}_0$ in the expression of $\boldsymbol{\psi}_n$ with their estimators, and the random variables $\boldsymbol{\zeta}$ and $\mathcal{E}_n$ with independent draws from their (approximated) distributions. Specifically, $\boldsymbol{\hat \psi}_{n,s}$ is given by: $$\boldsymbol{\hat \psi}_{n,s} = \mathbf{\hat A}_n^{-1}\mathbf{\hat V}_n^{1/2}\boldsymbol{\zeta}_s + \mathbf{\hat A}_n^{-1}\hat{\mathcal{E}}_{n,s},$$ where $\mathbf{\hat A}_n$ and $\mathbf{\hat V}_n$ are respectively consistent estimators of $\mathbf A_0$ and $\mathbf{V}_0$, $\boldsymbol{\zeta}_1,~ \dots~\boldsymbol{\zeta}_{\kappa}$ are independent draws from $N(0, ~\boldsymbol{I}_{K_{\theta}})$, and $\hat{\mathcal{E}}_{n,1}$, \dots, $\hat{\mathcal{E}}_{n,\kappa}$ are independent draws from the (approximated) distribution of $\mathcal{E}_0$. By Assumption (ref), a consistent estimator of $\mathbf{A}_0$ is simply $\mathbf{\hat A}_n = -\frac{1}{n}\sum_{i = 1}^n \partial_{\boldsymbol{\theta}}\partial_{\boldsymbol{\theta}^\prime}q\big(\boldsymbol{\hat{\theta}}_n, ~y_{i}, ~\boldsymbol{x}_{i}, ~\boldsymbol{\hat{\beta}}_{n,i}\big)$. We will later discuss how to obtain $\mathbf{\hat V}_n$ and $\hat{\mathcal{E}}_{n,s}$.
The sample $\mathcal{S}_{\kappa}$ plays a crucial role. It can be used to construct confidence intervals for $\boldsymbol{\theta}_0$. Without loss of generality, assume that $\boldsymbol{\theta}_0$ is a scalar. Let $T_{\alpha}$ be the $\alpha$-quantile of the sample $\{\boldsymbol{\hat\theta}_n - \boldsymbol{\hat \psi}_{n,s}/\sqrt{n}: s = 1, ~\dots,~ \kappa\}$. Then, $[T_{\frac{\alpha}{2}}, ~ T_{1 - \frac{\alpha}{2}}]$ is a consistent estimator of the $(1 - \alpha)$ CI of $\boldsymbol{\theta}_0$, in the sense that $\lim_{\kappa\to\infty}\operatorname{plim} \mathbb{P}(\boldsymbol{\theta}_0 \in [T_{\frac{\alpha}{2}}, ~ T_{1 - \frac{\alpha}{2}}]) = 1 - \alpha$.\footnote{In practice, the integer $\kappa$ must be sufficiently large to ensure that the approximation error due to the number of simulations is negligible. Unlike resampling methods, increasing $\kappa$ does not introduce numerical issues, as we do not require many estimators of $\boldsymbol{\theta}_0$ from various samples.}
Furthermore, following Theorem (ref), we can obtain a consistent estimator of the asymptotic variance of $\Delta_n$, by replacing the variance of $\mathcal{E}_0$ with the sample variance of $\hat{\mathcal{E}}_{n,s}$. This leads to the following estimator for the asymptotic variance of $\boldsymbol{\hat{\theta}}_n$:
where $\boldsymbol{\hat \Sigma}_n^{\kappa} = \hat{\mathbf{V}}_{n} + \frac{1}{\kappa - 1} \sum_{s = 1}^{\kappa} (\hat{\mathcal{E}}_{n,s} - \boldsymbol{\hat \Omega}_{n}^{\kappa}) (\hat{\mathcal{E}}_{n,s} - \boldsymbol{\hat \Omega}_{n}^{\kappa})^{\prime}$ and $\boldsymbol{\hat \Omega}_{n}^{\kappa} = \frac{1}{\kappa}\sum_{s = 1}^{\kappa} \hat{\mathcal{E}}_{n,s}$.
To obtain $\mathbf{\hat V}_n$ and $\hat{\mathcal{E}}_{n,s}$, we first need to compute $\mathbf{V}_n$ and $\mathcal{E}_{n}$, which are the conditional variance and conditional expectation of the IF, given $\mathbf{\hat{B}}_n$. In general, computing $\mathcal{E}_{n}$ is straightforward by replacing $\mathbf{y}_n$ in $\boldsymbol{\dot q}_{n}(\mathbf{y}_n,~\mathbf{\hat{B}}_n)$ with its specification (see Example (ref)). In this exercise, exogenous variables, such as $\mathbf{X}_n$, can be treated as nonstochastic, as is often done in practice. This simplification is innocuous and analogous to inferring $\boldsymbol{\theta}_0$ conditional on $\mathbf X_n$. Even if $\mathcal{E}_n$ does not have a closed-form expression, we can employ a large-sample approximation. For instance, $\mathcal{E}_n$ can be approximated by $(1/\sqrt{n})\sum_{i = 1}^n \boldsymbol{\dot q}_{n,i}(y_i,~\boldsymbol{\hat{\beta}}_{n,i})$.
Two scenarios may arise regarding the conditional variance $\mathbf{V}_n$. First, $\mathbf{\hat{B}}_n$ and $\mathbf{y}_n$ may be independent, conditional on the exogenous variables in the model. This occurs when the error terms in both stages are independent. An example is when the first stage involves predicting unobserved exogenous variables using auxiliary models that are independent of the second stage's error terms chernozhukov2018double, breza2020using, lubold2023identifying, boucher2020estimating. In such cases, $\mathbf{V}_n$ can also be computed by substituting $\mathbf{y}_n$ with its specification. The calculations here are similar to those in a single-step approach. As for the $\mathcal{E}_{n}$, we can use a large-sample approximation if $\mathbf{V}_n$ does not have a closed form. For example, we can estimate $\mathbf{V}_n$ using a heteroskedasticity and autocorrelation consistent (HAC) covariance matrix andrews1991heteroskedasticity.
The second scenario, which is more challenging, arises when the error terms in both stages are not independent. This renders intricate the dependence between $\boldsymbol{\dot q}_{n,i}(y_i,~\boldsymbol{\hat{\beta}}_{n,i})$'s conditional on $\mathbf{\hat{B}}_n$, making it difficult to estimate $\mathbf V_n$. This situation is typical in IV approaches and estimation methods that use $\mathbf{y}_n$ in the first stage dufaysselective.
One important context where we can address the problem includes linear IV methods. The SS of this approach consists of regressing $\mathbf y_n$ on $\mathbf{\hat X}_n$, where $\mathbf{\hat X}_n$ is the prediction of $\mathbf{X}_n$ using an OLS regression of $\mathbf{X}_n$ on an instrument set. To avoid the complication of treating the correlation between $\mathbf{\hat X}_n$ and the error term of the SS, we instead regress $\mathbf{\hat y}_n$ on $\mathbf{\hat X}_n$ in the SS, where $\mathbf{\hat y}_n$ is the prediction of $\mathbf y_n$ using an OLS regression of $\mathbf y_n$ on the same instrument set. The IF becomes $\boldsymbol{\dot q}_{n}(\mathbf{\hat y}_n,~\mathbf{\hat{B}}_n)$ and the FS estimator includes both $\mathbf{\hat y}_n$ and $\mathbf{\hat X}_n$. Conditional on $\mathbf{\hat y}_n$ and $\mathbf{\hat{B}}$, the IF is nonstochastic; i.e., $\mathbf{V}_n = 0$ and $\mathcal{E}_{n} = \boldsymbol{\dot q}_{n}(\mathbf{y}_n,~\mathbf{\hat{B}}_n)$. We conduct a simulation study with linear IV models (see DGP A and B) and related technical details are provided in OA (ref).
Once we have $\mathbf{V}_n$ and $\mathcal{E}_n$, we can obtain $\mathbf{\hat V}_n$ and $\hat{\mathcal{E}}_{n,s}$. We recall that $\mathbf{\hat V}_n$ is an estimator of $\mathbf{V}_0 = \operatorname{plim} \mathbf{V}_n$, whereas $\hat{\mathcal{E}}_{n,s}$ are independent variables that share the asymptotic distribution of $\mathcal{E}_n$. Since $\mathcal{E}_{n}$ and $\mathbf V_n$ are conditional moments, given $\mathbf{\hat B}_{n}$, they are functions of $\mathbf{\hat B}_{n}$. They can also depend on $\boldsymbol{\theta}_0$ and $\mathbf{B}_0 = (\boldsymbol{\beta}_{0,1},~\dots,~\boldsymbol{\beta}_{0,n})^{\prime}$ because the specification of $\mathbf y_n$ depends on $\boldsymbol{\theta}_0$ and $\mathbf B_0$. In the rest of this section, we thus use the notations $\mathbf V_n(\mathbf{\hat B}_{n}, ~\boldsymbol{\theta}_0, ~ \mathbf{B}_{0}) \equiv \mathbf V_n$ and $\mathcal{E}_n(\mathbf{\hat B}_{n}, ~\boldsymbol{\theta}_0, ~ \mathbf{B}_{0}) \equiv \mathcal{E}_{n}$ to indicate that $\mathcal{E}_{n}$ and $\mathbf V_n$ are functions of $\mathbf{\hat B}_{n}$, $\boldsymbol{\theta}_0$, and $\mathbf{B}_0$.
The estimator of $\mathbf V_0$ can be obtained by replacing $\boldsymbol{\theta}_0$ and $\mathbf{B}_{0}$ in $\mathbf V_n(\mathbf{\hat B}_{n}, ~\boldsymbol{\theta}_0, ~ \mathbf{B}_{0})$ with consistent estimators; i.e., $\mathbf{\hat V}_n = \mathbf V_n(\mathbf{\hat B}_{n}, ~\boldsymbol{\hat \theta}_n, ~ \mathbf{\hat B}_{n})$. In contrast, we cannot obtain $\hat{\mathcal{E}}_{n,s}$ simply by replacing $\boldsymbol{\theta}_0$ and $\mathbf{B}_{0}$ in $\mathcal E_n(\mathbf{\hat B}_{n}, ~\boldsymbol{\theta}_0, ~ \mathbf{B}_{0})$ with consistent estimators. Indeed, what makes $\mathcal E_n(\mathbf{\hat B}_{n}, ~\boldsymbol{\theta}_0, ~ \mathbf{B}_{0})$ random is $\mathbf{\hat B}_{n}$. Consequently, an independent variable with the same asymptotic distribution as $\mathcal E_n(\mathbf{\hat B}_{n}, ~\boldsymbol{\theta}_0, ~ \mathbf{B}_{0})$ must be $\mathcal{E}_n(\mathbf{\bar B}_{n,s}, ~\boldsymbol{\theta}_0, ~ \mathbf{B}_{0})$, where $\mathbf{\bar B}_{n,s}$ and $\mathbf{\hat B}_{n}$ have the same asymptotic distribution. By replacing $\boldsymbol{\theta}_0$ and $\mathbf{B}_{0}$ in $\mathcal{E}_n(\mathbf{\bar B}_{n,s}, ~\boldsymbol{\theta}_0, ~ \mathbf{B}_{0})$ with their estimators, we obtain $\hat{\mathcal{E}}_{n,s} = \mathcal{E}_n(\mathbf{\bar B}_{n,s}, ~\boldsymbol{\hat \theta}_n, ~ \mathbf{\hat B}_{n})$. In practice, we simulate $\mathbf{\bar B}_{n,s}$ from the estimator of the distribution of $\mathbf{\hat B}_{n}$.
Our approach requires the practitioner to possess a consistent estimator of the joint distribution of $\boldsymbol{\hat \beta}_{n,1}$, \dots, $\boldsymbol{\hat \beta}_{n,n}$. This estimator is obtained in the first stage for a broad class of models. For estimators of type $\boldsymbol{\hat\beta}_{n,i} = f(\boldsymbol{z}_i, \boldsymbol{\hat\gamma}_n)$, which encompass semiparametric methods, a simulation from an estimator of the distribution of $\mathbf{\hat{B}}_n$ is $\mathbf{\bar{B}}_{n,s} = (\boldsymbol{\bar\beta}_{n,1}^{(s)},~\dots,~\boldsymbol{\bar\beta}_{n,n}^{(s)})^{\prime}$, where $\boldsymbol{\bar\beta}_{n,i}^{(s)} = f(\boldsymbol{z}_i, \boldsymbol{\bar\gamma}_{n}^{(s)})$ and $\boldsymbol{\bar\gamma}_{n}^{(s)}$ is simulated from an estimator of the asymptotic distribution of $\boldsymbol{\hat\gamma}_n$. From a frequentist perspective, the estimator of the distribution of $\boldsymbol{\hat \gamma}_n$ is typically derived through an asymptotic analysis (e.g., a normal distribution centered at $\boldsymbol{\hat\gamma}_n$ with some covariance matrix). In the Bayesian paradigm, we use the posterior distribution and obtain simulations by employing a Gibbs sampler or Metropolis-Hastings.
As pointed out above, $\mathbf V_n$ may not be tractable in some complex models when $\mathbf{\hat{B}}_n$ and the error term in the second stage are not independent. In these cases, we cannot construct the sample $\mathcal{S}_{\kappa}$. Instead, we apply Corollary (ref); i.e., we infer $\Delta_n$ only when its asymptotic distribution is normal. However, since we can generally compute (or approximate) $\mathcal E_n$, then we can estimate the asymptotic mean of $\Delta_n$, thereby debasing $\boldsymbol{\hat \theta}_n$ (see DGP D in our simulation study). We present our debiasing technique in the next section.
Plug-in estimators may exhibit significant bias when the first-stage sampling error is substantial. This issue can arise when many covariates are included in the first-stage estimation (e.g., IV approach with many instruments) or when the number of observations in the first stage is low compared to $n$. Although $\boldsymbol{\hat\theta}_n$ can still be consistent in these cases, the limiting distribution of $\Delta_n$ may not have a zero mean belloni2014inference, belloni2017program, cattaneo2019two. In this section, we discuss how our method can be applied to address this issue.
Theorem (ref) implies that the asymptotic mean of $\Delta_n$ is $\mathbb{E}(\Delta_0)= \mathbf{A}_0^{-1}\mathbb{E}(\mathcal{E}_0)$. Our approach accommodates scenarios where $\mathbb{E}(\mathcal{E}_0)$ is not zero. Condition ((ref)) of Assumption (ref) only states that $\mathcal{E}_n$ has a limiting distribution. The condition $\mathbb{E}(\mathcal{E}_0)\neq 0$ suggests that the plug-in estimator may exhibit significant finite sample bias. The good news is that we can estimate $\mathbf{A}_0$ and $\mathbb{E}(\mathcal{E}_0)$. Using these estimates, we propose a debiased estimator and establish its asymptotic distribution. We consider the estimator that is given by:
where $\boldsymbol{\hat\Omega}_{n}^{\kappa} = \frac{1}{\kappa}\sum_{s = 1}^{\kappa} \hat{\mathcal{E}}_{n,s}$ is an estimator of $\mathbb{E}(\mathcal{E}_0)$ as defined in Equation (ref). The following theorem establishes the consistency of $\boldsymbol{\theta}_{n,\kappa}^{\ast}$ and its limiting distribution.
As mentioned earlier, $\mathcal{E}_n(\mathbf{\bar B}_{n,s}, ~\boldsymbol{\theta}_0, ~ \mathbf{B}_{0})$ and $\mathcal{E}_n(\mathbf{\hat B}_{n}, ~\boldsymbol{\theta}_0, ~ \mathbf{B}_{0})$ have the same asymptotic distribution. By the LLN, $\frac{1}{\kappa}\sum_{s = 1} ^{\kappa}\mathcal{E}_n(\mathbf{\bar B}_{n,s}, ~\boldsymbol{\theta}_0, ~ \mathbf{B}_{0})$ converges in probability to $\mathbb{E}(\mathcal{E}_{0})$ as $n$ and $\kappa$ grow to infinity.\footnote{This requires that $\mathbb{E}(\lvert \mathcal{E}_n(\mathbf{\bar B}_{n,s}, ~\boldsymbol{\theta}_0, ~ \mathbf{B}_{0})\rvert^{\nu}) < \infty$ for some $\nu > 2$. Under this condition, the LLN implies that $\frac{1}{\kappa}\sum_{s = 1} ^{\kappa}\mathcal{E}_n(\mathbf{\bar B}_{n,s}, ~\boldsymbol{\theta}_0, ~ \mathbf{B}_{0})$ converges in probability to $\mathbb{E}(\mathcal{E}_n(\mathbf{\bar B}_{n,s}, ~\boldsymbol{\theta}_0, ~ \mathbf{B}_{0}))$ as $\kappa \to \infty$. Since $\mathcal{E}_n(\mathbf{\bar B}_{n,s}, ~\boldsymbol{\theta}_0, ~ \mathbf{B}_{0})$ converges in distribution to $\mathcal{E}_0$ (Assumption (ref), Condition ((ref))), it follows that $\mathbb{E}(\mathcal{E}_n(\mathbf{\bar B}_{n,s}, ~\boldsymbol{\theta}_0, ~ \mathbf{B}_{0}))$ converges to $\mathbb{E}(\mathcal{E}_{0})$ as $n\to\infty$. This holds because $\mathbb{E}(\lvert \mathcal{E}_n(\mathbf{\bar B}_{n,s}, ~\boldsymbol{\theta}_0, ~ \mathbf{B}_{0})\rvert^{\nu}) < \infty$ chung2001course.} Thus, by assuming that $\frac{1}{\kappa} \sum_{s = 1}^{\kappa} \hat{\mathcal{E}}_{n,s}$ converges in probability to $\mathbb{E}(\mathcal{E}_{0})$ as $n$ and $\kappa$ grow to infinity in Theorem (ref), we indirectly impose that $\mathcal{E}_n(\mathbf{\bar B}_{n,s}, ~\boldsymbol{\theta}_0, ~ \mathbf{B}_{0})$ is smooth in $\boldsymbol{\theta}_0$ and $\mathbf{B}_{0}$, ensuring that both $\frac{1}{\kappa}\sum_{s = 1} ^{\kappa}\mathcal{E}_n(\mathbf{\bar B}_{n,s}, ~\boldsymbol{\theta}_0, ~ \mathbf{B}_{0})$ and $\hat{\mathcal{E}}_{n,s}$ to have the same limit. This condition is similar to Assumption (ref) when $\boldsymbol{\theta}_0$ and $\mathbf{B}_{0}$ in $\frac{1}{n}\sum_{i = 1}^n \partial_{\boldsymbol{\theta}}\partial_{\boldsymbol{\theta}^\prime}q(\boldsymbol{\theta}_0, ~y_{i}, ~\boldsymbol{x}_{i}, ~\boldsymbol{\beta}_{0,i})$ are replaced with consistent estimators (see OA (ref)).
As in the case of the standard plug-in estimator, we can construct the sample $\mathcal{S}^{\ast}_{\kappa} =\{\boldsymbol{\hat \psi}_{n,s}^{\ast}: s = 1, ~\dots, ~ \kappa\}$ analogous to $\mathcal{S}_{\kappa}$ to estimate the CDF $F^{\ast}$. The variables $\boldsymbol{\hat \psi}_{n,s}^{\ast}$ are given by $\boldsymbol{\hat\psi}_{n,s}^{\ast} = (\mathbf{\hat A}_n^{\ast})^{-1}(\mathbf{\hat V}_{n}^{\ast})^{1/2}\boldsymbol{\zeta}_{s} + (\mathbf{\hat A}_n^{\ast})^{-1}(\hat{\mathcal{E}}_{n,s}^{\ast} -\boldsymbol{\hat \Omega}_{n}^{\ast\kappa})$, where $\boldsymbol{\hat \Omega}_{n}^{\ast\kappa} = \frac{1}{\kappa}\sum_{s = 1}^{\kappa} \hat{\mathcal{E}}_{n,s}^{\ast}$. The variables $\mathbf{\hat A}_n^{\ast}$, $\hat{\mathcal{E}}_{n,s}^{\ast}$, and $\mathbf{\hat V}_{n}^{\ast}$ are respectively defined as $\mathbf{\hat A}_n$, $\hat{\mathcal{E}}_{n,s}$, and $\mathbf{\hat V}_{n,}$ using $\boldsymbol{\theta}_{n,\kappa}^{\ast}$ and not $\boldsymbol{\hat\theta}_n$.
The asymptotic variance given in Statement ((ref)) of Theorem (ref) is the same as that of Theorem (ref). This implies that we can use the debiased estimator even though $\mathbf V_n$ cannot be estimated but the asymptotic distribution of $\sqrt{n}(\boldsymbol{\theta}_{n,\kappa}^{\ast} - \boldsymbol{\theta}_0)$ is normal (i.e., if $\mathcal{E}_0$ is normally distributed). In this case, it is not necessary to construct the sample $\mathcal{S}^{\ast}_{\kappa}$ since the asymptotic variance $\mathbf{A}_0^{-1}(\mathbf{V}_0 + \mathbb{V}(\mathcal{E}_0))\mathbf{A}_0^{-1}$ can be directly estimated, as is done in standard inference approaches newey1984method, newey1994asymptotic, ackerberg2012practical, chen2015sieve.
One issue regarding the debiased estimator is that the estimate of the finite sample bias of $\boldsymbol{\hat\theta}_{n}$ may itself be biased if the estimator of the FS distribution is biased. This situation can arise in small samples when the FS parameters are bounded, and their estimators exhibit large variances. In such cases, draws from the FS distribution may hit their bounds, leading to a biased estimate of $\mathbb{E}(\mathcal{E}_{0})$. We illustrate this issue in our simulation study through a Copula-GARCH model. To mitigate the bias in the estimate of $\mathbb{E}(\mathcal{E}_{0})$, we replace the sample mean $\boldsymbol{\hat\Omega}_{n}^{\kappa}$ in Equation (ref) with the sample median of $\{\hat{\mathcal{E}}_{n,s}, ~s= 1, ~\dots, ~\kappa\}$. The median correction yields better performance by avoiding the problem of outliers in $\hat{\mathcal{E}}_{n,s}$. Theorem (ref) remains valid under the median correction if the distribution of $\mathcal{E}_{0}$ is symmetric.
In this section, we conduct a Monte Carlo study to assess the finite sample performance of the proposed CDF estimator and the debiased estimator.
We study four DGPs, denoted as DGP A--D. The sample size $n$ takes values in $\{250, ~500, ~1000, \break 2000\}$. DGP A is a treatment effect model with endogeneity. The model is defined as follows: $$y_i = \theta_0d_i + \varepsilon_i, \quad d_i = \mathbbm{1}\{z_i > 0.5(\varepsilon_i + 1.2)\}, \quad z_i \sim \text{Uniform}[0, ~1], \quad \varepsilon_i \sim \text{Uniform}[-1, 1],$$
where $d_i$ is a treatment status indicator ($d_i = 1$ if $i$ is treated), $z_i$ is an instrument for the treatment, and $\theta_0 = 1$. The treatment is endogenous given that it is correlated with the error term $\varepsilon_i$. The practitioner observes an i.i.d sample of $(y_i, ~d_i, ~z_i)$. We estimate $\theta_0$ by the IV method. In the first stage, we predict $\mathbb{E}(y_i|z_i)$ using an OLS regression of $y_i$ on $z_i$. We also predict $\mathbb{E}(d_i|z_i)$ using an OLS regression of $d_i$ on $z_i$. In the second stage, we regress the prediction of $\mathbb{E}(y_i|z_i)$ on the prediction of $\mathbb{E}(d_i|z_i)$.
DGP B is similar to DGP A with the difference that many regressors are involved in the first stage. The practitioner has access to $k_n = O(\sqrt{n})$ instruments, $z_{1,i}, ~\dots, ~z_{k_n,i}$. We maintain the specification of $y_i$ from DGP A but change the treatment status as follows: $$d_i = \mathbbm{1}\{0.2 + z_{1,i} + z_{2,i} + z_{3,i} + z_{4,i} > 0.5(\varepsilon_i + 1.2)\}, \quad z_{1,i}, ~\dots, ~z_{k_n,i} \overset{i.i.d}{\sim} \text{Uniform}[0, ~0.2].$$ Only four instruments, $z_{1,i}, ~\dots, ~z_{4,i}$, are relevant for $d_i$, while the others are superfluous variables that are independent of $d_i$. Yet, we include all $k_n$ instruments in the FS regressions, creating a scenario with many weak instruments cattaneo2019two. We consider the cases where $k_n = \big\lfloor 2\sqrt{n} \big\rceil$ and $k_n = \big\lfloor 4\sqrt{n} \big\rceil$, where $\lfloor . \rceil$ is the rounding to the nearest integer.
DGP C is a Poisson model with a latent covariate that is defined as: $$y_i \sim \text{Poisson}(\exp(\theta_{0,1} + \theta_{0,2}p_i)), \quad p_i = \sin^2(\pi z_i), \quad z_i \sim \text{Uniform}[0, ~10], \quad d_i \sim \text{Bernoulli}(p_i),$$
where $p_i$ is an unobserved probability and $\boldsymbol{\theta}_0 = (\theta_{0,1}, ~\theta_{0,2})^{\prime} = (-0.8, ~2)^{\prime}$. We consider an i.i.d. sample comprising data points $(y_i, ~z_i, ~d_i)$. The practitioner observes the pairs $(y_i, ~z_i)$ for all $i$ but only observes $d_i$ for a representative subsample of size $n^{\ast} = \lfloor n^{\alpha_n} \rceil$. The parameter $\alpha_n$ takes values in $\{1, ~0.985, ~0.945, ~0.91\}$, in the same order as $n$, i.e., $\alpha_n = 1$ if $n = 250$, $\alpha_n = 0.985$ if $n = 500$, and so forth. In the FS, we estimate $p_i = \mathbb{E}(d_i|z_i)$ using a semiparametric regression of $d_i$ on $z_i$ in the subsample of size $n^{\ast}$ where $d_i$ is observed. We rely on a piecewise cubic spline regression hastie2017generalized. The regression results can be used to estimate $p_i$ for all $i$ in the full sample, as we observe $z_i$ for all $i$. The second stage is a standard Poisson regression after replacing $p_i$ with its estimator.
DGP D is a multivariate time-series model similar to the model used in the simulation study by gonccalves2022bootstrapping. We consider $k_n$ returns $y_{1,i}$, \dots, $y_{k_n, i}$, where $i$ is time and $k_n \geq 2$. Each $y_{p,i}$, for $p = 2, ~\dots,~ k_n$, follows an AR(1)-GARCH(1, 1) defined as: $$y_{p,i} = \phi_{p,0} + \phi_{p,1}y_{p,i-1}+ \sigma_{p,i}\varepsilon_{p,i}, \quad \sigma_{p,i}^2 = \beta_{p,0}+\beta_{p,1}\sigma_{p,i-1}^2\varepsilon_{p,i-1}^2+\beta_{p,2}\sigma_{p,i-1}^2,$$ where $\phi_{p,0} = 0$, $\phi_{p,i-1} = 0.4$, $\beta_{p,0} = 0.05$, $\beta_{p,1} = 0.05$, $\beta_{p,2} = 0.9$, and $\varepsilon_{p,i}$ follows a standardized Student distribution of degree-of-freedom $\nu_p = 6$. The number of returns, $k_n$, increases with the sample size and takes values in $\{2, ~3, ~5, ~8\}$, in the same order as $n$, i.e., $k_n = 2$ if $n = 250$, $k_n = 3$ if $n = 500$, and so forth. We account for the correlation between the returns using the Clayton copula nelsen2006introduction. The joint density function of $y_i = (y_{1,i}, ~\dots,~ y_{p,i})^{\prime}$ conditional on $\mathcal{F}^{i-1}$ (information set at $i-1$) is given by $c_i(G_{1,i}(\boldsymbol{\beta}_{0,1}), ~\dots,~ G_{k_n,i}(\boldsymbol{\beta}_{0,k_n}), ~\theta_0)$, where $\boldsymbol{\beta}_{0,p} = (\phi_{p,0}, ~\phi_{p,1}, ~\beta_{p,0}, ~\beta_{p,1}, ~\beta_{p,2}, ~\nu_p)^{\prime}$, $G_{p,i}(\boldsymbol{\beta}_{0,p})$ is the CDF of $y_{p,i}$ conditional on $\mathcal{F}^{i-1}$, and $c_i$ is the PDF of $k_n$-dimensional Clayton copula of parameter $\theta_0 = 4$. The practitioner observes the sample $y_1,~\dots,~y_n$. We rely on a multiple-stage estimation strategy to estimate $\theta_0$. In the first $k_n$ stages, we separately estimate each $\boldsymbol{\beta}_{0,p}$ by applying an AR(1)-GARCH(1, 1) model to the sample $y_{p,1}, ~\dots, y_{p,n}$. In the last stage, we estimate $\theta_0$ by maximum likelihood (ML) after replacing $\boldsymbol{\beta}_{0,p}$ in the density function of $y_i$ with its estimator.
This section presents the simulation results. We perform 10,000 simulations and set $\kappa$ to 1,000. We begin by examining the estimates of the CDF of $\Delta_n$ (see Figures (ref) and (ref)). We use the sample distribution of $\Delta_n$ resulting from the simulations as the benchmark for comparing our estimates. The corresponding CDF is represented by the curve $F_0$.
For DGPs A, B, and C, the conditional variance $\mathbf V_n$ is tractable. We thus use the simulation method described in Section (ref) to estimate the asymptotic CDF of $\Delta_n$. The average of the estimated CDFs (for the 10,000 simulations) is represented by the curve $\hat F_n$. The curve $\hat H_n$ displays the average CDF obtained using the assumption that the asymptotic distribution of $\Delta_n$ is normal with a zero mean. The variance of this normal distribution is estimated by our simulation method (Equation (ref)). In contrast, for DGP D, the conditional variance $\mathbf V_n$ is not tractable. We thus rely on Corollary (ref) and directly estimate the asymptotic variance of $\Delta_n$ using a standard inference method for sequential extremum estimators newey1994large. Both curves $\hat F_n$ and $\hat H_n$ are normal CDFs. The mean of the normal distribution for the curve $\hat F_n$ is estimated by our simulation method (see Section (ref)), while $\hat H_n$ corresponds to a normal CDF with a zero mean. For each estimated CDF, we present in parentheses the $L_1$-Wasserstein distance to the sample CDF $F_0$.\footnote{The $L_1$-Wasserstein distance between a CDF $\hat R_n$ and $F_0$ is $\lVert \hat R_n - F_0 \rVert_w = {\displaystyle\varint}_{\mathbb R}\lvert \hat R_n(t) - F_0(t) \rvert dt$.}
For DGP A, both the $\hat H_n$ and $\hat F_n$ approximations yield strong performance, with the $\hat H_n$ providing a better fit of the actual CDF according to the Wasserstein distance. The result is not surprising since the first- and second-stage estimators are finite-dimensional estimators of type M. Thus $\Delta_n$ is asymptotically normally distributed with a zero mean murphy2002estimation. Importantly, even when the asymptotic normality is verified, accurately computing the variance of a plug-in estimator can be intricate. Equation (ref), which is used for $\hat H_n$, accurately approximates this variance.
Including many superfluous variables in the first stage of an IV approach can lead to biased estimates cattaneo2019two. We can observe this result with DGP B given that the true distribution of $\Delta_n$ is not centered at zero. Yet, our inference method captures this bias, whereas $\hat H_n$, which corresponds to standard inference approaches, does not. The bias of our estimated CDF is larger when $k_n = \big\lfloor 4\sqrt{n} \big\rceil$, but vanishes as $n$ grows.\footnote{We construct 95% CIs for $\boldsymbol{\theta}_0$ using $\hat F_n$ and $\hat H_n$, and evaluate their coverage rates. The results, which are presented in OA (ref), show $\hat H_n$ leads to over-rejection of the null hypothesis that $\boldsymbol{\theta}_0$ is equal to its actual value.}
In the case of DGP C, the size of the FS sample does not grow at the same rate as $n$. This makes most classical inference approaches inapplicable. Our simulation approach performs well for both $\sqrt{n}(\hat{\theta}_{n,1} - \theta_{0,1})$ and $\sqrt{n}(\hat{\theta}_{n,2} - \theta_{0,2})$, where $\hat{\theta}_{n,1}$ and $\hat{\theta}_{n,2}$ represent the respective estimators for $\theta_{0,1}$ and $\theta_{0,2}$. Because of the slow convergence rate in the first stage, the CDFs of $\sqrt{n}(\hat{\theta}_{n,1} - \theta_{0,1})$ and $\sqrt{n}(\hat{\theta}_{n,2} - \theta_{0,2})$ are not centered at zero and our approach captures this feature. Conversely, the normal approximations perform poorly.
Simulation results for DGP D show that our approach performs well but is somewhat less accurate in small samples. This discrepancy arises because inference for GARCH models requires a large number of observations. Moreover, as the number of returns increases with $n$, $\sqrt{n}(\log(\hat{\theta}_n) - \log(\theta_0))$ exhibits bias. Our approach captures this bias, whereas the normal approximations perform poorly.
We now turn to the finite sample performance of the debiased estimator (see Table (ref)). For DGP A, where the classical plug-in estimator does not exhibit finite sample bias (because both stages result in a finite-dimensional extremum estimator), the debiased estimator closely aligns with the classical estimator. For DGP B, our debiased estimator significantly reduces the bias of the classical estimator. For example, with $2\sqrt{n}$ instrumental variables involved in the first stage, the bias of the classical plug-in estimator is substantial at $-0.101$ for $n = 250$. In contrast, the debiased estimator reduces this bias to $-0.017$. Not surprisingly, the bias reduction performs less well with $4\sqrt{n}$ instrumental variables but decreases the bias by more than 70%.
The results are similar for DGP C. The biases of the estimates for $\theta_{0,1}$ and $\theta_{0,2}$ are respectively $0.114$ and $-0.184$ when $n = 250$ and $0.033$ and $-0.053$ when $n = 2000$. In contrast, these biases are negligible for the debiased estimator.
The bias correction performs less well for DGP D. When $n = 2,000$, the debiased estimator's bias is ten times lower in absolute value (0.008 for the debiased estimator and $-$0.086 for the classical estimator). However, the correction is less effective in smaller samples. This result aligns with the discussion in Section (ref), which notes that the estimate of finite sample bias may be substantially biased when inequality constraints are imposed on the parameters in the first stage, especially if their estimates exhibit high variance. Constraints on GARCH model parameters include $\phi_{0,1}^2 < 1$ and $0< \beta_{p,1} + \beta_{p,3} < 1$. If the distribution from which we simulate in the first stage has high variance, the constraints may tighten for many draws, leading to an incorrect bias approximation. For instance, in the second stage, the standard deviation of the debiased estimator is 2.516, approximately 7 times higher than that of the classical estimator. In such situations, median correction is preferable. For $n = 250$, the debiased estimator using the median correction shows a lower bias of 0.088, which is 3 times smaller than the bias of the classical estimator.\footnote{In OA (ref), we present estimates of the asymptotic CDF of $\Delta_{n,\kappa}^{\ast}:=\sqrt{n}(\boldsymbol{\theta}_{n,\kappa}^{\ast} - \theta_0)$. Unlike the case of the classical plug-in estimator, the true sample CDFs are asymptotically centered at zero.}
In this section, we revisit the empirical analysis conducted by fortin2015peer on peer effects in adolescent fast-food consumption habits. Given the potential externalities that are associated with fast-food consumption and its link to overweight issues among adolescents, there may be justification for introducing a consumption tax on fast food. The optimal tax level hinges on the social multiplier of eating habits, emphasizing the need for an accurate measure of peer effects. We improve the popular IV estimator of peer effects by expanding the set of instruments, including many weak instruments. We reduce the bias of the estimate using our approach.
This section presents the model used in this application. We consider a set of $R$ schools, where the number of students in the $r$-th school is denoted by $n_r$. Students within the same school interact. The network in the $r$-th school is represented by an adjacency matrix $\mathbf{G}_r = [g_{r,ij}]_{\substack{i = 1, ~\dots,~ n_r \\j = 1, ~\dots,~ n_r}}$, where $g_{r,ij} = 1$ if student $j$ is a friend of student $i$ and $g_{r,ij} = 0$ otherwise. We restrict friendships to the same school; students from different schools cannot be friends. Moreover, self-friendships are not allowed, in the sense that $g_{r,ii} = 0$ for all $r$ and $i$. We consider the following linear-in-means peer effect model:
where $y_{r,i}$ is the weekly fast-food consumption frequency of student $i$ (reported frequency in days of fast-food restaurant visits in the past week), $\boldsymbol{x}_{r,i}$ is a vector of student $i$'s observable characteristics, $\varepsilon_{r,i}$ is an error term assumed to be independent of $\boldsymbol{x}_{r,i}$ and $\mathbf{G}_r$, and $n_{r,i} = \sum_{j = 1}^{n_r}g_{r,ij}$ is the number of friends of student $i$. The parameter $\theta_{0,1}$ captures peer effects, which measure the influence of an increase in the average friend's fast-food consumption frequency on one's fast-food consumption frequency.\footnote{The uniqueness of equilibrium in this model requires $\lvert \theta_{0,1} \rvert < 1$. That is, students do not increase their consumption frequency greater than the increase in their average friends’ consumption frequency bramoulle2009identification.} The parameter $\boldsymbol{\theta}_{0,2}$ reflects the effect of student characteristics, whereas $\boldsymbol{\theta}_{0,3}$ captures contextual effects, i.e., the influence of the average observable characteristics among friends. The parameter $\alpha_{0,r}$ accounts for unobserved effects of school characteristics, such as school location, regional taxes, and pricing policies.
The average friend's fast-food consumption frequency, which is measured by $\bar y_{r,i} = \frac{\sum_{j = 1}^{n_r}g_{r,ij} y_{r,j}}{n_{r,i}}$, is endogenous. Consequently, the classical OLS estimates of the parameters in (ref) are likely to be inconsistent. Fortunately, one does not need to seek instruments for $\bar y_{r,i}$ elsewhere; they can be generated from the model. For the sake of clarity, we rewrite Equation (ref) in a matrix form for school $r$. Let $\mathbf{y}_r = (y_{r,1}, ~\dots, ~y_{r,n_s})^{\prime}$, $\mathbf{X}_r = (\boldsymbol{x}_{r, 1}, ~\dots, ~ \boldsymbol{x}_{r, n_r})^{\prime}$, and $\boldsymbol{\varepsilon}_r = (\varepsilon_{r,1}, ~\dots, ~\varepsilon_{r,n_s})^{\prime}$.\footnote{The notations $\mathbf{y}_r$ and $\mathbf{X}_r$ are only used in this section and must not be confused with $\mathbf{y}_n$ and $\mathbf{X}_n$ used elsewhere.} Let also the row-normalized adjacency matrix $\mathbf{\tilde G}_r = [\tilde g_{r,ij}]_{\substack{i = 1, ~\dots,~ n_r \\j = 1, ~\dots,~ n_r}}$, where $\tilde g_{r,ij} = 1/\sum_{j = 1}^{n_r} g_{r,ij}$ if $j$ is an $i$'s friend and $\tilde g_{r,ij} = 0$ otherwise. The linear-in-means peer effect model at the school level is:
where $\mathbf{1}_{n_r}$ is an $n_r$-dimensional vector of ones. By premultiplying the terms of Equation (ref) by $\mathbf{\tilde G}_r$, we can observe that $\mathbf{\tilde G}_r^2 \mathbf{X}_r$ is correlated with the endogenous variable $\mathbf{\tilde G}_r\mathbf{y}_r$ if $\boldsymbol{\theta}_{0,3}\ne 0$. Since $\mathbf{\tilde G}_r^2 \mathbf{X}_r$ is not an explanatory variable in Equation (ref), it can thus serve as an excluded instrument kelejian1998generalized, bramoulle2009identification. This instrument is interpreted as the average friends of their average friends of the characteristics $\boldsymbol{x}_{r,i}$.
However, the instrument $\mathbf{\tilde G}_r^2 \mathbf{X}_r$ might suffer from weakness if $\boldsymbol{\theta}_{0,3} \approx 0$. To address this concern, we can show from Equation (ref) that:
Therefore, $\mathbb{E}(\mathbf{\tilde G}_r\mathbf{y}_r|\mathbf{X}_r, \mathbf{\tilde G}_r)$ can be used as an instrument. This approach is optimal since $\mathbb{E}(\mathbf{\tilde G}_r\mathbf{y}_r|\mathbf{X}_r, \mathbf{\tilde G}_r)$ fully captures exogenous variations of the endogenous variable $\mathbf{\tilde G}_r\mathbf{y}_r$. In practice, employing this instrument entails a TS IV approach. A first IV method with $\mathbf{\tilde G}_r^2 \mathbf{X}_r$ as an instrument is used to estimate $\boldsymbol{\theta}_0 = (\alpha_{0,r}, ~\theta_{0,1}, ~\boldsymbol{\theta}_{0,2}^{\prime}, ~\boldsymbol{\theta}_{0,3}^{\prime})^{\prime}$. This estimate can also be employed to approximate $\mathbb{E}(\mathbf{\tilde G}_r\mathbf{y}_r|\mathbf{X}_r, \mathbf{\tilde G}_r)$ by replacing $\boldsymbol{\theta}_0$ in Equation (ref) with its estimate. A second IV method is performed with the estimate of $\mathbb{E}(\mathbf{\tilde G}_r\mathbf{y}_r|\mathbf{X}_r, \mathbf{\tilde G}_r)$ as an instrument to estimate $\boldsymbol{\theta}_0$.
While the optimal IV approach has been thoroughly considered in the literature, it may not entirely resolve the issue of weak instruments. Specifically, if the instrument $\mathbf{\tilde G}_r^2 \mathbf{X}_r$ that is used in the first IV approach is weak, the estimation of $\boldsymbol{\theta}_0$ can be biased, leading to a biased estimate for $\mathbb{E}(\mathbf{\tilde G}_r\mathbf{y}_r|\mathbf{X}_r, \mathbf{\tilde G}_r)$. To circumvent this problem, we propose a new approach that consists of expanding the set of instruments for $\mathbf{\tilde G}_r\mathbf{y}_r$. As $(\mathbf{I}_{n_r} - \theta_{0,1}\mathbf{\tilde G}_r)^{-1} = \sum_{p = 0}^{\infty}\theta_{0,1}^p\mathbf{\tilde G}_r^p$, it can be shown from Equation (ref) that: $$\textstyle\mathbb{E}(\mathbf{\tilde G}_r\mathbf{y}_r|\mathbf{X}_r, \mathbf{\tilde G}_r) = \alpha_{0,r}\mathbf{\tilde G}_r(\mathbf{I}_{n_r} - \theta_{0,1}\mathbf{\tilde G}_r)^{-1}\mathbf{1}_{n_r} + \mathbf{\tilde G}_r\mathbf{X}_r\boldsymbol{\theta}_{0,2} + \sum_{p = 0}^{\infty}\theta_{0,1}^p\mathbf{\tilde G}_r^{2+p}\mathbf{X}_r(\theta_{0,1}\boldsymbol{\theta}_{0,2} + \boldsymbol{\theta}_{0,3}).$$ This suggests the use of $\mathbf{\tilde G}_r^{2+p}\mathbf{X}_r$, for $p=0,~1,~2, ~\dots, k_{\max}$ as instruments, where $k_{\max}$ can be as large as possible for the matrix of instruments to be full rank. This set of instruments can be interpreted as averages of $\boldsymbol{x}_{r,i}$ among close- and long-distance friends. We control for 25 student characteristics in $\boldsymbol{x}_{r,i}$ (see below) and set $k_{\max} = 9$. This leads to 250 excluded instruments for $\mathbf{\tilde G}_r\mathbf{y}_r$. Despite this large number of instruments, our inference method can be used to construct CIs for the parameters of the model. Moreover, following Theorem (ref), we correct for the finite sample bias of the resulting IV estimator.
We use data from the National Longitudinal Study of Adolescent to Adult Health (Add Health) survey. The purpose of this survey was to investigate how various social contexts (families, friends, peers, schools, neighborhoods, and communities) affect adolescents' health and risk behaviors. The survey provides nationally representative and detailed information on adolescents in grades 7--12 from 144 schools during the 1994--1995 school year in the United States (US). Approximately 90,000 students were asked to complete a brief questionnaire covering demographics, family backgrounds, academic performance, and health-related behaviors, as well as friendship links (best friends within the same school, up to 5 females and up to 5 males). From this main sample, an in-home sample (core sample) of about 20,000 students was randomly selected. These students participated in a more extensive questionnaire featuring detailed questions. This subsample has been followed in subsequent waves of the survey.
We use the Wave II dataset, which encompasses most of the variables that are relevant to this study. This wave targets the subsample of 20,000 students tracked over time. However, only 16 schools, comprising approximately 3,000 students, were completely surveyed in Wave II. To avoid the issue of sampled networks, we focus on the sample that was derived from these 16 schools, where we can observe the entire list of nominated best friends, and all nominated friends are also surveyed. In addition to the Wave II dataset, we gather information on each student's race and their mother's background from Wave I.
Our dependent variable is the weekly fast-food consumption frequency, measured by the reported frequency (in days) of fast-food restaurant visits in the past week. The final sample consists of 2,735 students, with an equal distribution between boys and girls. We control for 25 observable characteristics in $\mathbf{X}_r$ such as students' gender, grade, race, weekly allowance, and parents' education and occupation. On average, students report consuming fast food 2.35 days per week. The average age of students at the time of Wave II data collection is 16.62 years. Additional details on the data summary can be found in Table (ref) in OA (ref).
Figure (ref) displays estimates of peer effects using our approach and alternative methods, including the OLS estimator, the classical IV (CIV) estimator, the optimal IV (OIV) estimator, our IV estimator with many instruments (IV-MI), and the corresponding debiased IV estimator with many instruments (DIV-MI).\footnote{See full results, including coefficients of control variables in OA (ref).} The OLS approach overlooks the endogeneity issue, whereas the CIV method uses $\mathbf{\tilde G}_r^2 \mathbf{X}_r$ as an instrument. The OIV estimator employs the estimate of $\mathbb{E}(\mathbf{\tilde G}_r\mathbf{y}_r|\mathbf{X}_r, \mathbf{\tilde G}_r)$ as an instrument, replacing unknown parameters in Equation (ref) with their CIV estimates.
The OLS estimate indicates that the peer effect parameter is significant. The estimate decreases from 0.192 to 0.150 when we control for school-fixed effects. In contrast, the CIV estimator has a large variance, indicating that the coefficient is not statistically significant. This imprecision is a consequence of the weakness of the instruments mikusheva2022inference, leading to a biased estimator toward the OLS one. For instance, although it is known that the model suffers from an endogeneity problem, the Hausman-Wu endogeneity test (not reported here) indicates that the OLS and CIV estimators are not significantly different.
As discussed earlier, this issue also invalidates the OIV approach since biased CIV estimates are used to estimate the optimal instrument. Notably, we observe that the 95% confidence interval of the OIV is even larger than that of CIV. While the estimator becomes more precise when we control for school-fixed effects, the results still indicate that peer effects are not significant.
These findings align with the results of fortin2015peer. Their OIV estimate of the peer effect parameter is 0.110 with a standard error approximated at 0.395, indicating non-significance. Additionally, they implemented a quasi-maximum likelihood (QML) method, estimating peer effects at 0.129, but the coefficient is significant only at the 10% level.
After expanding the pool of instruments, the IV estimator estimate with many instruments reveals significant peer effects. The estimate decreases from 0.276 to 0.208 after accounting for school-fixed effects. Moreover, the results highlight evidence of finite sample bias, as the confidence intervals are not centered on the estimates. This bias is a consequence of the numerous instrumental variables in the first stage. After correcting for this bias, the estimates slightly increase to 0.300 and 0.218, respectively. The increase is smaller when controlling for school-fixed effects.
In summary, our results highlight the presence of peer effects in adolescent fast-food consumption habits. These results have two important implications. First, key players in the network can play a crucial role as channels for influencing adolescent habits. The more key players are influenced, the greater the potential spread of this influence in the network ballester2006s, zenou2016key. Second, the social multiplier becomes crucial in determining the impact of a policy on adolescent fast-food consumption frequency. The social multiplier coefficient, given by $1/(1 - \theta_{0,1})$, is estimated at 1.279 for the model with fixed effects, considering the finite sample bias. This implies that the effect of a tax increase on fast-food consumption frequency, in an environment where adolescents do not interact with each other, must be multiplied by 1.279 when they interact agarwal2021thy.
This paper proposes a simulation-based approach to estimate the asymptotic CDF of two-stage estimators. We consider a broad class of estimators in the first stage and an extremum estimator in the second stage.
A key challenge in inference using two-stage estimation methods is that the asymptotic distribution of the plug-in estimator is influenced by first-stage sampling error. We tackle this issue by disentangling the first-stage sampling error from the model error in the second stage. Notably, we consider the possibility that the limiting distribution of $\sqrt{n}(\boldsymbol{\hat\theta}_n - \boldsymbol{\theta}_0)$ may not be normal and not centered at zero. This flexibility enables bias reduction for the plug-in estimator, particularly in cases where first-stage sampling error induces significant bias in the second stage. We introduce a debiased plug-in estimator and demonstrate that its limiting distribution has a zero mean. We conduct an extensive simulation study confirming the finite sample performance of our debiased estimator.
We leverage the proposed approach to study peer effects on adolescents' fast-food consumption habits. We employ an instrumental variable (IV) approach and address the issue of biased estimates that can arise when instruments are weak. By incorporating the characteristics of both close and distant friends, we expand the IV set to include a larger number of weak instruments. We correct the finite sample bias of the IV estimator and provide valid inference.
This paper contributes to the emerging literature on inference methods when standard regularity conditions are violated. The proposed approach is straightforward to implement and does not impose restrictions on the class of first-stage estimators. Our method is also suitable for complex models. Unlike resampling methods, it avoids the need for multiple computations of the estimators, except when resampling is required to obtain the asymptotic distribution in the first stage.