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.
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.
Inference for Fixed Effects Estimators when Panels are Unbalanced
\thispagestyle{empty}
abstractWe derive the asymptotic properties of two-way fixed effects M-estimators with missing observations in an asymptotic framework in which the numbers of cross-sectional units and time periods grow jointly. We allow the selection process to be deterministic (conditional on the unobserved effects and initial conditions), stochastic, or mixed, and we impose only a conditional mean restriction. The uncorrected estimators are asymptotically normal but not centered at zero, suffering from incidental parameter and feedback biases. Feedback bias can be induced by predetermined regressors in the outcome equation and by a predetermined selection process. We propose debiased estimators that handle both sources without requiring knowledge of which regressors or selection components are predetermined.\\[1em]
JEL Classification: C13, C23\\
Keywords: panel data, unbalanced panel data, dynamic model, two-way fixed effects, incidental parameter problem, asymptotic bias correction.
\onehalfspacing
\setcounter{page}{1}
Introduction
Although unbalanced panel data are common in empirical work, the theoretical econometrics literature focuses mostly on balanced panels. Missing observations may arise because units enter and exit the sample depending on their unobserved heterogeneity, or because nonresponse depends on past outcomes. Both mechanisms matter for inference.
We provide the first complete theoretical analysis of two-way fixed effects M-estimators in a large-$N,T$ asymptotic framework that accommodates missing observations due to sample selection and attrition. Throughout the paper, we refer to the process that generates the missing observations as the selection process. Whereas the literature typically considers the selection process to be either deterministic (conditional on the unobserved effects and initial conditions) or stochastic (see, e.g., Chapter 17 of h2022), we introduce a mixed process that encompasses both pure cases. For the stochastic component, we build on Chapter 19 of w2010 and impose a general conditional mean restriction rather than a full conditional independence assumption. We derive the asymptotic distribution of the fixed effects estimators and characterize their bias and variance. In contrast to balanced panels, both bias and variance depend on the selection process. The uncorrected estimators are asymptotically normal but not centered at zero. We derive debiased estimators whose asymptotic distributions are centered at zero.
Our contributions also have important implications for applied work. Predetermined regressors introduce feedback bias that must be corrected for valid inference. This holds not only for predetermined regressors in the outcome equation, such as a lagged outcome, but also for those in the selection process. Feedback bias can therefore arise even when the outcome equation depends solely on strictly exogenous regressors. A key advantage of our bias correction is that it requires no knowledge of which regressors are predetermined, whether in the outcome equation or in the selection process. Furthermore, the selection process introduces an additional layer of heterogeneity, leading to two consequences. First, when the design depends on unobserved effects, the heterogeneity carries over to the subsamples on which jackknife corrections are based, invalidating them, as our simulations confirm. Our analytical correction does not require homogeneity. Second, the additional heterogeneity can generate heteroskedasticity, which conventional heteroskedasticity-robust standard errors can address.
\noindentRelated literature. Our paper is related to the large-$N,T$ literature on fixed effects estimators. This strand uses bias correction methods to address the inference problems caused by incidental parameters, first identified by ns1948 and n1981. Prominent debiasing methods include analytical and jackknife corrections (see, e.g.,\ hn2004), split-panel jackknife corrections (see, e.g.,\ dj2015), and bootstrap corrections (see, e.g.,\ ks2016). Most contributions focus on one-way fixed effects estimators. For panel models with two-way unobserved effects, as we consider here, the literature is sparser. fw2016 develop the first asymptotic theory of fixed effects estimators for such models and derive analytical and split-panel jackknife bias corrections. Other work has focused on more specialized settings, such as logit or network models, or strictly exogenous regressors (see g2017, yjfl2018, jo2019, d2019, h2026). We build on the analysis of fw2016 because it allows for a broad class of panel and network models with both strictly exogenous and predetermined regressors.
The problem of missing observations has received little attention in the large-$N,T$ literature on fixed effects estimators. fw2018 conjecture the asymptotic distributions and bias expressions of fixed effects estimators under a conditional independence assumption regarding the selection process. We relax this assumption to a conditional mean restriction, allowing for a more general selection process. We also show that one of their conjectured regularity assumptions is insufficient and must be strengthened to ensure the invertibility of the expected incidental parameter Hessian. Overall, we confirm their conjecture that bias-corrected estimators can be formed in a straightforward manner even when observations are missing. For one-way fixed effects estimators, dj2015 briefly show how the profile-likelihood correction can be adapted. Two further contributions on missing observations in large-$N,T$ settings are sww2026 and cs2026interactive. Both focus on linear models with interactive fixed effects, while we consider a broader class of models that includes nonlinear models with additive fixed effects.
\noindentOutline. Section (ref) introduces the model and the fixed effects estimators. Section (ref) presents the asymptotic theory. Section (ref) reports the simulation results. Section (ref) concludes.
\noindentNotation. Throughout the paper, $\mathbb{P}\left(\cdot\right)$ and $\mathbb{E}\left[ \cdot \right]$ denote probability and expectation. A superscript “0” on a parameter denotes its true population value. We write a.\,s.\ and wpa1 for almost surely and with probability approaching one, respectively. We use $o_{P}(1)$ to denote a sequence of random variables that converges in probability to zero and $\mathcal{O}_{P}(1)$ to denote a sequence that is bounded in probability. For two positive sequences $a_{NT}$ and $b_{NT}$, we denote by $a_{NT} \asymp b_{NT}$ that the two sequences are of the same order of magnitude. We use $\xrightarrow{d}$ and $\xrightarrow{p}$ to denote convergence in distribution and in probability, respectively. Unless otherwise stated, all stochastic statements are understood as almost-sure statements conditional on $\Phi$, the sigma-algebra generated by the unobserved effects and initial conditions.
Model and Estimator
We observe panel data $\{(y_{it}, x_{it}, s_{it}) \colon i \in \{1, \ldots, N\}, \, t \in \{1, \ldots, T\}\}$, where $y_{it}$ is an outcome variable, $x_{it}$ is a $K$-dimensional vector of predetermined regressors, and $s_{it}$ is a selection indicator that equals $s_{it} = 1$ if $y_{it}$ and every component of $x_{it}$ are observed, and $s_{it} = 0$ otherwise. The panel is unbalanced whenever $s_{it} = 0$ for some $(i, t)$.
We factor the selection indicator as $s_{it} = d_{it} r_{it}$, where $d_{it} \in \{0, 1\}$ is a $\Phi$-measurable design indicator and $r_{it} \in \{0, 1\}$ is a response indicator.
definition[Selection Process]
The selection process is deterministic if $r_{it} = 1$ for all $i, t, N, T$, stochastic if $d_{it} = 1$ for all $i, t, N, T$, and mixed otherwise.
Deterministic selection arises from the sampling design when, for example, entry and exit dates depend on unobserved effects. Stochastic selection arises from nonresponse or attrition, which may depend on past outcomes. The mixed case combines the two. Our simulation experiments in Section (ref) study such a mixed design.
We consider the following semiparametric unobserved effects model for $i \in \{1, \ldots, N\}$ and $t \in \{1, \ldots, T\}$,
equation[equation omitted — 152 chars of source]
where $\mathcal{C}_{i}^{t} \coloneqq \sigma(\{(x_{it^{\prime}}^{\prime}, r_{it^{\prime}}) \colon t^{\prime} \in \{1, \ldots, t\}\})$, $g(\pi)$ is a known link function, and $\pi_{it}(\beta, \mu) \coloneqq x_{it}^{\prime} \beta + \mu$ is the linear index. Conditioning on $\Phi$ is implicit throughout, so $d_{it}$ does not appear in $\mathcal{C}_{i}^{t}$. The regressors and the response indicator are predetermined. Hence, both may be functions of lagged outcomes but not of contemporaneous ones. Here, $\beta$ is a $K$-dimensional vector of model parameters, $\phi$ is an $L$-dimensional vector of nuisance parameters (with $L \coloneqq N + T$), and $\mu$ is a scalar. Let $\phi \coloneqq (\alpha^{\prime}, \gamma^{\prime})^{\prime}$, so that $\mu_{it}(\phi) \coloneqq \alpha_{i} + \gamma_{t}$. In economic applications, $\alpha_{i}$ and $\gamma_{t}$ are referred to as unobserved individual and time effects, which capture individual heterogeneity and common shocks, respectively. We impose no restrictions on the relationship between the regressors and the unobserved effects, and we do not assume a parametric distribution for the latter.
We follow a fixed effects approach and estimate $\phi$ jointly with $\beta$. We take $(\beta^{0}, \phi^{0})$ to be the unique minimizer of the population problem,
equation[equation omitted — 265 chars of source]
where the objective function,
equation[equation omitted — 258 chars of source]
consists of a criterion function $\psi_{it}(\pi) \coloneqq \psi(y_{it}, \pi)$ and a penalty term $\mathcal{P}_{NT}(\phi)$. The penalty imposes the normalization $\sum_{i = 1}^{N} \alpha_{i} = \sum_{t = 1}^{T} \gamma_{t}$, which removes the invariance of $\mu_{it}(\phi)$ to the shift $(\alpha_{i}, \gamma_{t}) \mapsto (\alpha_{i} + c, \gamma_{t} - c)$ for $c \in \operatorname{\mathbb{R}}$. We write $\mu_{it}^{0} \coloneqq \mu_{it}(\phi^{0}) = \alpha_{i}^{0} + \gamma_{t}^{0}$.
We tie the criterion to the model by requiring the (pseudo-)score of $\psi_{it}(\pi)$ to be affine in $y_{it}$,
equation[equation omitted — 125 chars of source]
where $w(\pi)$ is a known weight. The same link $g(\pi)$ therefore enters both the conditional mean in (ref) and the criterion function in (ref). Integrating gives $\psi_{it}(\pi) = \int^{\pi} w(z)\,(g(z) - y_{it}) \, \mathrm{d} z$, so $g(\pi)$ and $w(\pi)$ determine the criterion function up to a $\pi$-independent term. For example, $g(\pi) = \pi$ and $w(\pi) = 1$ yield $\psi_{it}(\pi) = (y_{it} - \pi)^{2}/2$, i.e., OLS. The affine form is a natural choice when only the conditional mean is restricted. We focus on objective functions that are strictly convex in a neighborhood of the true parameter values and sufficiently smooth in all parameters. This encompasses a wide range of estimators, including many popular (pseudo-)maximum-likelihood and (non)linear least-squares estimators, but it excludes non-smooth criterion functions, such as the check function of quantile regression.
We estimate $\beta^{0}$ and $\phi^{0}$ by minimizing the sample analog of (ref). The M-estimator is
equation[equation omitted — 251 chars of source]
Asymptotic Theory
Assumptions
Let $X$ denote the $NT \times K$ matrix of regressors, where $x_{it}^{\prime}$ is its $it$-th row, and let $\mathcal{T}_{i} \coloneqq \{t \in \{1, \ldots, T\} \colon s_{it} = 1\}$ denote the set of periods in which unit $i$ is observed. Let $(\mathrm{d}^{\mathrm{r}} \psi(\beta, \phi))_{it} \coloneqq \partial^{\mathrm{r}} \psi_{it}(\pi_{it}(\beta, \mu_{it}(\phi))) / \left(\partial \pi_{it}\right)^{\mathrm{r}}$ denote the $\mathrm{r}$-th derivative of the criterion function with respect to the linear index, and let $(\mathrm{d}^{\mathrm{r}}_{\mathcal{C}} \psi)_{it} \coloneqq \mathbb{E}\big[(\mathrm{d}^{\mathrm{r}} \psi)_{it} \mid \mathcal{C}_{i}^{t}\big]$ denote the corresponding conditional expectation. When the derivatives are evaluated at the true parameter values, we suppress their arguments.
We also define population weighted least-squares projections. For each $k \in \{1, \ldots, K\}$, let
align[align omitted — 522 chars of source]
where $x_{it, k}$ denotes the $k$-th element of $x_{it}$. The resulting $NT \times K$ matrix of fitted values is $\mathfrak{X}$, with $it$-th row $\mathfrak{x}_{it} \coloneqq (\mu_{it}(\xi_{1}^{0}), \ldots, \mu_{it}(\xi_{K}^{0}))^{\prime}$. The fitted values are therefore the weighted least-squares projection of $\mathbb{E}\left[ r_{it} \, x_{it} (\mathrm{d}^{2}_{\mathcal{C}} \psi)_{it} \right] / \mathbb{E}\left[ r_{it} (\mathrm{d}^{2}_{\mathcal{C}} \psi)_{it} \right]$ onto the span of the additive effects, with weights $\mathbb{E}\left[ s_{it} (\mathrm{d}^{2}_{\mathcal{C}} \psi)_{it} \right]$. We define the population residuals as $\ddot{X} \coloneqq X - \mathfrak{X}$, whose $it$-th row is $\ddot{x}_{it}^{\prime}$.
We impose the following assumption.
assumption[Sampling and Regularity Conditions]
Let $z_{it} \coloneqq (y_{it}, x_{it}^{\prime}, r_{it})^{\prime}$, $\nu \in (0, 1)$, $\delta \coloneqq 2 \kappa + \nu$ for some integer $\kappa \geq 3$, and $\varphi > \delta (\delta - \nu)/(2 \nu)$. Furthermore, let $\varepsilon > 0$ and let $\Theta^{0}(\varepsilon)$ be a subset of $\operatorname{\mathbb{R}}^{K + 1}$ that contains an $\varepsilon$-neighborhood of $(\beta^{0}, \mu_{it}^{0})$ for all $i, t, N, T$.
\begin{enumerate}[(i)]
• Asymptotics: We consider joint limits of sequences in which both panel dimensions diverge proportionally: $N, T \rightarrow \infty$ with $N / T \rightarrow \tau \in (0, \infty)$.
• Sampling: Conditional on $\Phi$, $\{\{z_{it}\}_{t = 1}^{T} \colon i \in \{1, \ldots, N\}\}$ is independent across $i$, and, for each $i$, $\{z_{it}\}_{t = 1}^{T}$ is $\alpha$-mixing with mixing coefficients satisfying $\sup_{i} a_{i}(m) = \mathcal{O}(m^{- \varphi})$ a.\,s.\ as $m \rightarrow \infty$, where $\mathcal{A}_{i}^{t}$ is the sigma-algebra generated by $(z_{it}, z_{i(t - 1)}, \ldots)$, $\mathcal{B}_{i}^{t}$ is the sigma-algebra generated by $(z_{it}, z_{i(t + 1)}, \ldots)$, and
\begin{equation*}
a_{i}(m) \coloneqq \sup_{t} \sup_{A \in \mathcal{A}_{i}^{t}, B \in \mathcal{B}_{i}^{t + m}}
\left\lvert \mathbb{P}\left(A \cap B\right) - \mathbb{P}\left(A\right) \mathbb{P}\left(B\right) \right\rvert \quad a.\,s.
\end{equation*}
• Model: We assume that for all $i, t, N, T$, $\mathbb{E}\left[ y_{it} \mid \mathcal{C}_{i}^{t} \right] = g(\pi_{it}(\beta^{0}, \mu_{it}^{0}))$. The unobserved effects $\phi^{0}$ are normalized to $\sum_{i = 1}^{N} \alpha_{i}^{0} = \sum_{t = 1}^{T} \gamma_{t}^{0}$.
• Smoothness and Moments: We assume that $(\beta, \mu) \rightarrow \psi_{it}(\pi_{it}(\beta, \mu))$ is four times continuously differentiable over $\Theta^{0}(\varepsilon)$ a.\,s. Each element of $x_{it}$, and the partial derivatives of $\psi_{it}(\pi_{it}(\beta, \mu))$ with respect to the elements of $(\beta, \mu)$ up to fourth order, are bounded in absolute value uniformly over $(\beta, \mu) \in \Theta^{0}(\varepsilon)$ by a function $\Psi((y_{it}, x_{it}^{\prime})) > 0$ a.\,s. In addition, $\sup_{it} \mathbb{E}\left[ (\Psi((y_{it}, x_{it}^{\prime})))^{\delta} \right]$ is a.\,s.\ uniformly bounded over $N, T$.
• Convexity: There exists a constant $c_{H} > 0$ such that $\inf_{\{(i, t) \colon d_{it} = 1\}} \mathbb{E}\left[ r_{it} (\mathrm{d}^{2}_{\mathcal{C}} \psi)_{it} \right] \geq c_{H}$ a.\,s.\ uniformly over $i, t, N, T$. Furthermore, there exists a constant $c_{W} > 0$ such that
\begin{equation*}
\underset{\{w \in \mathbb{R}^{K} \colon \lVert w \rVert_{2} = 1\}}{\min} \; \frac{1}{NT} \sum_{i = 1}^{N} \sum_{t = 1}^{T} \mathbb{E}\left[ s_{it} (\mathrm{d}^{2}_{\mathcal{C}} \psi)_{it} (\ddot{x}_{it}^{\prime} w)^{2} \right] \geq c_{W} \quad \text{a.\,s.}
\end{equation*}
• \textit{Design}: There exist constants $c_{U}, c_{O} \in (0, 1]$ such that $\inf_{i} T^{- 1} \sum_{t = 1}^{T} d_{it} \geq c_{U}$ and $\inf_{t \neq t^{\prime}} N^{- 1} \sum_{i = 1}^{N} d_{it} d_{it^{\prime}} \geq c_{O}$ a.\,s.\ uniformly over $N, T$.
\end{enumerate}
remark[Assumption (ref)]
We comment on each part of Assumption (ref) in turn.
\begin{enumerate}[(i)]
• The asymptotic framework, in which the cross-sectional and time dimensions diverge at proportional rates, is the same as in hk2011 and fw2016.
• We restrict dependence in the sample $\{(y_{it}, x_{it}, r_{it})\}$ in two ways: independence across units conditional on $\Phi$, and $\alpha$-mixing (strong mixing) within units, with coefficients that decay at a polynomial rate. Importantly, we do not require $r_{it}$ to be independent of $(y_{it}, x_{it})$ conditional on $\Phi$. That is, we do not assume missing-at-random, which this literature commonly invokes to treat the missingness as ignorable. We follow hk2011 and fw2016 and use $\alpha$-mixing because it is the weakest among the standard mixing conditions in econometrics and is preserved under measurable transformations. Any other form of weak dependence under which the regularity conditions underlying our results remain satisfied could replace it.
• We impose a conditional mean restriction to identify the model parameters. The regressors and the response indicator may be predetermined. Both may be arbitrarily related to past outcomes, but not to contemporaneous or future ones. Together with the affine score (ref), the conditional mean restriction implies that the score is conditionally mean-zero, $\mathbb{E}\left[ (\mathrm{d}^{1} \psi)_{it} \mid \mathcal{C}_{i}^{t} \right] = 0$ for all $i, t, N, T$. We use this repeatedly below. Beyond the first moment, the conditional distribution of $y_{it}$ is left unrestricted. Imposing further restrictions on higher-order conditional moments, as in fw2016, would allow one to exploit Bartlett identities (see b1953) to simplify the bias and variance expressions. This can improve the finite-sample performance of the debiased estimator, as noted by f2009. The normalization ensures unique identification.
• We adopt smoothness and moment conditions of the same form as in fw2016. The dominating function $\Psi((y_{it}, x_{it}^{\prime}))$ uniformly bounds both the regressors and all derivatives of the criterion function up to fourth order. Requiring $\sup_{it} \mathbb{E}\left[ (\Psi((y_{it}, x_{it}^{\prime})))^{\delta} \right]$ to be a.\,s.\ uniformly bounded over $N, T$ ensures that these envelopes have a finite $\delta$-th moment. Compared to fw2016, we allow $\delta \geq 6 + \nu$ instead of fixing $\delta = 8 + \nu$. The choice of $\delta$ affects, for example, the bandwidth of the debiased estimator defined below.
• The convexity condition has two parts. The first part requires the response-weighted expected second derivative of the criterion function with respect to $\pi_{it}$ to be a.\,s.\ bounded away from zero, uniformly over the $(i, t)$-pairs with $d_{it} = 1$ and over $N, T$. This ensures local strict convexity in $\phi$ and is needed for the asymptotic expansions. Together with (iv), the first part implies that $\inf_{\{(i, t) \colon d_{it} = 1\}} \mathbb{E}\left[ r_{it} \right] > 0$ a.\,s.\ uniformly over $i, t, N, T$, i.e., the $(i, t)$-pairs with $d_{it} = 1$ have positive response probabilities. The second part is a generalized non-collinearity condition. Together, the two parts imply strict convexity of the expected objective function over the relevant part of the parameter space.
• We place two conditions on the $\Phi$-measurable design component of the selection process. The first requires that, for each individual, the fraction of periods with $d_{it} = 1$ be bounded away from zero, uniformly over $N, T$. The second requires that any two periods $t \neq t^{\prime}$ share a fraction of units with $d_{it} = d_{it^{\prime}} = 1$ that is bounded away from zero, uniformly over $N, T$. The second condition can be interpreted as an overlap, or connectivity, condition. Without it, the design could split the $(i, t)$-pairs into blocks of periods with no units in common, producing a block-specific identification problem between $\alpha$ and $\gamma$. Because $d_{it}$ may itself depend on the unobserved effects or initial conditions, we impose both conditions almost surely. Since $\sum_{i = 1}^{N} d_{it} \geq \sum_{i = 1}^{N} d_{it} d_{it^{\prime}}$ for any $t \neq t^{\prime}$, the second condition implies $\inf_{t} N^{- 1} \sum_{i = 1}^{N} d_{it} \geq c_{O}$ for $T \geq 2$.
\end{enumerate}
Asymptotic Distribution
We write $\overline{\mathbb{E}}\left[ \cdot \right]$ for the probability limit of the enclosed expression as $N, T \rightarrow \infty$. By construction, $\mathbb{E}\left[ s_{it} \ddot{x}_{it} (\mathrm{d}^{2}_{\mathcal{C}} \psi)_{it} \right]$ is the vector of weighted residuals from the weighted least-squares problem (ref). Consequently,
equation[equation omitted — 368 chars of source]
theorem[Asymptotic Distribution of $\hat{\beta}$]
Let Assumption (ref) hold. Then,
\begin{equation*}
\sqrt{N T} \, (\hat{\beta} - \beta^{0}) \xrightarrow{d} \operatorname{\mathcal{N}}\big(- \tau^{\frac{1}{2}} \overline{W}_{\infty}^{- 1} \overline{B}_{\alpha, \infty} - \tau^{- \frac{1}{2}} \overline{W}_{\infty}^{- 1} \overline{B}_{\gamma, \infty}, \overline{W}_{\infty}^{- 1} \overline{\Sigma}_{\infty} \overline{W}_{\infty}^{- 1}\big) \, ,
\end{equation*}
where
\begin{align*}
& \overline{B}_{\alpha, \infty} \coloneqq \overline{\mathbb{E}}\left[ - \frac{1}{N} \sum_{i = 1}^{N} \frac{\sum_{t = 1}^{T} \sum_{t^{\prime} = t}^{T} \mathbb{E}\left[ s_{it^{\prime}} \ddot{x}_{it^{\prime}} (\mathrm{d}^{2} \psi)_{it^{\prime}} s_{it} (\mathrm{d}^{1} \psi)_{it} \right]}{\sum_{t = 1}^{T} \mathbb{E}\left[ s_{it} (\mathrm{d}^{2}_{\mathcal{C}} \psi)_{it} \right]} \right. \, + \\
& \quad \left. \frac{1}{2 \, N} \sum_{i = 1}^{N} \frac{\big\{\sum_{t = 1}^{T} \mathbb{E}\left[ s_{it} \ddot{x}_{it} (\mathrm{d}^{3}_{\mathcal{C}} \psi)_{it} \right]\big\} \Big\{\sum_{t = 1}^{T} \mathbb{E}\left[ s_{it} \big\{(\mathrm{d}^{1} \psi)_{it}\big\}^{2} \right]\Big\}}{\big\{\sum_{t = 1}^{T} \mathbb{E}\left[ s_{it} (\mathrm{d}^{2}_{\mathcal{C}} \psi)_{it} \right]\big\}^{2}} \right] \, , \\
& \overline{B}_{\gamma, \infty} \coloneqq \overline{\mathbb{E}}\left[- \frac{1}{T} \sum_{t = 1}^{T} \frac{\sum_{i = 1}^{N} \mathbb{E}\left[ s_{it} \ddot{x}_{it} (\mathrm{d}^{2} \psi)_{it} (\mathrm{d}^{1} \psi)_{it} \right]}{\sum_{i = 1}^{N} \mathbb{E}\left[ s_{it} (\mathrm{d}^{2}_{\mathcal{C}} \psi)_{it} \right]} \right. \, + \\
& \quad \left. \frac{1}{2 \, T} \sum_{t = 1}^{T} \frac{\big\{\sum_{i = 1}^{N} \mathbb{E}\left[ s_{it} \ddot{x}_{it} (\mathrm{d}^{3}_{\mathcal{C}} \psi)_{it} \right]\big\} \sum_{i = 1}^{N} \mathbb{E}\left[ s_{it} \big\{(\mathrm{d}^{1} \psi)_{it}\big\}^{2} \right]}{\big\{\sum_{i = 1}^{N} \mathbb{E}\left[ s_{it} (\mathrm{d}^{2}_{\mathcal{C}} \psi)_{it} \right]\big\}^{2}} \right] \, , \\
& \overline{W}_{\infty} \coloneqq \overline{\mathbb{E}}\left[ \frac{1}{N T} \sum_{i = 1}^{N} \sum_{t = 1}^{T} \mathbb{E}\left[ s_{it} (\mathrm{d}^{2}_{\mathcal{C}} \psi)_{it} \, \ddot{x}_{it} \, \ddot{x}_{it}^{\prime} \right] \right] \, , \\
& \overline{\Sigma}_{\infty} \coloneqq \overline{\mathbb{E}}\left[ \frac{1}{N T} \sum_{i = 1}^{N} \sum_{t = 1}^{T} \mathbb{E}\left[ s_{it} \big\{(\mathrm{d}^{1} \psi)_{it}\big\}^{2} \ddot{x}_{it} \, \ddot{x}_{it}^{\prime} \right] \right] \, ,
\end{align*}
assuming these limits exist.
remark[Theorem (ref)]
Our results nest those of fw2016 for balanced panels. To keep the discussion concise, we mainly comment on features specific to unbalanced panels.
\begin{enumerate}[(i)]
• As in fw2016, there are two bias terms, $\overline{B}_{\alpha, \infty}$ and $\overline{B}_{\gamma, \infty}$, arising from the estimation of $\alpha$ and $\gamma$, respectively. $\tau^{1 / 2} \overline{B}_{\alpha, \infty} / \sqrt{NT}$ is of order $T^{- 1}$, whereas $\tau^{- 1 / 2} \overline{B}_{\gamma, \infty} / \sqrt{NT}$ is of order $N^{- 1}$. The two terms are nearly symmetric in the indices $i$ and $t$. The exception is the double sum in the numerator of $\overline{B}_{\alpha, \infty}$, which has no counterpart in $\overline{B}_{\gamma, \infty}$ because of the conditional cross-sectional independence assumption. We therefore refer to
\begin{equation*}
\overline{B}_{\alpha, \infty}^{\star} \coloneqq \overline{\mathbb{E}}\left[ - \frac{1}{N} \sum_{i = 1}^{N} \frac{\sum_{t = 1}^{T} \sum_{t^{\prime} = t + 1}^{T} \mathbb{E}\left[ s_{it^{\prime}} \ddot{x}_{it^{\prime}} (\mathrm{d}^{2} \psi)_{it^{\prime}} s_{it} (\mathrm{d}^{1} \psi)_{it} \right]}{\sum_{t = 1}^{T} \mathbb{E}\left[ s_{it} (\mathrm{d}^{2}_{\mathcal{C}} \psi)_{it} \right]}\right]
\end{equation*}
as feedback bias. This is essentially a n1981-type bias. We use the term feedback bias (see, e.g., b2025) to stress that it can arise under general dynamic feedback, not only in dynamic models.
• For OLS and (pseudo-)Poisson maximum likelihood (PPML) estimators, the bias expressions simplify. Specific to these estimators is that $(\mathrm{d}^{2} \psi)_{it}$ does not depend on $y_{it}$ and is therefore $\mathcal{C}_{i}^{t}$-measurable. Since $s_{it}$ and $\ddot{x}_{it}$ are $\mathcal{C}_{i}^{t}$-measurable as well, the conditionally mean-zero score of Remark (ref) (iii) implies
\begin{equation*}
\mathbb{E}\left[ s_{it} \ddot{x}_{it} (\mathrm{d}^{2} \psi)_{it} (\mathrm{d}^{1} \psi)_{it} \mid \mathcal{C}_{i}^{t} \right] = \mathbf{0}_{K} \quad for all i, t, N, T \, .
\end{equation*}
In addition, $\sum_{i = 1}^{N} \mathbb{E}\left[ s_{it} \ddot{x}_{it} (\mathrm{d}^{3}_{\mathcal{C}} \psi)_{it} \right] = \mathbf{0}_{K}$ and $\sum_{t = 1}^{T} \mathbb{E}\left[ s_{it} \ddot{x}_{it} (\mathrm{d}^{3}_{\mathcal{C}} \psi)_{it} \right] = \mathbf{0}_{K}$ for all $i, t, N, T$. For OLS because $(\mathrm{d}^{2} \psi)_{it} = 1$ and $(\mathrm{d}^{3} \psi)_{it} = 0$. For PPML because $(\mathrm{d}^{3}_{\mathcal{C}} \psi)_{it} = (\mathrm{d}^{2}_{\mathcal{C}} \psi)_{it}$ so that both sums reduce to the normal equations (ref) of the weighted least-squares problem (ref). Consequently, $\overline{B}_{\alpha, \infty} = \overline{B}_{\alpha, \infty}^{\star}$ and $\overline{B}_{\gamma, \infty} = \mathbf{0}_{K}$.
• The selection process enters the asymptotic distribution as an additional source of heterogeneity. It can render subsamples heterogeneous and thereby invalidate jackknife approaches such as those of fw2016 and h2026_jackknife, which rely on forming homogeneous subsamples. We illustrate this in Section (ref).
• Imposing further restrictions on higher-order conditional moments, as in fw2016, would allow one to simplify the bias and variance components using Bartlett identities. For example, by the second Bartlett identity, $\mathbb{E}\left[ s_{it} \{(\mathrm{d}^{1} \psi)_{it}\}^{2} \right] = \mathbb{E}\left[ s_{it} (\mathrm{d}^{2}_{\mathcal{C}} \psi)_{it} \right]$ and $\mathbb{E}\big[s_{it} \{(\mathrm{d}^{1} \psi)_{it}\}^{2} \ddot{x}_{it} \, \ddot{x}_{it}^{\prime}\big] = \mathbb{E}\big[s_{it} (\mathrm{d}^{2}_{\mathcal{C}} \psi)_{it} \ddot{x}_{it} \, \ddot{x}_{it}^{\prime}\big]$ for all $i, t, N, T$. Hence, $\overline{\Sigma}_{\infty} = \overline{W}_{\infty}$, $\sum_{i = 1}^{N} \mathbb{E}\left[ s_{it} \{(\mathrm{d}^{1} \psi)_{it}\}^{2} \right] = \sum_{i = 1}^{N} \mathbb{E}\left[ s_{it} (\mathrm{d}^{2}_{\mathcal{C}} \psi)_{it} \right]$, and $\sum_{t = 1}^{T} \mathbb{E}\left[ s_{it} \{(\mathrm{d}^{1} \psi)_{it}\}^{2} \right] = \linebreak \sum_{t = 1}^{T} \mathbb{E}\left[ s_{it} (\mathrm{d}^{2}_{\mathcal{C}} \psi)_{it} \right]$ for all $i, t, N, T$. Furthermore, if $y_{it} \in \{0, 1\}$ for all $i, t, N, T$ and $(\hat{\beta}, \hat{\phi})$ are estimated by maximum likelihood, the simplifications implied by the second Bartlett identity follow automatically from the conditional mean assumption (see Remark 7 in cs2026).
• When the regressors and the selection process are strictly exogenous, $\overline{B}_{\alpha, \infty}^{\star} = \mathbf{0}_{K}$, and $\overline{B}_{\alpha, \infty}$ simplifies accordingly. In unbalanced panels, feedback bias can remain when the selection process is predetermined, even if the regressors are strictly exogenous.
\end{enumerate}
Bias Correction
The asymptotic distribution of $\hat{\beta}$ is not centered at zero, which invalidates standard inference. We therefore propose an analytical bias correction. It adjusts $\hat{\beta}$ using estimates of the two bias terms, which are constructed as sample analogs of the expressions in Theorem (ref), with the true parameter values replaced by the corresponding fixed effects estimates. We also briefly discuss alternative inference procedures based on the split-panel jackknife of fw2016 and the jackknife $t_{q}$-statistic of h2026_jackknife.
Let $(\widehat{\mathrm{d}^{\mathrm{r}} \psi})_{it}$ and $(\widehat{\mathrm{d}^{\mathrm{r}}_{\mathcal{C}} \psi})_{it}$ denote the sample analogs of $(\mathrm{d}^{\mathrm{r}} \psi)_{it}$ and $(\mathrm{d}^{\mathrm{r}}_{\mathcal{C}} \psi)_{it}$. Moreover, for each $k \in \{1, \ldots, K\}$, let
equation[equation omitted — 350 chars of source]
be the sample analog of (ref). The resulting $NT \times K$ matrix of fitted values is $\widehat{\mathfrak{X}}$, with the $it$-th row being $\hat{\mathfrak{x}}_{it} \coloneqq (\mu_{it}(\hat{\xi}_{1}), \ldots, \mu_{it}(\hat{\xi}_{K}))^{\prime}$. We set $\hat{\ddot{X}} \coloneqq X - \widehat{\mathfrak{X}}$ so that $\hat{\ddot{x}}_{it} = x_{it} - \hat{\mathfrak{x}}_{it}$, and we denote the debiased estimator by $\tilde{\beta}$.
Theorem (ref) establishes the asymptotic distribution of $\tilde{\beta}$ and the consistency of the variance estimators. To estimate $\overline{B}_{\alpha, \infty}$, we adapt the truncated spectral-density estimator of \textcites{hk2007}{hk2011} and the finite-sample adjustment of fw2016. We denote the corresponding bandwidth parameter by $h$.
theorem[Asymptotic Distribution of $\tilde{\beta}$]
Let Assumption (ref) hold. Then,
\begin{equation*}
\widehat{\Sigma} \xrightarrow{p} \overline{\Sigma}_{\infty} \quad and \quad \widehat{W}^{- 1} \xrightarrow{p} \overline{W}_{\infty}^{- 1} \, .
\end{equation*}
If, in addition, $h \asymp T^{\varsigma}$ for some $\varsigma \in (0, (\kappa - 1) / (2 \kappa))$, then,
\begin{equation*}
\sqrt{N T} \, (\tilde{\beta} - \beta^{0}) \xrightarrow{d} \operatorname{\mathcal{N}}\big(0, \overline{W}_{\infty}^{- 1} \overline{\Sigma}_{\infty} \overline{W}_{\infty}^{- 1}\big) \, ,
\end{equation*}
where $\tilde{\beta} = \hat{\beta} + \widehat{W}^{- 1} (T^{- 1} \, \widehat{B}_{\alpha} + N^{- 1} \, \widehat{B}_{\gamma})$ with
\begin{align*}
& \widehat{B}_{\alpha} \coloneqq - \frac{1}{N} \sum_{i = 1}^{N} \frac{\sum_{m = 0}^{h} \omega_{i}(m) \sum_{t = m + 1}^{T} s_{it} \hat{\ddot{x}}_{it} (\widehat{\mathrm{d}^{2} \psi})_{it} s_{i(t - m)} (\widehat{\mathrm{d}^{1} \psi})_{i(t - m)}}{\sum_{t = 1}^{T} s_{it} (\widehat{\mathrm{d}^{2}_{\mathcal{C}} \psi})_{it}} \, + \\
& \quad \frac{1}{2 \, N} \sum_{i = 1}^{N} \frac{\big\{\sum_{t = 1}^{T} s_{it} \hat{\ddot{x}}_{it} (\widehat{\mathrm{d}^{3}_{\mathcal{C}} \psi})_{it} \big\} \Big\{\sum_{t = 1}^{T} s_{it} \big\{(\widehat{\mathrm{d}^{1} \psi})_{it}\big\}^{2}\Big\}}{\big\{\sum_{t = 1}^{T} s_{it} (\widehat{\mathrm{d}^{2}_{\mathcal{C}} \psi})_{it}\big\}^{2}} \, , \\
& \widehat{B}_{\gamma} \coloneqq - \frac{1}{T} \sum_{t = 1}^{T} \frac{\sum_{i = 1}^{N} s_{it} \hat{\ddot{x}}_{it} (\widehat{\mathrm{d}^{2} \psi})_{it} (\widehat{\mathrm{d}^{1} \psi})_{it}}{\sum_{i = 1}^{N} s_{it} (\widehat{\mathrm{d}^{2}_{\mathcal{C}} \psi})_{it}} + \frac{1}{2 \, T} \sum_{t = 1}^{T} \frac{\big\{\sum_{i = 1}^{N} s_{it} \hat{\ddot{x}}_{it} (\widehat{\mathrm{d}^{3}_{\mathcal{C}} \psi})_{it} \big\} \sum_{i = 1}^{N} s_{it} \big\{(\widehat{\mathrm{d}^{1} \psi})_{it}\big\}^{2}}{\big\{\sum_{i = 1}^{N} s_{it} (\widehat{\mathrm{d}^{2}_{\mathcal{C}} \psi})_{it}\big\}^{2}} \, , \\
& \widehat{W} \coloneqq \frac{1}{N T} \sum_{i = 1}^{N} \sum_{t = 1}^{T} s_{it} (\widehat{\mathrm{d}^{2}_{\mathcal{C}} \psi})_{it} \, \hat{\ddot{x}}_{it} \, \hat{\ddot{x}}_{it}^{\prime} \, , \quad \widehat{\Sigma} \coloneqq \frac{1}{N T} \sum_{i = 1}^{N} \sum_{t = 1}^{T} s_{it} \big\{(\widehat{\mathrm{d}^{1} \psi})_{it}\big\}^{2} \, \hat{\ddot{x}}_{it} \, \hat{\ddot{x}}_{it}^{\prime} \, ,
\end{align*}
and $\omega_{i}(m) = \lvert \mathcal{T}_{i} \rvert / \sum_{t = m + 1}^{T} s_{it} s_{i(t - m)}$.
remark[Theorem (ref)]
Again, our results nest those of fw2016 for balanced panels. The analytical bias correction does not inflate the variance, i.e., the uncorrected and the debiased estimator have the same asymptotic variance. The correction factor $\omega_{i}(m)$ is adapted to allow for gaps in individual time series. For balanced panels it collapses to the factor proposed in fw2016, $\omega_{i}(m) = T / (T - m)$ for all $i, N$. Consistency of $\widehat{B}_{\alpha}$ requires $\varsigma \in (0, (\kappa - 1) / (2 \kappa))$. Within this range the choice of bandwidth is asymptotically immaterial. The rate-optimal exponent $\varsigma^{\ast} = (\kappa - 1) / (6 \kappa \delta)$ nonetheless lies near the lower end, favoring small exponents and hence rather small bandwidths relative to $T$. This is in line with fw2016, who recommend bandwidths no larger than four and who also recommend reporting results for several bandwidths. Increasing $\kappa$ widens the admissible range of exponents, with its upper limit approaching $1 / 2$, yet drives $\varsigma^{\ast}$ toward zero (recall $\delta = 2 \kappa + \nu$). Finally, the researcher need not know whether the regressors or the selection process are predetermined. The proposed bias correction is valid in either case.
remark[Jackknife Inference]
Alternatively, one may conduct inference using the split-panel jackknife of fw2016 or the jackknife $t_{q}$-statistic of h2026_jackknife, at the cost of an additional unconditional homogeneity assumption: $\{(y_{it}, x_{it}^{\prime}, s_{it}, \alpha_{i}, \gamma_{t}) \colon i \in \{1, \ldots, N\}, t \in \{1, \ldots, T\}\}$ is identically distributed across $i$ and strictly stationary across $t$ for all $N, T$. Unlike the rest of the paper, this assumption is not conditional on $\Phi$. A design $d_{it}$ that is non-random given $\Phi$ is unconditionally random whenever it depends on unobserved effects or initial conditions. Thus, $s_{it}$ is also unconditionally random. The assumption is sufficient to ensure that the biases are homogeneous across all subsamples. It rules out, for example, structural breaks or time trends in the data-generating processes of the observed variables (including the selection process) and the unobserved effects. Section (ref) examines the consequences of violating this assumption. The resampling strategies proposed in fw2016 and h2026_jackknife carry over unchanged to unbalanced panels.
Simulation Experiments
We study the finite-sample behavior of the uncorrected and debiased fixed effects estimators for the parameters of a dynamic probit model in unbalanced panels. We adapt the data-generating process of fw2016 for dynamic probit models to unbalanced panels. For debiasing, we consider analytical bias corrections with $h \in \{0, 1, 2\}$ and jackknife corrections, forming subsamples in the same way as for balanced panels. We report the bias relative to the true parameter value (in percent), the coverage rate of confidence intervals with a 95% nominal level, and the average length of these intervals.
The first $101$ periods, indexed by $t \in \{-100, \ldots, 0\}$, serve as a burn-in for $\gamma_{t}$, $y_{it}$, $x_{it}$, and $r_{it}$. Estimation uses $t \in \{1, \ldots, T\}$ only. We initialize the processes at $t = -100$ by drawing $y_{i(-100)}$ and $r_{i(-100)}$ from a Bernoulli distribution with a success probability of $1 / 2$ each, and $x_{i(-100)}$ from a standard normal distribution. For $i \in \{1, \ldots, N\}$ and $t \in \{-99, \ldots, T\}$,
align*[align* omitted — 384 chars of source]
where $\alpha_{i}, \gamma_{t} \sim \operatorname{\text{iid.}\;} \operatorname{\mathcal{N}}(0, 1 / 16)$, $u_{it} \sim \operatorname{\text{iid.}\;} \operatorname{\mathcal{N}}(0, 1)$, $v_{it} \sim \operatorname{\text{iid.}\;} \operatorname{\mathcal{N}}(0, 1 / 2)$, and $w_{it} \sim \operatorname{\text{iid.}\;} \operatorname{\mathcal{U}}(0, 1)$.
The design component is governed by an entry period $e_{i}$ and an exit period $f_{i}$, along with $d_{it} = \operatorname{\mathbbold{1}} \{ e_{i} \leq t \leq f_{i} \}$. Units with $4 \lvert \alpha_{i} \rvert > 1$ are present throughout, i.e., $e_{i} = 1$ and $f_{i} = T$. The remaining units enter later and stay for a fraction of the post-entry periods, with $e_{i} = \min\{ 1 + \lfloor T_{\text{half}}(1 - l_{i}) \rfloor, T_{\text{half}} \}$ and $f_{i} = \min\{ T, e_{i} + \lceil 7 \, (T - e_{i}) / 10 \rceil \}$, where $T_{\text{half}} = \lfloor T / 2 \rfloor$, $l_{i} = \{ F_{\mathcal{N}}(4 \alpha_{i}) - F_{\mathcal{N}}(-1) \} / \{ F_{\mathcal{N}}(1) - F_{\mathcal{N}}(-1) \}$, and $F_{\mathcal{N}}(\cdot)$ denotes the standard normal cumulative distribution function. The selection process is mixed $s_{it} = d_{it} r_{it}$. The design is a deterministic function of $\alpha_{i}$ and is therefore $\Phi$-measurable, whereas the response component is stochastic and predetermined.
We set $\beta_{y} = 0.5$ and $\beta_{x} = 1$, and generate samples with $N = 200$ and $T \in \{10, 20, 40\}$. The intercept $\rho_{0}$ of the response equation is calibrated so that, on average, half of the observations are missing: $\rho_{0} = 0.095$ for $T = 10$, $\rho_{0} = 0.260$ for $T = 20$, and $\rho_{0} = 0.350$ for $T = 40$. All results are based on $10{,}000$ simulated samples for each configuration.
The design component violates the unconditional homogeneity assumption underlying the jackknife corrections (see Remark (ref)). Because $\alpha_{i}$ is i.i.d. across $i$, and $e_{i}$ and $f_{i}$ depend on $i$ only through $\alpha_{i}$, the vector $(y_{it}, x_{it}^{\prime}, s_{it}, \alpha_{i}, \gamma_{t})$ remains identically distributed across $i$ for all $N$. It is not, however, strictly stationary across $t$. Since $f_{i} < T$ for every unit with $4 \lvert \alpha_{i} \rvert \leq 1$, $\mathbb{P}\left(d_{it} = 1\right)$ falls by $t = T$ to the share of units present throughout, $2 \, \{1 - F_{\mathcal{N}}(1)\} \approx 0.32$. Because $s_{it} = d_{it} r_{it}$, the selection probability $\mathbb{P}\left(s_{it} = 1\right)$ likewise peaks around $T_{\text{half}}$ and declines thereafter to a value strictly below $0.32$. Consequently, the subsamples formed across $i$, which correct $\overline{B}_{\gamma, \infty}$, remain homogeneous, whereas the subsamples formed across $t$, which correct $\overline{B}_{\alpha, \infty}$, do not.
table[table omitted — 3,381 chars of source]
Table (ref) reports the results. The uncorrected estimator is severely biased, particularly for $\beta_{y}$. The bias in $\beta_{x}$ is smaller but still distorts coverage. The analytical correction with a positive bandwidth removes most of the bias and restores coverage close to the nominal level. The correction with $h = 0$ ignores the feedback bias, which predominates for $\hat{\beta}_{y}$. Because the regressor $x_{it}$ is strictly exogenous, $h = 0$ removes most of the bias in $\hat{\beta}_{x}$. Nonetheless, positive bandwidths perform better, reflecting the feedback bias arising from both the lagged outcome and the predetermined selection process. Performance is similar across positive bandwidths once $T$ is moderately large. For $T = 10$, the smaller bandwidth ($h = 1$) clearly dominates. Both patterns are consistent with Remark (ref), which favors small bandwidths. The two jackknife corrections share the same point estimate and differ only in their standard errors. They overcorrect the bias in $\hat{\beta}_{y}$ because the design renders the subsamples formed across $t$ heterogeneous. The split-panel jackknife of fw2016 often undercovers, while the $t_{2}$-statistic of h2026_jackknife is more robust, achieving nominal coverage at the cost of substantially wider intervals. Across all estimators, bias and interval length shrink as $T$ increases, as predicted by asymptotic theory.
Concluding Remarks
We have studied fixed effects M-estimation for unbalanced panels under selection that is $\Phi$-measurable in its design component and stochastic in its response component. We do not assume missing-at-random as common in the literature. The response may depend on past outcomes, regressors, and unobserved effects. Feedback bias can arise from a predetermined selection process, even when the regressors are strictly exogenous. Moreover, selection that depends on unobserved effects can render subsamples heterogeneous, thereby invalidating jackknife corrections. Our analytical correction does not rely on such homogeneity. It is valid whether or not the regressors and the selection process are predetermined, and the researcher need not know which case applies. For applied work, we recommend a positive but small bandwidth and suggest reporting results for several bandwidths.
Extending the analysis to estimators with non-smooth criterion functions, such as quantile regression, is a natural direction for future work.
\printbibliography