Closed-form estimation and inference for panels with attrition and refreshment samples

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.

37,354 characters · 7 sections · 20 citation commands

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

abstract\onehalfspacing It has long been established that, if a panel dataset suffers from attrition, auxiliary (refreshment) sampling restores full identification under additional assumptions that still allow for nontrivial attrition mechanisms. Such identification results rely on implausible assumptions about the attrition process or lead to theoretically and computationally challenging estimation procedures. We propose an alternative identifying assumption that, despite its nonparametric nature, suggests a simple estimation algorithm based on a transformation of the empirical cumulative distribution function of the data. This estimation procedure requires neither tuning parameters nor optimization in the first step, i.e., it has a closed form. We prove that our estimator is consistent and asymptotically normal and demonstrate its good performance in simulations. We provide an empirical illustration with income data from the Understanding America Study. JEL Classification: C23 Keywords panel data, survey, attrition, selection, refreshment sample, two-step GMM

\onehalfspacing \frenchspacing

Introduction

Attrition in panel data is a widespread phenomenon. Units tracked over time may drop out of the sample for several reasons, including self-selection, non-survival, and increasing survey burden. When attrition is nonrandom, the ensuing bias in structural estimates is difficult to handle both theoretically and computationally. Complete debiasing is impossible without either strong (parametric or nonparametric) assumptions on the attrition process, such as the selection on observables, or the availability of auxiliary data. The latter often comes in the form of refreshment samples, i.e., extra random samples from the population in the drop-out period, see, e.g., deng2013handling,taylor2020evaluating,watson2021refreshment. Refreshment samples are available in many widely used survey panels such as the RAND Malaysian Family Life Survey (MFLS), the Medical Expenditure Panel Survey (MEPS), and the Current Population Survey (CPS). In a seminal paper, hirano2001combining proved that, with refreshment samples, identification is restored under a quasi-separability assumption on the attrition process, which they call additive nonignorability. Subsequent papers developed estimation and inference procedures relying on assumptions of varying strength. For example, bhattacharya2008inference proposed a sieve-based semiparametric, asymptotically normal estimator under additive nonignorability, while hoonhout2019nonignorable derived a two-step GMM procedure for the case of multi-wave panels. Recently, franguridi2024estimation suggested a computationally tractable estimation procedure using iterative proportional fitting (raking), and franguridi2024robust provided a debiased version of the raking estimator along with an influence function-based estimator of its asymptotic variance. Alternative approaches to estimation and inference with refreshment samples or other auxiliary data include hellerstein1999imposing,nevo2003using,d2010new,si2015semi,sadinle2019sequentially,franguridi2025inference,franguridi2025generalized, among others.

The aforementioned estimation techniques either require strong assumptions on the attrition process or lead to computationally challenging multistep estimators. Our main contribution is an estimation procedure that relies on transformations of the empirical cumulative distribution function of the data, and that requires neither functional optimization nor the choice of tuning parameters and works for both discrete and continuous data; we call this procedure closed-form. This procedure is particularly attractive in the presence of high-dimensional covariates due to its immunity to the curse of dimensionality and admits bootstrap inference. These advantages come at the cost of a slightly nonstandard quasi-separability assumption on the attrition process. We consider it a fair price for the simplicity of theoretical analysis, computational feasibility, and absence of tuning parameters.

The rest of the paper is organized as follows. (ref) introduces the model and derives the key identification result. (ref) presents our closed-form estimator, its asymptotic analysis, and the construction of confidence intervals. (ref) illustrates the performance of the estimator and the confidence intervals in Monte Carlo simulations. (ref) presents an empirical application. (ref) concludes.

Framework and identification

We follow the setup of hirano2001combining.

Let $Z_{it} = (X_{it}',Y_{it})' \in \mathbb{R}^{d}$ denote the stacked vector of covariates and outcomes for unit $i$ at time $t = 1,2$. We observe $Z_{it}$ for a random sample of units $i=1,\dots,n_1$ at time $t=1$. There is no initial nonresponse. However, at time $t=2$, the units may drop out of the sample. Let $W_i$ be the indicator of unit $i$ staying in the sample and suppose, without loss of generality, that the stayers are units $i=1,\dots,n_2$. In addition to this incomplete panel, we observe an auxiliary (refreshment) sample $Z_{i2}^r$, $i=1,\dots,n_r$, from an unconditional distribution of $Z_{i2}$. Let $F$ be the cumulative distribution function (CDF) of $Z = (Z_1,Z_2)$ (where we drop the unit subscript). We are interested in estimating and conducting inference for a parameter $\theta \in \Theta \subset \mathbb{R}^{d_\theta}$ defined by the moment conditions

align[align omitted — 128 chars of source]

where $m: \mathbb{R}^{2d} \times \Theta \to \mathbb{R}^{d_m}$ is a known moment function. We assume for simplicity that $d_m = d_\theta$, but our results can be generalized to an arbitrary number of moments.

This framework is very general and includes the estimands of interest in both linear and nonlinear panel data models. Although we focus on the classical case of the target parameter defined by the moment conditions (ref), the results of this paper hold as long as $\theta=\theta(F)$ is a Hadamard differentiable functional of $F$. We illustrate the broad applicability of our setup with a series of examples.

exam[linear regression with two-way fixed effects] Consider the standard linear regression with two-way fixed effects, \begin{align*} y_{it} = \alpha_i + f_t + x_{it}'\theta + \varepsilon_{it}, \quad i=1,\dots, n, \,\, t=1,2. \end{align*} A huge literature is devoted to the identification of the slope coefficients $\theta$ under various assumptions. When the covariates are strictly exogenous, one can use the within transformation $\ddot{\zeta}_{it} = \zeta_{it} - \frac{1}{n} \sum_{i=1}^n \zeta_{it} - \frac{1}{2}(\zeta_{i1}+\zeta_{i2}) + \frac{1}{2n} \sum_{i=1}^n (\zeta_{i1}+\zeta_{i2})$ to identify $\theta$ as the OLS coefficient in the transformed regression, \begin{align*} \operatorname{\mathbb{E}} \left( \ddot{y}_{it} - \ddot{x}_{it}'\theta \right) \ddot{x}_{it} = 0. \end{align*} This is the moment condition of the form (ref), and hence can be estimated under attrition when a refreshment sample is available in the second period. When $x_{it}$ contains lagged outcomes, at least three periods are needed to estimate the above model (for example, via the Arellano-Bond GMM). Our methodology can be extended to handle more than two periods along the lines of hoonhout2019nonignorable, but we leave it for future work.

exam[difference-in-differences] Consider the classical framework in which outcomes $y_{it}$ are tracked for individuals $i$ over periods $t=1,2$, and some individuals are treated in period $2$, which is denoted by $d_{i2}=1$. The standard diff-in-diff estimand is \begin{align*} DID := \operatorname{\mathbb{E}} \left[ y_{i2}-y_{i1} \,\vert\, d_{i2}=1 \right] - \operatorname{\mathbb{E}} \left[ y_{i2}-y_{i1} \,\vert\, d_{i1}=0 \right], \end{align*} which, under the parallel trends assumption, is equal to the average treatment effect on the treated, ATT. Formally, using the potential outcome notation, \begin{align*} ATT &:= \operatorname{\mathbb{E}}[y_{i2}(1) - y_{i2}(0) \,\vert\, d_{i2}=1 ] = DID. \end{align*} In practice, we estimate the ATT via the two-way fixed effects regression \begin{align*} y_{it} = \alpha_i + f_t + \beta d_{it} + \varepsilon_{it}. \end{align*} This model can be estimated under attrition in the second period and the availability of a refreshment sample as discussed in the previous example by defining $z_1=(y_1,d_1), z_2=(y_2,d_2)$.

exam[quantile treatment effects] Consider the panel data with $T = 3$ periods, where the treatment only occurs in the last period. The data are a random sample from $(y_1,y_2,y_3,d)$, where $y_t$ is the observed outcome at time $t$ and $d \in \{0,1\}$ is the indicator of treatment. Suppose that the object of interest is the quantile treatment effect on the treated \begin{align*} QTT(\tau) = F_{y_3(1)|d=1}^{-1}(\tau) - F_{y_3(0)|d=1}^{-1}(\tau). \end{align*} The first term is identified directly from the data. For the second term, callaway2019quantile show that, under their assumptions of distributional parallel trends and copula stability, \begin{align*} F_{y_3(0)|d=1}(y) = \operatorname{\mathbb{P}} & \left[ F_{\Delta y_3|d=0}^{-1}\left( F_{\Delta y_{2}|d=1}(\Delta y_{2}) \right) \le y - F_{y_{2}|d=1}^{-1} \left(F_{y_{1}|d=1}(y_{1}) \right) \,\vert\, d=1\right]. \end{align*} The distributions $F_{\Delta y_3|d=0}$, $F_{\Delta y_2|d=1}$, $F_{y_2|d=1}$, and $F_{y_1|d=1}$ are identified directly from the data, and hence the second term in the QTT is also identified. Therefore, we can write the QTT as a (complicated) nonlinear functional of the joint distribution of $(y_1,y_2,y_3,d)$. Suppose some units may drop out from the sample in period $t=3$, but a refreshment sample is available. Then the QTT model fits our framework with the first two periods combined into one, i.e., with $z_1=(y_1,y_2)$ and $z_2 = (y_3,d)$.

Now, we introduce our main identifying assumption for point identification of the joint distribution of $Z_1$ and $Z_2$. Let $F^w$, $F_1$, and $F_2$ be the CDFs of $(Z_1,Z_2)|W=1$, $Z_1$, and $Z_2$, respectively. These distributions can be readily estimated from the balanced panel (retaining stayers only), the first-period sample, and the refreshment sample, respectively. The key identity relating the target distribution $F$ to the data is

align[align omitted — 151 chars of source]

The weight $\operatorname{\mathbb{P}}(W=1) / \operatorname{\mathbb{P}}(W=1|Z_1\le z_1,Z_2\le z_2)$ is not identified without further restrictions. To see why, notice that this object is an unrestricted function of the joint distribution of $Z_1,Z_2$, while the information available for its identification is only the two marginal CDFs $F_1$ and $F_2$. To close the gap, we impose a separability assumption on the weight.\footnote{If no assumptions are imposed on the attrition process, any distribution with marginals $F_1$ and $F_2$ is consistent with the data. Hence, in most cases, the partial identification approach will not lead to informative bounds on the structural parameter.}

assumption$\operatorname{\mathbb{P}}(W=1|Z_1\le z_1,Z_2\le z_2) = G(k_1(z_1)+k_2(z_2))$ for a known, differentiable, strictly increasing function $G: \mathbb{R} \to (0,\infty)$ and some unknown functions $k_1: \mathbb{R}^{d} \to \mathbb{R}$, $k_2:\mathbb{R}^{d} \to \mathbb{R}$.\footnote{We believe that choosing the link function $G$ is not a critical decision in practice, as documented, e.g., in little1991models,hirano2001combining. However, formally comparing the classes of selection mechanisms captured by (ref) when varying the link function is outside of the scope of this paper.}

This assumption is compatible with the missing-completely-at-random condition ($k_1=k_2=const$). It leads to an explicit identification of $k_1$ and $k_2$, see (ref) below. Besides, it neither implies nor is implied by the analogous assumption on $\operatorname{\mathbb{P}}(W=1|Z_1=z,Z_2=z)$ in hirano2001combining, viz.

align[align omitted — 131 chars of source]

for some known link function $\tilde G$ and unrestricted functions $\tilde k_1,\tilde k_2$. The following example illustrates this point.

examSuppose $Z_1,Z_2 \in [0, 1]^2$ with density $f(z_1, z_2) = z_1 + z_2 $ and the conditional probability of staying is given by \begin{align*} \operatorname{\mathbb{P}}(W=1|Z_1=z_1,Z_2=z_2) = az_1^2 + b z_1 z_2 + a z_2^2. \end{align*} for some constants $a,b$. Using the formula \begin{align*} \operatorname{\mathbb{P}}(W=1|Z_1\le z_1,Z_2\le z_2)& = \frac{1}{F(z_1,z_2)} \int_{-\infty}^{z_1} \int_{-\infty}^{z_2} \operatorname{\mathbb{P}}(W=1|Z_1=t_1,Z_2=t_2) f(t_1,t_2) \, dt_1 dt_2=\\ & =\frac{3az_1^3+2(a+b)z_1^2z_2+2(a+b)z_1z_2^2+3az_z^3}{6(z_1+z_2)}, \end{align*} it is easy to show that (a) when $a=2/11$, $b=7/11$, (ref) holds with $G(x)=x^2/11$, but (ref) does not hold; (b) when $a=1/2$, $b=0$, (ref) does not hold, while (ref) holds; (c) when $a=0$, $b=1$, both (ref) and (ref) hold. Put differently, even when the alternative assumption (ref) is imposed, our assumption is still valid for a class of DGPs with nontrivial attrition processes.

The next example shows that the left-hand side of (ref) can be a valid CDF under (ref).

examSuppose $Z_1,Z_2$ are scalar random variables and $G(x)=\exp(x)$. We have \[ F(z_1,z_2) \propto \frac{F^w(z_1,z_2)}{\exp(k_1(z_1)+k_2(z_2))}, \] or, taking the mixed partial derivative, \[ f(z_1,z_2) \propto \frac{(1-k_1'(z_1)D_{z_1}F^w )(1- k_2'(z_2)D_{z_2}F^w ) + (f^w-1)}{\exp(k_1(z_1)+k_2(z_2))}. \] This is a proper density if the numerator is nonnegative. Since $D_{z_1}F^w$, $D_{z_2} F^w$, and $f^w$ are nonnegative functions, a simple sufficient condition is for $k_1$ and $k_2$ to be decreasing. For a more concrete example, consider $Z_1,Z_2 |W=1 \sim \text{ iid U}[0,1]$ and $k_1(z_1)=a + c_1 z_1$, $k_2(z_2)=b+c_2z_2$. Then the aforementioned condition holds whenever $c_1,c_2 \le 1$. Normalization is achieved by an appropriate choice of $a$ and $b$.

(ref) leads to a closed-form identification result for the unknown function $k_1(z_1)+k_2(z_2)$ which informs an estimation procedure based on empirical CDFs. This is in contrast to hirano2001combining that identifies $\tilde k_1(z_1)+\tilde k_2(z_2)$ implicitly as a solution to a nonlinear functional optimization problem which makes the estimation procedure hard to analyze and implement, see franguridi2024robust.

theoremUnder (ref), the target CDF can be written as \begin{align} F(z_1,z_2) = \Phi\left(\operatorname{\mathbb{P}}(W=1),F_1(z_1),F_2(z_2),F_1^w(z_1),F_2^w(z_2),F^w(z_1,z_2) \right), \end{align} where the function $\Phi$ is defined by \begin{align} \Phi(p,F_1,F_2,F_1^w,F_2^w,F^w) = \frac{p F^w}{G\left(G^{-1}\left(\frac{p F_1^w}{F_1} \right) + G^{-1}\left(\frac{p F_2^w}{F_2} \right) - G^{-1}\left(p\right) \right)}. \end{align}

proofDenote $p=\operatorname{\mathbb{P}}(W=1)$. Equation (ref) implies \begin{align*} G(k_1(z_1)+k_2(z_2)) = \frac{pF(z_1,z_2|W=1)}{F(z_1,z_2)}. \end{align*} Since $G$ is strictly increasing, the inverse function $G^{-1}$ exists. Plugging in $z_1=\infty$ and/or $z_2=\infty$, we get \begin{align*} k_1(z_1) + k_2(\infty) &= G^{-1}\left(\frac{p F_1(z_1|W=1)}{F_1(z_1)} \right), \\ k_1(\infty) + k_2(z_2) &= G^{-1}\left(\frac{p F_2(z_2|W=1)}{F_2(z_2)} \right), \\ k_1(\infty)+k_2(\infty) &= G^{-1}\left(p\right). \end{align*} Normalizing (say) $k_2(\infty)=0$, we obtain, for all $z_1,z_2$, \begin{align*} k_1(z_1) &= G^{-1}\left(\frac{p F_1(z_1|W=1)}{F_1(z_1)} \right), \\ k_2(z_2) &= G^{-1}\left(\frac{p F_2(z_2|W=1)}{F_2(z_2)} \right) - G^{-1}\left(p\right). \end{align*} Substituting these expressions into (ref) completes the proof.

Estimation and inference

(ref) suggests a two-step estimation procedure. In the first step, we estimate $F$ by plugging in the empirical CDF estimators of $F_1,F_2,F_1^w,F_2^w,F^w$. In the second step, we substitute this estimator into the moment conditions (ref) and compute $\hat\theta$ that sets these moment conditions to zero.

To describe this procedure formally, define

align[align omitted — 421 chars of source]

Then the plug-in estimator of $F$ is

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

We now use this estimator to calculate the sample analog of the moment condition (ref) as the integral with respect to the distribution induced by $\hat F$. Because its arguments are empirical CDFs, $\hat F$ is a piecewise constant function. Let $\hat{\mathcal{Z}} \subset \mathbb{R}^{2d}$ be the finite set of its discontinuity points. In our case, \[ \hat{\mathcal{Z}} := \hat{\mathcal{Z}}_1 \times (\hat{\mathcal{Z}}_2 \cup \hat{\mathcal{Z}}_2^r), \] where $\hat{\mathcal{Z}}_1=\{z_{1i}, \, i\in n_1\}$, $\hat{\mathcal{Z}}_2=\{z_{2i}, \, i\in n_2\}$, and $\hat{\mathcal{Z}}_2^r=\{z_{2i}^r, \, i\in n_r\}$ are the respective empirical supports. Then the induced “probability” (“jump size”) of $\hat F$ at any $\zeta=(\zeta_1,\dots,\zeta_{2d})\in\hat{\mathcal{Z}}$ is\footnote{If $\hat F$ were the empirical CDF corresponding to $n$ distinct points $\hat{\mathcal{Z}}=\{\zeta_1,\dots,\zeta_n\}$, then the formula would yield $\hat f(x_i)=1/n$ for all $i=1,\dots,n$.}

align[align omitted — 201 chars of source]

where $h_k$ is any positive scalar that is less than the minimal gap between the observations of the $k$-th component of $\zeta$, $0<h_k < \min_{i,j:\, \zeta_{k,i} \neq \zeta_{k,j}}|\zeta_{k,i}-\zeta_{k,j}|$.\footnote{The definition of $\hat f$ does not depend on the choice of $h_k$, and so $h_k$ is not a tuning parameter.} This is a discrete version of the formula

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

which links a CDF $F$ to its PDF $f$ in case the former is differentiable. We emphasize, however, that the formula (ref) works regardless of whether the data are discrete or continuous. Notice that, because $\hat F$ does not have to be a valid empirical CDF, $\hat f$ may take negative values in finite samples.\footnote{Rearrangement of $\hat{F}(z_1,z_2)$ similar to chernozhukov2009improving does not guarantee that the rearranged function is a proper CDF.} This does not affect the asymptotic properties of the resulting estimator.

The complete estimation procedure is as follows.

Algorithm.

enumerate• Calculate $\hat p$ and the empirical CDFs $\hat F_1, \hat F_2, \hat F_1^w, \hat F_2^w, \hat F^w$; • Plug in to obtain $\hat F = \Phi(\hat p, \hat F_1, \hat F_2, \hat F_1^w, \hat F_2^w, \hat F^w)$; • Calculate the jump size $\hat f(z_{1},z_{2})$ at points $(z_{1},z_{2}) \in \hat{\mathcal{Z}}$ according to formula (ref); • Set $\hat\theta$ such that\footnote{Another way to implement the second step is to draw a random sample from the discrete distribution defined by (a trimmed version of) $\hat f$ and estimate $\theta$ by the conventional GMM based on this sample.} \begin{align} \int m(z_1,z_2,\hat\theta) \, d\hat F(z_1,z_2) = \sum_{(z_1,z_2) \in \hat{\mathcal{Z}}} m(z_{1},z_{2},\hat\theta) \hat f(z_{1},z_{2})=0. \end{align}

This procedure has several advantages. First, it does not require choosing any tuning parameters despite relying on a semiparametric (ref). Second, the first-step estimator has an explicit formula and only depends on empirical CDFs, which makes it suitable for high-dimensional data. Finally, because the first-step estimator is a smooth transformation of empirical CDFs, which jointly converge to a known Gaussian process by Donsker's theorem, this procedure allows for the delta method bootstrap, as the following theorem shows.

To prove consistency of $\hat\theta$, we follow newey1994large. First, under two different sets of assumptions on $m$, we show convergence in probability of $\int m(z;\theta)\, d\hat{F}(z)$ to $\int m(z;\theta)\, dF(z)$ uniformly over $\theta$. This is an analog of the uniform law of large numbers (ULLN) with sample averaging replaced by integration with respect to an estimated CDF.

assumption\begin{enumerate} • The parameter $\theta$ belongs to a compact set $\Theta \subset \mathbb{R}^{d_{\theta}}$; • The true value $\theta$ is a unique solution of (ref) in $\Theta$. \end{enumerate}

assumption• The true distribution $F$ has bounded support.

assumptionThe moment function satisfies the following: \begin{enumerate} • $m(Z;\theta)$ is a.s. continuous in $\theta\in\Theta$ • There exists $M_1>0$ such that $\Vert m(\cdot ;\theta)\Vert_{Lip}\le M_1$ for all $\theta\in\Theta$, where $$\Vert m(\cdot;\theta)\Vert_{Lip}\equiv \sup_{Z\ne\tilde{Z}} \frac{ \Vert m(Z;\theta)-m(\tilde{Z};\theta)\Vert_{1}}{\Vert Z-\tilde{Z}\Vert_{1}}.$$ \end{enumerate}

(ref) imposes smoothness on the moment function $m(Z,\theta)$, viz. $m(Z,\theta)$ is Lipschitz continuous w.r.t. $Z$ for any fixed $\theta\in\Theta$ with a uniformly bounded Lipschitz constant.

lemSuppose (ref) hold. Then $M(\theta) = \int m(z,\theta) \, dF(z) $ is continuous and \begin{equation} \sup_{\theta\in\Theta}\left\| \int m(z;\theta)d\hat{F}(z) - \int m(z;\theta)dF(z)\right\|\overset{p}{\to} 0. \end{equation}

Proof can be found in (ref).

To allow for non-smooth moment functions, such as $m(z,\theta)=1(z\le\theta)$, one can apply (ref).

assumptionThe moment function satisfies the following: \begin{enumerate} • $m(z;\theta)$ is of bounded variation as a function of $z$ for each $\theta\in\Theta$; • For every $\theta\in\Theta$, the set $\{z\,\vert\,\lim_{\gamma\to\theta}m(z;\gamma) = m(z;\theta)\}$ has probability 1 w.r.t. $F$; • There exists a function $d(z)$ such that $\Vert m(z;\theta) \Vert \le d(z)$ for all $\theta\in\Theta$ and $\int d(z)\, dF(z)<\infty;$ • Function $u(z;\theta, d) = \sup_{\vert\gamma-\theta\vert\le d}\vert m(z;\gamma)-m(z;\theta)\vert$ is of bounded variation as a function of $z$ for each $\theta\in\Theta$ and for each $d\le \varepsilon$ for some $\varepsilon>0.$ \end{enumerate}

(ref) is similar to the classical statement of the uniform law of large numbers, see Lemma 2.4 in newey1994large, with the additional high-level requirement on the function $u(z;\theta, d)$ to be of bounded variation.

lemUnder (ref), the conclusions of (ref) hold.

Proof can be found in (ref). It follows from Lemma 1 of tauchen1985diagnostic with an adjustment to the fact that $\hat{F}$ is not an empirical CDF.

theorem[Consistency] If the assumptions of (ref) or (ref) hold, then $$\hat{\theta}\overset{p}{\to} \theta.$$

(ref) follows from (ref), (ref), and an argument similar to the one in Theorem 2.1 of newey1994large.

Finally, we establish $\sqrt{n}$-consistency and asymptotic normality of $\hat\theta$ and the validity of the nonparametric bootstrap for inference on $\theta$. To this end, we impose the following assumption that ensures the Hadamard differentiability of $\theta$ as a functional of the distribution of the data.

assumption\begin{enumerate} • There exists $\varepsilon>0$ and a neighborhood $\tilde \Theta$ of $\theta$ such that \[ M(\tilde \theta,\tilde F) := \int m(z,\tilde \theta) \, d \tilde F(z) \] is continuously differentiable in $\tilde \theta \in \tilde \Theta$ for all $\tilde F$ such that $\sup_{z \in \mathbb{R}^{2d}}|\tilde F(z)-F(z)| < \varepsilon$. • The Jacobian $\frac{\partial}{\partial \theta} M(\theta,F)$ exists and is nonsingular. \end{enumerate}

(ref).(ref) holds if $m(z,\tilde \theta)$ is continuously differentiable in $\tilde \theta$ for $F$-a.e. $z$, without further restrictions on the distribution $F$. It also holds in the “minimum distance” case $m(z,\tilde \theta)=\tilde m(z)-\tilde \theta$. (ref).(ref) is standard and posits local identification.

theorem[Inference] Suppose that (ref) holds. Then $\sqrt{n}(\hat\theta - \theta)$ converges weakly to a normal distribution that can be consistently estimated by the nonparametric bootstrap.

Proof can be found in (ref).

Monte Carlo simulation

table[table omitted — 1,346 chars of source]

table[table omitted — 1,282 chars of source]

In this section, we illustrate the performance of our two-step estimator and the associated bootstrap confidence intervals in simulations. We abstract away from both conditioning on the covariates and nonlinear moments, as in the examples above, and employ two simple data-generating processes (DGPs), one with discrete variables and one with continuous variables. Both DGPs employ the logit link function $G(x) = (1+e^{-x})^{-1}$ for the conditional selection probability.\footnote{We explore the sensitivity of our methodology to misspecification of the link function in (ref).}

Discrete data

Our first DGP is a discrete Markov process with scalar outcomes. We posit that $Z_1$ has the uniform distribution over $\{1,\dots,m\}$, where $m \in \{5,10,20\}$, and $Z_2$ is determined from $Z_1$ given a transition matrix with positive elements. The attrition functions are

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

where the constants $c_1,c_2$ and the transition matrix are chosen so that the unconditional attrition rate is $30\%$. The target parameter is $\theta(m) = \operatorname{\mathbb{P}}_m(Z_2=1|Z_1=1)$ with true values

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

Since the parameter of interest is a probability that decreases with the number of support points, we present the bias, root mean-squared error (RMSE), and mean absolute deviation (MAE) as shares of the true value of $\theta(m)$.

Continuous data

Our second DGP is a copula model for continuous variables $Z_1$ and $Z_2$. The (target) joint distribution is modeled as a Gaussian copula

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

where $\Phi_\nu$ is a bivariate Gaussian CDF with zero means, unit variances, and correlation $\nu \in \{0.2,0.5,0.8\}$, and $\Phi$ is the CDF of $N(0,1)$. The functions in the attrition mechanism are

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

where the constants $c_0,c_1,c_2$ are chosen so that the overall attrition rate is about $70\%$. The target parameter is $\theta(\nu)=\operatorname{\mathbb{E}}_{\nu}[Z_1 Z_2]$ with the true values

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

Tables (ref) and (ref) present the simulation results. We report bias, root mean squared error (RMSE), and mean absolute error (MAE) of our estimator and of the “naive” estimator that ignores attrition and only uses the balanced panel. We also calculate empirical coverages of the bootstrap confidence interval for nominal confidence levels $1-\alpha \in \{0.90,0.95,0.99\}$. We do not report the empirical coverages for the naive estimator due to its inconsistency.

Unsurprisingly, the naive estimator exhibits larger biases across various specifications and sample sizes, while our (consistent) estimator performs well, with both the bias and the RMSE decreasing with the sample size. The confidence interval coverage is close to nominal for both the discrete DGP and the continuous DGP.

Empirical illustration

We illustrate our methodology in a simple model of the change of income over the life course. We use survey data from the Understanding America Study (UAS) conducted by the uas2018, a large household panel collected and maintained by the USC Center for Economic and Social Research. We estimate the following model,

align[align omitted — 147 chars of source]

where $\text{income}_{it}$ is the inverse hyperbolic sine of the household income.\footnote{This transformation approximates the logarithm of income for sufficiently large values of income while retaining zero income values.}

We use UAS wave 14 (year 2018) as the first period and UAS wave 15 (year 2020) as the second period.\footnote{Available at https://uasdata.usc.edu/page/Comprehensive+File+And+Panel+Dataset} The number of households responding in the first period is $8145$, out of which $6432$ households also respond in the second period, with an attrition rate $21\%$. There are 3983 households sampled in the second period that are not part of the sample in the first period. These households serve as the refreshment sample. We set the link function to be logistic, $G(x) = 1/(1+e^{-x})$.

(ref) contains the estimates and the bootstrap standard errors for our procedure using the refreshment sample and for the naive procedure using only the balanced panel. Neither of the estimates is significant at the 5% nominal level. However, the effect size is substantial: at the average level of income of the working population in 2018 (which is \$60,465), an additional year of life is associated with an average increase of $\sinh(\operatorname{asinh}(\$60,465)+\hat\theta_1+\hat\theta_2)-\$60,465 =\$24,824$ in annual household income.

table[table omitted — 420 chars of source]

Conclusion

In most panel datasets encountered in practice, some units are only observed up to a non-terminal period. Such attrition may lead to significant bias in the estimates and invalid inference. Fortunately, the availability of refreshment samples allows one to remove the bias and restore valid inference under reasonably weak assumptions on the attrition process.

In this paper, we introduce one such assumption that leads to simple, computationally feasible estimators that admit bootstrap inference. We hope this paper may serve as the first step towards developing feasible estimation and inference algorithms for data structures with nonrandom missingness and auxiliary information.

Topics for further research include an extension to multiple periods, more general missingness patterns, or high-dimensional covariates; dealing with initial non-response; theoretical and computational analysis of behavior under misspecification; establishing a semiparametric efficiency bound; and deriving optimal estimators.

Acknowledgements

We appreciate valuable comments and suggestions from Isaiah Andrews, Tim Armstrong, Arie Kapteyn, Sergey Lototsky, Kirill Ponomarev, John Pepper, Geert Ridder, Fedor Sandomirskiy, Alexander Shapoval, seminar participants at Princeton University and University of Rochester, and conference participants at 16th Greater New York Metropolitan Area Econometrics Colloquium, 2024 North American Summer Meeting of the Econometric Society, Western Economic Association International 99th Annual Conference, and Econometric Society Summer School in Dynamic Structural Econometrics. All errors and omissions are our own.