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.
117,967 characters · 17 sections · 64 citation commands
Bounds on Average Effects in Discrete Choice Panel Data Models
\thispagestyle{empty}
Keywords: Panel data, discrete choice, average effects, set identification, outer bounds, female labor force participation.
JEL Codes: C01, C23, C25
Panel data models with individual-specific effects make it possible to control for unobserved heterogeneity and confounding due to omitted variables that are constant over time. Nonlinear models are required to correctly describe discrete outcomes, and the main complication in such nonlinear panel models is the unknown distribution of unobserved heterogeneity, which constitutes an infinite-dimensional parameter. The fixed effects approach leaves this distribution unspecified, eliminating misspecification concerns (as opposed to the correlated random effects approach which models this distribution parametrically). However, lack of sufficient time-series variation in short panels means that this unknown distribution remains set-identified. An important consequence of this is a general lack of point-identification of average effects. While it is theoretically possible to recover the sharp identified set for average effects, in empirically relevant panel dimensions this often becomes an infeasible task due to a curse of dimensionality. This is a serious issue because average effects are typically the ultimate object of interest, especially from the policy perspective. In this paper, considering a general semiparametric setting, we propose alternative outer bounds which are simple to obtain and remain free of the curse of dimensionality in empirically relevant settings.
Formally, let $Y_i=(Y_{i1},\ldots,Y_{iT})$ be the vector of observed outcomes for individual $i=1,\ldots,n$, where $T$ is the number of time periods and $n$ is the number of cross-sectional units. Throughout, we assume that $n\to\infty$ but $T$ remains fixed. The semiparametric panel models we consider in this paper describe the distribution of $Y_i$ conditional on a vector of observed conditioning variables $Z_i$ as
Here, $f(y_i \, | \,z_i,a_i;\beta)$ is the distribution of $Y_i$ conditional on $Z_i$ and the (vector of) unobserved individual effects $A_i$, and it is assumed to be known up to the finite dimensional parameter $\beta$. The distribution of $A_i$ conditional on $Z_i$, given by $\pi(a_i|z_i)$, is left unrestricted. Both $\beta$ and $\pi = \pi(a_i|z_i)$ are unknown. The vector of conditioning variables usually consists of observed covariates $(X_{i1},\ldots,X_{iT})$ and/or initial conditions $(Y_{i0},Y_{i,-1},\ldots)$. Given the true distribution of $Y_i$ conditional on $Z_i$, the identified set for the model parameters consists of all pairs $(\beta,\pi)$ that satisfy (ref).
In empirical research, the ultimate object of interest is generally an average effect of the form
where $m(Z_i,A_i,\beta)$ is some function of interest. The exact choice of $m(\cdot,\cdot,\cdot)$ may vary from application to application, leading to different definitions of $\overline{m}$; see, among others, Chamberlain(84), BlundellPowell(03),BlundellPowell(04), AltonjiMatzkin(05), Wooldridge(05),Wooldridge(05chp), BesterHansen(09), GrahamPowell(12), HoderleinWhite(12).\footnote{A different quantity of interest, which we will not consider, is the quantile structural function of ImbensNewey(09).} AbrevayaHsu21 provide a detailed discussion of different average effects used in the literature.
The average effect in (ref) can be rewritten as
which clearly depends on $(\beta,\pi)$. In discrete choice models those model parameters (and in particular $\pi$) are usually only partially-identified, implying that $\overline m$ is also typically only partially-identified.
HonoreTamer06 and CFHN13 provide methods for obtaining the identified set when covariates are discrete. More recently, there has been an increased interest in the identification and estimation of average effects in various settings; see, e.g., AC21, DDL21, LiuPoirierShiu21, BM22, BMS22, and DGK24.\footnote{Lack of point-identification of $\pi(a_i|z_i)$ does not invariably lead to set-identification of $\overline m$. An interesting contribution in this vein is by AC21 who obtain point-identification of the average effect with respect to the lagged dependent variable in a dynamic logit model. However, such case-specific results usually remain an exception. A different route is to obtain point-identification of average effects under additional restrictions on the data generating process as in LiuPoirierShiu21. In contrast to these approaches, our aim is to provide a method which applies to an arbitrary function $m(Z_i,A_i,\beta)$ in a generic semiparametric framework. }
Unfortunately, obtaining the sharp identified set is often practically infeasible for sample sizes typically encountered in applications, due to a curse of dimensionality. This is because obtaining the sharp identified set for $\overline m$ typically requires estimates of the conditional probabilities $f_{Y|Z}(y_i|z_i)$. Since $Z_i$ usually contains a time-vector of (multiple) covariates, the curse of dimensionality is obvious for continuous covariates. However, even with discrete covariates the number of conditional probabilities that would need to be estimated is usually large. Suppose, for example, $Y_{it}, X_{it} \in \{0,1\}$, and that $Z_i=(X_{i1},\ldots,X_{iT})$. This implies $2^{2T}$ different conditional probabilities $f_{Y|Z}(y_i|z_i)$; for, say, $T=5$ this yields $1,024$ conditional probabilities. Estimation of objects of such numbers would require a much larger cross-sectional sample size than available in the majority of applications.\footnote{ Some general inference frameworks, like the ones described in chen2011sensitivity, are in principle applicable to models of the form (ref) and (ref). However, we are not aware of any framework that addresses the main statistical challenge that we are facing -- namely that $Z_i$ is high-dimensional and every realization of $Z_i$ is unique in standard panel applications. }
Motivated by this issue, we propose alternative bounds on the average effect $\overline m$ which can be feasibly obtained in realistic data settings. Our proposal is based on finding appropriate functions $L(Z_i,Y_i,\beta)$ and $U(Z_i,Y_i,\beta)$ such that
We show that suitable functions $L(\cdot,\cdot,\cdot)$ and $U(\cdot,\cdot,\cdot)$ can be obtained by solving an appropriate linear program for each realized value of $Z_i$. Asymptotically valid lower and upper bounds are then given by
respectively, for some appropriate estimator $\widehat \beta$, assuming that $\beta$ is point-identified. We prove the validity of the proposed bounds and provide asymptotically valid inference methods on $\overline m$. Our approach allows for discrete, as well as continuous covariates. We also provide computationally feasible methods for obtaining the suggested bounds. Importantly, these do not require searching over the space of possible distributions for $\pi(a_i|z_i)$, but only over the domain of $A_i$ itself. Consequently, implementation of our method is computationally straightforward and fast.
Our proposal differs from the existing literature in several ways. Firstly, we do not propose a different approach to obtaining the sharp identified set for $\overline m$; rather, we obtain outer bounds on this set. This has the virtue of avoiding the curse of dimensionality associated with the conditioning variable $Z_i$. Indeed, our outer bounds can be feasibly obtained at standard sample sizes even if the vector of conditioning variables $Z_i$ is continuous, or high-dimensional, or takes on many different values within the sample. Secondly, given our general semiparametric setting, the proposed method can easily be applied to different models (and functions $m(Z_i,A_i,\beta)$) of interest, such as the static logit, dynamic logit or the more complicated random coefficient logit models.
DDL21 propose an alternative method to achieve inference on $\overline m$. Their paper initially focuses on inference on the sharp identified set, but they also consider “outer bounds” (different from ours) that avoid non-parametric estimation of intermediate objects, similar in spirit to our results here. However, their approach currently only applies to static logit and ordered logit models (and for several choices of average effects), while in this paper we consider general models of the form (ref) (and more general average effects of the form (ref)).
Throughout the paper, we consider the case where $\beta$ is point-identified. However, our approach can easily be extended to models with partially-identified $\beta$, and we suggest two different extensions in the Supplementary Appendix. We however also note that methods for point-estimation of $\beta$ are well-established in the literature, and these methods are regularly used by applied researchers. Indeed, for essentially every type of discrete outcome variable (e.g. binary, count data, ordered choice, multinomial choice, \ldots) there exist appropriate specifications for $f(y_i \, | \,z_i,a_i;\beta)$ that allow point identification and $\sqrt{n}$-consistent estimation of $\beta$ by the conditional likelihood method. In static models, this approach relies on the availability of a sufficient statistic for $A_i$ (conditional on $Z_i$), which is satisfied in exponential-family models.\footnote{To provide a non-exhaustive list of examples, see, e.g., Rasch(61), Andersen(70), Chamberlain(80), Chamberlain(85) for binary choice logit; lancaster2000incidental, BluGriWin2002 for count data Poisson; and das_panel_1999, Baetschmann2015, Muris2017 for ordered choice logit models (using binarization).} In dynamic panel models, one can similarly find specifications for $f(y_i \, | \,z_i,a_i;\beta)$ such that estimation of $\beta$ via the generalized method of moments is possible.\footnote{See, for example, honore2020dynamic, kitazawa2021transformations for dynamic binary choice logit; BluGriWin2002 for dynamic count data Poisson; and honore2021dynamicOrdered for dynamic ordered choice. HonoreKyriazidou(00) and bartolucci2010dynamic also consider estimation of $\beta$ in dynamic binary choice panel models.} More generally, the functional differencing method of bonhomme2012functional can be viewed as a unifying framework for point-estimation of $\beta$ in both static and dynamic panel models of the form (ref).
Notice also that there are interesting models that do not require estimation of any common parameters $\beta$. A prominent example is the binary choice random coefficient model, which allows for richer forms of heterogeneity than the classical fixed effects specification; see Example (ref) below and our discussion in Section (ref). An alternative approach to such models uses finite discrete mixtures BC07, BC10, BC14, for which BC13 establish identification conditions in terms of the number of time periods and mixture components. Our framework accommodates both continuous and discrete specifications for the distribution of unobserved heterogeneity.
The rest of the paper is organized as follows: The main idea of our approach is introduced in Section (ref). Section (ref) presents the general construction of our bounds, including the linear programs used to obtain them. Section (ref) provides further discussion of the bounds, including an illustrative example and comparisons to the sharp identified set. Section (ref) addresses inference when common parameters must be estimated, providing two approaches for constructing asymptotically valid confidence intervals. Sections (ref) and (ref) present simulation evidence and an empirical application to female labor force participation, respectively. Section (ref) concludes. Proofs and additional results are provided in the Supplementary Appendix.
We observe discrete outcomes $Y_i \in {\cal Y}$, and conditioning variables $Z_i \in {\cal Z}$ for a cross-sectional sample of units $i=1,\ldots,n$. Unobserved heterogeneity is modeled through an unobserved latent variable $A_i \in {\cal A}$. The probability of observing $Y_i = y$ conditional on $Z_i=z$ and $A_i = a$ is given by $f\left(y\, |\, z,a; \beta_0 \right)$ where $\beta_0 \in \mathcal{B} \subset \mathbb{R}^{\dim \beta}$ and $f : \mathcal{Y} \times \mathcal{Z} \times \mathcal{A} \times \mathcal{B} \rightarrow [0,1]$ is a known function. The joint distribution of the conditioning variables $Z_i$ and $A_i$ is left unspecified. We focus on panel data models, where $Y_i = (Y_{i1},\ldots,Y_{iT})$ is a vector of outcomes $Y_{it} \in {\cal Y}_t$. The vector of conditioning variables $Z_i$ can, for example, be equal to $X_i = (X_{i1},...,X_{iT})$ in static models, or to $Z_i=(X_i,Y_{i0})$ in dynamic models where $Y_{i0}$ is the initial condition from time period $t=0$. We assume throughout that the covariates $X_i$ are strictly exogenous, meaning that $(X_{i1},\ldots,X_{iT})$ is independent of $(\varepsilon_{i1},\ldots,\varepsilon_{iT})$ conditional on $A_i$.\footnote{This rules out predetermined (but not strictly exogenous) covariates such as lagged values of other endogenous variables. However, lagged values of the dependent variable $Y_{it}$ itself can be accommodated, as in our dynamic model examples, because these enter through the conditioning set $Z_i$ rather than as covariates $X_i$.} In dynamic models, we assume that the initial condition $Y_{i0}$ is observed. Our goal is to provide inference methods on average effects of the form
where $m : \mathcal{Z} \times \mathcal{A} \times \mathcal{B} \rightarrow \mathbb{R}$ is a known function.
To focus on the main features and intuition behind our proposed approach, in this section we abstract away from estimation of $\beta_0$ and assume that it is known. In Section (ref) we will consider the case where $\beta_0$ is unknown but point-identified. The case of set-identified $\beta$ (along with a simulation analysis for the probit model) is considered in Sections (ref) and (ref) in the Supplementary Appendix. The random coefficient model in Example (ref) below provides an interesting case where no estimation of $\beta_0$ is necessary, because the model does not feature any such common parameter. In that case, the results in this section are already fully sufficient for inference on $\overline m$.
While our approach is general enough to accommodate different panel models of interest (including dynamic ones), for illustration purposes we focus on two running examples.
Our proposal for inference on $\overline m$ is based on the simple idea that suitable non-random functions $L, U : {\cal Z} \times {\cal Y} \times {\cal B} \rightarrow [b_{\min}, b_{\max}]$ which satisfy,
can be used to obtain asymptotically valid bounds on $\overline{m}$. To see how, notice that when evaluated at $\beta_0$, the condition in (ref) is equivalent to
which, by the Law of Iterated Expectations, implies that
This suggests that asymptotically valid bounds on $\overline m$ are given by
To formally show this, we impose the following regularity conditions.
Assumption (ref)(i) demands cross-sectional sampling. Assumption (ref)(ii) imposes correct specification of our parametric model for $Y_i$ conditional on $Z_i$ and $A_i$. This assumption also implies that all covariates contained in $Z_i$ are strictly exogenous (as opposed to pre-determined), as mentioned before. Assumption (ref)(iii) requires uniform bounds on the functions $m\left(z,a,\beta \right)$ that define the average effect of interest $\overline m$. This holds for typical choices for $\overline{m}$ such as those in Examples (ref) and (ref), and it can easily be confirmed for any given $m(z,a,\beta)$.\footnote{In principle a weaker condition such as $b_{\min} \leq \mathbb E[m\left(Z_i,A_i,\beta_0 \right) |A_i=a] \leq b_{\max}$ might also be used here, or bounds on second or higher-order moments of $m\left(Z_i,A_i,\beta_0 \right)$ are also conceivable, but in all the applications we consider in the paper the original Assumption (ref)(iii) holds, and we find it attractive that this assumption can be verified without knowing anything about the data generating process of $Z_i$ and $A_i$. More generally, Assumption (ref)(iii) could be replaced by any assumption that guarantees that ${\rm Var}\left[L(Z_i,Y_i, \beta_0 )\right]$, and ${\rm Var}\left[U(Z_i,Y_i, \beta_0 )\right]$ are finite in Theorem (ref). } Importantly, we do not put any restriction on the joint distribution of $Z_i$ and $A_i$. In particular, $Z_i$ can be discrete or continuous, $Z_i$ and $A_i$ can be arbitrarily related, and they do not have to be bounded or have bounded moments.
We now introduce our general construction of the bound functions $L(z,y,\beta )$ and $U(z,y,\beta )$. To concentrate solely on bound construction, in this section we still consider the case with known $\beta_0$. A full theory with estimated $\beta_0$ is provided in Section (ref). In terms of implementation, the construction methods remain the same for given $\beta$, independent of whether it is $\beta_0$ or its estimate.
In obtaining asymptotically valid bounds, the key requirement on the functions $L(z,y,\beta )$ and $U(z,y,\beta )$ is that they satisfy (ref) and that they are bounded. Of course, one wants the estimated bounds on $\overline m$ to be informative, in the sense that the interval in (ref) is as narrow as possible. At the same time, importantly, for given $z$ and $\beta$, $L(z,y,\beta )$ and $U(z,y,\beta )$ have to be chosen such that (ref) holds for all $a\in \mathcal A$. This can be reformulated as a standard optimization problem. Namely, for any given $z \in {\cal Z}$ and $\beta \in {\cal B}$ we can choose $L(z,y,\beta)=\ell(y)$ and $U(z,y,\beta)=u(y)$ as solutions to the following optimization problem with some appropriate objective function $Q(\ell(\cdot),u(\cdot),z,\beta)$,\footnote{ The solutions $L(z,y,\beta)=\ell(y)$ and $U(z,y,\beta)=u(y)$ may not be unique. But in a practical implementation some concrete solution will still be obtained by the specific linear solver used for implementation, and Theorem (ref) is still valid, since it only depends on the constraints being satisfied. }
In the current setting where we assume that $\beta_0$ is known, (ref) will be solved at $\beta=\beta_0$. When $\beta_0$ is estimated, (ref) will be solved at some estimate $\beta=\widehat{\beta}$. When no common parameter is estimated (as in Example (ref)), the objective function and the constraints will be free of $\beta$.
The restrictions of the program (ref) guarantee the conditions of Theorem (ref) and also impose that $\ell(y) \leq u(y)$. Consequently, any choice of the objective function $Q(\ell(\cdot),u(\cdot),z,\beta)$ yields valid bounds with $\widehat L \leq \widehat U$. It is important to stress that in order to construct the bounds $\widehat L$ and $\widehat U$ we only need to solve the program in (ref) once for every $i\in \{1,...,n\}$ at $z=Z_i$. Contrary to the sharp identified set, construction of our bounds does not involve conditional choice probabilities, and therefore remains free of the curse of dimensionality.
Display (ref) states our approach to obtaining bounds in its most general form, in the sense that the econometrician can choose any objective function $Q(\ell(\cdot),u(\cdot),z,\beta)$ that she sees fit. It is computationally attractive to consider objective functions which turn the optimization problem into a linear program, and we now discuss two intuitive choices of objective functions that are indeed linear in $\ell(\cdot)$ and $u(\cdot)$.\footnote{ In particular cases, it might be possible to obtain analytic expressions for the bound functions. But for the class of semi-parametric panel models and average effects introduced in Section (ref), it is unlikely that analytic expressions for the bounds can be obtained in general. The distinction between a numerical method and analytic expressions for the bound functions is analogous to the distinction between the functional differencing method in bonhomme2012functional and the analytical moment functions in honore2020dynamic for the purpose of inference on $\beta$. }
A linear program can be implemented by using the objective function
where $p(a|z)$ is some (potentially non-proper) “prior distribution”. In the absence of any additional information on $a$, such as the case considered in this paper, one can simply use $p(a|z)=1$. Indeed, throughout all the applications based on this baseline linear program, we use the flat prior $p(a|z)=1$. We note, however, that our bounds are valid for any choice of “prior”. When available, any extra information on the distribution of unobserved heterogeneity can be incorporated into the choice of $p(a|z)$, but we do not pursue this here.
If we are unwilling to specify a prior $p(a|z)$, then we can choose the objective function
where, instead of integrating over $a\in \mathcal{A}$ with a prior distribution, we choose the worst-case value of $a\in \mathcal{A}$ that maximizes the expected bounds $\sum_{y\in \mathcal{Y}}\left[ u(y)-\ell(y) \right] \,f(y \, | \, z,a;\beta )$. Hence, we call the ensuing approach the uniform linear program. To be precise, this objective function cannot be used directly to yield a linear program since it is not linear in $u(y)$ and $\ell(y) $; however, an equivalent representation of this problem as a linear program is obtained as follows:
In this linear program, the variable set is extended by $s \in \mathbb{R}$. When profiling out $s \in \mathbb{R}$ from this program one finds that for given $\ell,u \, : \, \mathcal{Y} \rightarrow \mathbb{R}$ the optimal $s$ is given by
which is identical to the objective function in (ref). Thus, solving the linear program in (ref) gives the desired bound functions $L(z,y,\beta)=\ell(y)$ and $U(z,y,\beta)=u(y)$ that correspond to choosing the objective function (ref) in our general program (ref).
Both the baseline and the uniform linear programs are valid options and will yield valid outer bounds. The uniform program focuses on the worst-case value of $a\in \mathcal A$ and is therefore likely to yield (slightly) wider outer bounds compared to the baseline linear program.\footnote{In preliminary analysis available upon request, we have considered both linear programs for the four models analyzed in the simulation study of Section (ref). Our results reveal for the static logit, dynamic logit and the random coefficient static logit models that while the uniform linear program led to some widening of the bounds, the change was not substantial. The only significant change was observed for the random coefficient dynamic logit model.} Nevertheless, our general recommendation is to use the uniform linear program, as it does not require specifying a prior $p(a|z)$ and therefore requires less input from the researcher.
In practice, implementation of the linear programs (ref) and (ref) requires choosing a grid $\mathcal{A}_g \subset \mathcal{A}$ to approximate the constraints. For logit-based models, computational efficiency can be improved by rewriting the constraints in terms of sufficient statistics. These and other implementational details are discussed in Section (ref) of the Supplementary Appendix.
The following example simply corresponds to the nonparametric bounds in CFHN13. It is therefore not representative of how we obtain the bounds in this paper in general, but we still find the example instructive, since it provides analytical expressions for bounds satisfying (ref).
We consider the static binary choice model of Example (ref) for the case where $X_{it}\in\{0,1\}$ is the only covariate and the error term $\varepsilon_{it}$ is stationary over time $t$. The average effect is given by (ref) with $x_1=1$ and $x_2=0$, that is,
where the time averaging is not needed due to stationarity.\footnote{In this case, $Z_i=X_i=(X_{i1},\ldots,X_{iT})$. Notice that, contrary to the general case, here $m(A_i,\beta_0)$ does not depend on $X_i$. This is because the average effect is calculated with respect to specific values of $X_{it}$ and there are no other covariates.} For $d \in \{0,1\}$, let
For $v(X_i,d)=1$ we define $\overline{Y}(Y_i,X_i,d) := \sum_{t \in {\cal T}(X_i,d)} Y_{it} / \left| {\cal T}(X_i,d) \right|$ to be the average of $Y_{it}$ over those time periods ${\cal T}(X_i,d)=\{ t \, : \, X_{it}=d\}$ where $X_{it}$ equals $d$. For $v(X_i,d)=0$ we simply let $\overline{Y}(Y_i,X_i,d) :=0$.\footnote{ Essentially, $\overline{Y}(Y_i,X_i,d)$ can be defined as any real number, given that its contribution to the bounds will be equal to zero whenever $v(X_i,d)=0$.} Valid outer bound functions are then given by
The stationarity assumption then guarantees that
which is exactly the condition (ref) that our bound functions are supposed to satisfy.\footnote{ Due to stationarity we have
while for $v(X_i,d)=0$, the above bounds $L(X_i,Y_i)$ and $U(X_i,Y_i)$ simply revert to the appropriate worst-case bounds (zero or one) that are possible for the unidentified expectations.}
Again, we want to point out that this example is not characteristic of our bounds more generally. In particular, here $L(X_i,Y_i)$ and $U(X_i,Y_i)$ do not depend on $\beta_0$, and neither the single-index structure $X_{it}\beta_0 + A_i + \varepsilon_{it}$ nor the parametric assumption on the error distribution are utilized to show validity of the bounds --- the bounds here are valid for any model $Y_{it}=g(X_{it},A_i,\varepsilon_{it})$, as long as the function $g(\cdot,\cdot,\cdot)$ is constant over $t$, and the conditional distribution of the shocks $\varepsilon_{it}$ is stationary over $t$.
From the corresponding discussion in CFHN13 we also know that, as $T \rightarrow \infty$, the width of these bounds, $ \mathbb{E}[U(X_i,Y_i) - L(X_i,Y_i)]$, shrinks proportionally to the probability of $X_{it}$ being constant over $t$. Under appropriate distributional assumptions on $X_{it}$ (e.g.\ $X_{it}$ independent across $t$ and random), this implies that the width of the bounds shrinks exponentially fast in $T$. While we do not explore large-$T$ analysis here, we suspect that similar results hold more generally for the bounds in this paper.
The key difference between our outer bounds, $ \mathbb{E} U(X_i,Y_i)$ and $\mathbb{E} L(X_i,Y_i)$, and the identified set for $\overline m$ is how they depend on the conditional choice probabilities, $P(Y_i|X_i)$. In particular, while our outer bounds are linear functions of choice probabilities, the upper and lower boundaries of the identified set are complicated nonlinear functions of $P(Y_i|X_i)$. The goal of this subsection is to briefly explain this difference and its consequences for inference on the average effects.
For simplicity, we stick to the static binary choice example with a single binary covariate discussed in the last subsection, and we assume that $\varepsilon_{it}$ has standard logistic distribution. Let $f(y|x,a;\beta_0)$ be the corresponding conditional distribution of $Y_i|X_i,A_i$. As long as we have some variation on the covariates across time, $\beta_0$ is point-identified in this model (see e.g.\ Chamberlain(85),Chamberlain(10)).
For $x \in \{0,1\}^T$, let $p(x):= \big[P(Y_i=y \, |\, X_i=x) \,: \, y \in \{0,1\}^T \big] $ be the $2^{T}$-vector of choice probabilities conditional on $X_i=x$, and define $\overline m(x) := \mathbb{E}\left[ m(A_i,\beta_0) \, \big| \, X_i=x \right]$. Next, let $\Pi(x,p(x))$ be the set of conditional distributions $A_i|X_i$ that are compatible with the choice probabilities $p(x)$: that is, we have $\pi(\cdot|x) \in \Pi(x,p(x))$ if and only if $P(Y_i=y \, |\, X_i=x)=\int_{\mathbb{R}} f(y|x,a;\beta_0) \pi(a|x) da$. Since $\beta_0$ and $p(x)$ are point-identified, the only ambiguity in the identification of $\overline m(x) $ is due to the unknown distribution of $A_i|X_i$. Then, defining
the identified set for $\overline m = \mathbb{E}\left[ \overline m(X_i) \right]$ is given by $\big[ \mathbb{E} L_{\rm id}(X_i,p(X_i)), \, \mathbb{E} U_{\rm id}(X_i,p(X_i))\big]$. All this is of course well-known. What we want to highlight here is that the above construction inevitably yields a complicated nonlinear dependence of the boundaries of the identified set on the observable choice probabilities $p(x)$ through $\Pi(x,p(x))$. In contrast, our bounds
are by construction linear functions of the vector of conditional choice probabilities $p(x)$.
This distinction between non-linearity (for the identified set) vs linearity (for our outer bounds) in $p(x)$ has a fundamental effect on inference: the sample analogs of our bounds, $\frac 1 n \sum_{i=1}^n L(X_i,Y_i)$ and $\frac 1 n \sum_{i=1}^n U(X_i,Y_i)$, avoid estimating $p(x)$ naturally. In contrast, we are not aware of any inference procedure on the sharp identified set that would avoid consistent estimation of $p(x)$.\footnote{ DDL21 present two different inference procedures for average effects in static panel logit models, one that relies on consistent estimation of $p(x)$, and one that does not. In the latter case, they also obtain certain outer bounds on the identified set that are different from our proposal. In Section (ref) in the Supplementary Appendix, we consider a brief comparison between their outer bounds and the ones proposed in this paper. } Especially when $p(x)$ is hard to estimate, the nonlinear dependence of the identified set on $p(x)$ can cause significant issues in inference. Hence, as already mentioned in the introduction, reliable inference on the identified set is problematic unless the sample size $n$ is much larger than the number of possible values for $(X_i,Y_i)$. Our outer bounds are by design immune to this.
To illustrate the points made here, we consider a brief simulation exercise. Let
where $\varepsilon_{it}\sim \mathrm{Logit}(0,1)$ and $x_{it}$ is discrete uniform with support $[0,\mathcal X -1]$. Then, $X_{it}$ can take on one of $|\mathcal X|$ equidistant values between 0 and 1. We consider $|\mathcal X|\in\{6,12\}$. The analysis for either case is based on 1000 replications of panels with $T=2$, $n=200$. The average effect of interest is as in (ref). For each replication, we obtain the estimated sharp identified set and our outer bounds based on the construction in Section (ref).\footnote{ The outer bounds presented here, based on the construction in Section (ref), provide much narrower bounds than the simple analytical expressions in (ref)-(ref). This is not surprising, given that our bounds utilize stronger model assumptions.} Then we report the 2.5% and 97.5% sample quantiles of these quantities across all replications.\footnote{More precisely, the reported $2.5\%$ sample percentile for the estimated identified sets corresponds to the $2.5\%$ sample percentile of the estimated lower bounds (of identified sets) across all replications. Similarly, the reported $97.5\%$ sample percentile corresponds to the $97.5\%$ sample percentile of the estimated upper bounds (of identified sets) across all replications. The sample percentiles for the outer bounds are obtained analogously. See also Section (ref) in the Supplementary Appendix for more information on the calculation of the identified set results.} Doing so enables us to compare the limits of the estimated confidence intervals, without estimating the confidence bands directly. Results are presented in Figure (ref). When $|\mathcal X | = 6$, the lower and upper $2.5\%$ percentiles of the estimated bounds of the identified set provide valid coverage. However, when $|\mathcal X|$ increases to 12, the same percentiles fail to include the average effect itself almost all the time. This reflects an underlying bias in the estimation of the sharp identified set. The outer bounds are immune to this issue, and still provide valid coverage. This example illustrates that, although the outer bounds are not sharp, they can be more reliable in inference compared to estimators of the sharp identified set itself. The results suggest, as expected, that issues arise as the cardinality of the support of the covariate increases. Therefore, the case with continuous $X_{it}$ will be subject to more pronounced issues.
We now extend this comparison to the population level for several specific logit-based binary choice models: the static logit and random coefficient logit models, and dynamic variants thereof. In all cases (except for the random coefficient dynamic logit model) we use the linear program in (ref). Here we compare our outer bounds to the population sharp identified set, estimation of which is challenging whenever the support of the conditioning variables $(Z_1,\ldots,Z_T)$ is not small relative to the sample size.\footnote{Similar to CFHN13 we obtain the sharp identified set by solving an appropriate linear program. See Section (ref) in the Supplementary Appendix for a more detailed discussion.} Our simulation results below present the identified sets and our outer bounds. In Section (ref) in the Supplementary Appendix, we also provide results concerning the widths of the outer bounds and the identified sets.
For static logit we consider both the discrete and continuous covariate cases, with the data generating processes (DGPs) given by
and
respectively. In both cases, $\varepsilon_{it} \sim {\rm Logit}(0,1)$. For the discrete covariate case we consider the average effect based on (ref) with $(x_1,x_2)=(1,0)$. The analysis for the continuous covariate case focuses on the average effect based on (ref). To focus solely on the difference between the bounds and the identified set, we set $\beta=\beta_0$ (a simulation analysis for obtaining bounds when $\beta$ is estimated will be provided in Section (ref)).
Results are presented in Figure (ref), where the reported outer bounds are the averages of the estimated bounds across 1000 replications of panels with $n=1000$.\footnote{ Since we are averaging over a large number of replications, the bounds reported in Figure (ref) are essentially equal to the population outer bounds $ \mathbb{E}[L(Z_i,Y_i,\beta_0)]$ and $ \mathbb{E}[U(Z_i,Y_i,\beta_0)]$, which justifies the comparison to the identified set. The same comment applies to the comparisons made in Figures (ref)-(ref). } The identified set and the outer bounds are obtained for $\beta_0 \in [-2,2]$. The support of $A_i$ is approximated by a grid of 100 equidistant points between $-5$ and 5. Several observations are in order. First, in all cases, the outer bounds mimic the behavior of the identified set. In particular, both the identified set and our bounds shrink to a point when $\beta_0=0$ but become wider as $|\beta_0|$ increases. At $T=5$ both the bounds and the identified set become almost a point for the majority of $\beta_0$ we consider. Also, both types of bounds yield the correct sign for the average effect. Second, the difference between the identified set and the outer bounds vanishes almost completely at $T=5$. This is an important result: as mentioned previously, obtaining the identified set in applications with moderate $T$ is practically infeasible due to the large number of conditional probabilities $P(Y=y|Z=z)$ one has to estimate, even when $Z$ is discrete. Our results show that the method proposed here stands out as a viable and computationally feasible alternative in such cases.
The random coefficient example is based on the DGP
where $\varepsilon_{it} \sim {\rm Logit } (0,1)$. Our interest is in identifying the average effect based on (ref). We note that Theorem (ref) fully applies here, as there are no structural parameters to be estimated.
Results are based on 1000 replications, and are presented in Figure (ref). We consider $n=1000$ and $T\in \{ 3,5,8,10 \}$ with $A_2 \in [-2,2]$. To approximate the supports of $A_{i,1}$ and $A_{2,i}$ we use grids of 50 equidistant points between $-5/5$ and $-7/7$, respectively. Not surprisingly, the presence of a random coefficient renders the average effect more difficult to identify. Indeed, for small $T$ even the sign of the average effect remains inconclusive for values of $A_2$ close to zero. More importantly, although the identified set becomes narrower as $T$ increases, it does not shrink to a point even when $T$ is 8 or 10. For reasons discussed before, obtaining the identified set for such large $T$ will in practice be infeasible. Simulation results confirm that our proposed method provides a reliable alternative. Indeed, the outer bounds are quite close to the identified set at $T=8,10$.
We next focus on the dynamic logit model with a continuous covariate. The DGP is
and we consider the average effect
where
AC21 have shown that in the given setting, the average effect in (ref) is point-identified (see their Proposition 3). The comparison in this part is then that between the outer bounds and the point-identified average effect. We investigate the behavior of the outer bounds in this case in panels of size $n=1000$ and $T\in \{ 4,6,8 \}$ with $\beta=1$ and $\gamma \in \left[ -2,2 \right]$. The support of $A_i$ is approximated by a grid of 50 equidistant points between $-5$ and $5$. The results, presented in Figure (ref), are based on 1000 replications and confirm that the outer bounds nearly point-identify the average effect, unless when $\gamma$ is large; however this issue tends to disappear as $T$ increases. This is not surprising since under a large $\gamma$, the term $\gamma Y_{i,t-1}$ will act similar to a fixed effect.
Finally, we consider the random coefficient dynamic logit model given by
For this exercise, we focus on the average effect
We consider 1000 replications where $n=1000$ and $T\in \left\{ 4,6,8,10 \right\}$, and vary $A_2$ between $-2$ and $2$. As in the static logit variant of this model, the supports of $A_{i,1}$ and $A_{2,i}$ are approximated by grids of 50 equidistant points between $-5/5$ and $-7/7$, respectively. Figure (ref) reveals that the identified set can be quite wide. This is in line with the earlier observations for the random coefficient static logit model. However, the identified set becomes wider as $A_2$ increases. This is similar to the asymmetry observed in the dynamic logit case. When $A_2$ is large, $Y_i$ is more likely to be a vector of 1s. Hence, again, the effect of $A_{i,2}Y_{i,t-1}$ is hard to distinguish from that of $A_{i,1}$. In results not reported here, we observed that the outer bounds tend to be conservative in this particular case when the uniform linear program (ref) is used. We therefore used the baseline linear program which utilizes (ref) in obtaining the bounds reported in Figure (ref). The resulting outer bounds perform well in tracking the identified set as $T$ increases.
We now consider the case where the common parameter vector $\beta_0$ has to be estimated. Our construction of the bound functions $L\left(z,y,\beta \right)$ and $U\left(z,y,\beta \right)$ remains essentially unchanged, but they are now evaluated at a consistent estimator of $\beta_0$, rather than the true $\beta_0$. The goal here is to provide asymptotic results that account for the noise in the estimation of $\beta_0$.
If the bound functions $L\left(z,y,\beta \right)$ and $U\left(z,y,\beta \right)$ were differentiable in $\beta$, then accounting for the estimation of $\beta_0$ when providing one-sided confidence intervals on the bounds $\mathbb{E}[L(Z_i,Y_i,\beta_0)]$ and $\mathbb{E}[U(Z_i,Y_i,\beta_0)]$ would be a straightforward application of the delta method. Unfortunately, because we obtain $L\left(z,y,\beta \right)$ and $U\left(z,y,\beta \right)$ as the solution to a linear program, it is generally not possible to verify any smoothness of those functions in $\beta$.\footnote{For example, in the static logit model, we know that the upper and lower bound functions are generally not unique for any given $\beta$, implying that they also cannot be continuous or smooth as functions of $\beta$. This shows that smoothness of $f(y|z,a,\beta)$ and $m(z,a,\beta)$ in $\beta$ does not imply smoothness of the bounds in $\beta$. } The convergence rate and inference results in this section therefore make no assumption whatsoever on the continuity or smoothness of the bound functions.\footnote{One could, alternatively, construct $L\left(z,y,\beta \right)$ and $U\left(z,y,\beta \right)$ such that they still satisfy the assumptions of Theorem (ref), but are also smooth in $\beta$ (e.g.\ in a particular model for a particular average effect of interest, one may simply find explicit analytic expressions for the bound functions). We leave the exploration of such possibilities to future work.}
Before discussing inference on $\overline m$, our first goal is to show that the population bounds $\mathbb{E}[L(Z_i,Y_i,\beta_0)]$ and $\mathbb{E}[U(Z_i,Y_i,\beta_0)]$ can be estimated at $\sqrt{n}$ rate, even if $\beta_0$ is estimated. For that purpose, we split the set of observations $\{1,\ldots,n\}$ into the disjoint subsets ${\cal I}_1 = \{1,\ldots,\lfloor n / 2 \rfloor \}$ and ${\cal I}_2 = \{ \lfloor n / 2 \rfloor +1,\ldots,n\}$. For any subset of observed units ${\cal I} \subset \{1,\ldots,n\}$ we denote by $Y_{({\cal I})}$ and $Z_{({\cal I})}$ the collection of all observations $Y_i$ and $Z_i$ with $i \in {\cal I}$. Furthermore, we define the function $\bar s \,: \, \{ 1,\ldots, n\} \rightarrow \{1,2\}$ by
For each $s \in \{1,2\}$ we have an estimator $\widehat \beta_s = \widehat \beta_s(Y_{({\cal I}_s)},Z_{({\cal I}_s)}) $ that only depends on the observed data $(Y_i,Z_i)$ for $i \in {\cal I}_s$. In other words, $\widehat \beta_{1}$ and $\widehat \beta_{2}$ are estimators of $\beta$ obtained using the first and second half-sample, respectively. Our estimates for the upper and lower bounds in (ref) then generalize to
Notice that the “cross-fitting” construction in (ref) ensures that for any $i$, $(Z_i,Y_i)$ and $\widehat \beta_{\bar s (i)}$ are always from two different half-samples, and therefore independent of each other. Consequently, conditional on the half-sample $\mathcal I_{\bar s(i)}$, $L(Z_i,Y_i,\widehat \beta_{\bar s (i)})$ and $U(Z_i,Y_i,\widehat \beta_{\bar s (i)})$ are independently distributed over $i$. In contrast, if the bound estimators were based on $\widehat \beta$ obtained from the full-sample, $L(Z_i,Y_i,\widehat \beta)$ and $U(Z_i,Y_i,\widehat \beta)$ would be arbitrarily dependent over $i$, ruling out a standard Law of Large Numbers. Along with reasonable assumptions on the behavior of $\widehat \beta_{\bar s(i)}$, as well as smoothness conditions on the functions $f(y|z,a;\beta)$ and $m(z,a,\beta)$ in $\beta$, the conditional independence is sufficient for proving the consistency of the bounds in (ref) for $\mathbb{E}[L(Z_i,Y_i,\beta_0)]$ and $\mathbb{E}[U(Z_i,Y_i,\beta_0)]$.
This theorem generalizes consistency of the outer bounds to the case of estimated $\beta_0$. The proof is straightforward and provided in the appendix. By contrast, obtaining inference results under estimated $\beta_0$ is more complicated due to the linear program yielding potentially non-smooth bound functions. This non-smoothness is not specific to our case: it arises generically in methods that construct estimating equations or bounds via linear programming. For example, bonhomme2012functional constructs moment functions for panel data models by solving linear programs, and faces the same challenge that small changes in parameters can cause discrete jumps in the solution. Similarly, the bounds in CFHN13 are obtained via linear programming, and their “perturbed bootstrap” inference method is specifically designed to circumvent the non-smoothness problem. The bottom line is that one cannot simply deploy the delta method to account for randomness introduced by the estimation of $\beta_0$, and so a different approach is needed. In the remainder of this section, we introduce two inference methods.
Our first inference method is inspired by the handling of common parameters in the “perturbed bootstrap” approach of CFHN13. The idea is to simply take the union of our “known $\beta_0$” confidence intervals in Section (ref) over a confidence set of the unknown $\beta_0$. For that purpose, define
We then have the following theorem.
Theorem (ref) provides a straightforward albeit potentially conservative way of obtaining confidence bands that incorporate the uncertainty due to estimation of $\beta_0$. This uncertainty is captured by $\gamma$ whereas $\alpha$ parameterizes the uncertainty due to estimation of the population outer bounds by sample averages. For a desired level of confidence $1-c$, one can trade off between these two sources of uncertainty by choosing $\alpha$ and $\gamma$ as desired. Another option is to find the narrowest confidence interval across all $(\alpha,\gamma)$ such that $c=\alpha+\gamma$. Notice that the infimum and supremum cannot be calculated exactly, so one has to do a grid search across a sufficiently large selection of $\beta\in \mathcal B _{1-\gamma}$. Especially when $\beta$ contains several parameters, this method can be demanding. Nevertheless, the attraction of Theorem (ref) is that as long as a valid confidence interval for $\beta_0$ can be constructed, inference on $\overline m$ requires only a straightforward application of the methods described in Section (ref). Fortunately, there is a large literature on obtaining valid confidence intervals on the common parameters $\beta_0$ in the type of panel data models with fixed effects that we consider here; see, for example, arellano2003discrete and ArellanoBonhomme(11) for reviews, as well as our discussion in the introduction.
Interestingly, the confidence set for $\beta_0$ in Theorem (ref) can also accommodate cases where $\beta_0$ is not point-identified, as long as a valid confidence set ${\cal B}_{1-\gamma}$ can be constructed. This is particularly relevant for models such as the probit, where the common parameters are generally set-identified when $T$ is fixed DGK24. In Section (ref) of the Supplementary Appendix, we discuss two approaches for obtaining outer bounds when $\beta_0$ is set-identified: the first uses a random coefficients specification, while the second directly incorporates the identified set for $\beta_0$ into bound construction. Section (ref) provides simulation evidence for the probit model.
As mentioned before, evaluating the infimum and supremum over $\beta \in {\cal B}_{1-\gamma}$ in Theorem (ref) can be challenging. As an alternative inference method, we therefore suggest modifying the linear program that is used to calculate the upper and lower bounds for $\overline m$ such that the uncertainty about $\beta_0$ is accounted for within the constraints of the linear program.
In Section (ref), the crucial requirement on our bound functions $L\left( z, y, \beta \right)$ and $U\left( z, y, \beta \right)$ was that they satisfy the inequalities in (ref) for a fixed value $\beta$. To account for the fact that the true $\beta_0$ is unknown, we now slightly generalize this idea. Given a {\it finite} set ${\cal B}_{\rm{sub}} \subset {\cal B}$ of possible values for $\beta$, we demand that the bound functions $L\left( z, y, {\cal B}_{\rm{sub}} \right)$ and $U\left( z, y, {\cal B}_{\rm{sub}} \right)$ satisfy the inequalities in (ref) for each value $\beta \in {\cal B}_{\rm{sub}}$, that is, we demand
As in (ref), we want the inequality in (ref) to hold for all $z \in {\cal Z}$ and $a \in {\cal A}$.\footnote{For consistency of notation, our previous bounds $L\left( z, y, \beta \right)$ and $U\left( z, y, \beta \right)$ could have been written as $L\left( z, y, \{\beta\} \right)$ and $U\left( z, y, \{\beta\} \right)$ to agree with (ref), but this is a minor mismatch of notation. }
Next, for each half-sample $s \in \{1,2\}$, let $\widehat {\cal B}_{s}$ be a set of points estimated only from observations $i \in {\cal I}_s$, such that the convex hull of $\widehat {\cal B}_{s}$, ${\rm Conv}(\widehat {\cal B}_{s})$, provides a $1-\gamma/2$ confidence set for $\beta_0$. For example, for a one-dimensional parameter $\beta$, we choose $\widehat {\cal B}_{s}=\{\widehat \beta_{{\rm low},s},\widehat \beta_{{\rm up},s}\}$ to consist of the lower and upper bounds of a confidence interval for $\beta_0$. Then, ${\rm Conv}(\widehat {\cal B}_{s})= [\widehat \beta_{{\rm low},s},\widehat \beta_{{\rm up},s}]$ is just a standard confidence interval in that case. More generally, we have to find a confidence set that can be generated as a convex hull of a finite number of points.\footnote{ If $\beta$ is higher-dimensional, then one simple choice for $\widehat {\cal B}_{s}$ would be the Cartesian product of one-dimensional confidence bounds for each component of $\beta$, using a Bonferroni correction to maintain the correct confidence level $1-\gamma/2$. However, this yields a set with $2^{\dim(\beta)}$ vertices, which may be computationally costly. A more efficient construction uses the cross-polytope, which requires only $2 \dim(\beta)$ vertices. Suppose that $\sqrt{n/2}(\widehat\beta_s - \beta_0) \xrightarrow{d} N(0, \Sigma_\beta)$ with consistent estimator $\widehat\Sigma_{\beta,s}$. Let $c_{1-\gamma/2}^{(d)}$ be the $(1-\gamma/2)$-quantile of $\|Z\|_1 = \sum_{j=1}^d |Z_j|$ where $Z \sim N(0, I_d)$. Then
where $e_j$ denotes the $j$-th standard basis vector in $\mathbb{R}^d$. The set $\widehat{\cal B}_{s}$ consists of $2d$ vertices whose convex hull provides an asymptotically valid $(1-\gamma/2)$ confidence set for $\beta_0$. The critical value $c_{1-\gamma/2}^{(d)}$ can be computed numerically; for example, $c_{0.975}^{(1)} \approx 2.24$, $c_{0.975}^{(2)} \approx 3.02$, and $c_{0.975}^{(3)} \approx 3.67$. } Let also ${\rm diam}\left( {\cal B}_{\rm{sub}} \right)$ be the diameter of the set ${\cal B}_{\rm{sub}} $. Finally, we define
We require the following additional assumptions for this inference method, which strengthen Assumption (ref)(ii) and also formalize the requirement that ${\rm Conv}(\widehat{ \mathcal B} _s)$ is a confidence interval.
Then, the following lemma shows that conditional on $ \beta_0 \in \widehat{\mathcal B}_{\bar s}$, $ L(Z_i,Y_i,\widehat{\mathcal B}_{\bar{s}})$ and $U(Z_i,Y_i,\widehat{\mathcal B}_{\bar{s}})$ provide valid bounds on $\overline m$ in expectation.
Once Lemma (ref) is obtained, then all that is left to do is to account for the sampling uncertainty when replacing the expected value over $ L(Z_i,Y_i,\widehat{\mathcal B}_{\bar{s}})$ and $ U(Z_i,Y_i,\widehat{\mathcal B}_{\bar{s}})$ by the sample averages in (ref), analogously to Theorem (ref).
Theorem (ref) demands that equation (ref) holds, but does not specify any explicit construction of the bound functions. Analogously, Theorem (ref) requires that equation (ref) hold, but again does not specify any explicit construction of the bounds. In order to actually construct the bounds we use the methods described earlier, but we replace the constraint (ref) by (ref). Specifically, the program in display (ref) then gets modified as follows: For any given $z \in {\cal Z}$ and any finite set ${\cal B}_{\rm{sub}} \subset {\cal B}$ with $\overline \beta = \left| {\cal B}_{\rm{sub}} \right|^{-1} \sum_{\beta \in {\cal B}_{\rm{sub}}} \beta$ we can choose $L(z,y,{\cal B}_{\rm{sub}})=\ell(y)$ and $U(z,y,{\cal B}_{\rm{sub}})=u(y)$ as solutions to the following optimization problem:\footnote{Here, the choice of evaluating the objective at $\overline\beta$ is for convenience. Alternative choices, such as $\sum_{\beta \in {\cal B}_{\rm{sub}}} Q(\ell, u, z, \beta)$, are equally valid since the constraints (not the objective function) guarantee validity of the bounds.}
By choosing the objective function $Q(\ell(\cdot),u(\cdot),z,\overline \beta)$ as in (ref) or (ref), we again have to solve a linear program to obtain the bounds.
In this part we investigate the small sample behavior of the proposed bounds and confidence bands. We focus on the static logit and random coefficient logit models. The setting largely follows Section (ref). In particular, we use the DGPs in (ref) and (ref)-(ref) with a single discrete covariate, and focus on the same average effects. The main difference is that we now estimate $\beta_0$ in the static logit model, and also provide confidence bands. The results, presented in Figures (ref)-(ref), provide the population average effect, and the cross-replication averages of estimated bounds and 95% confidence bands.
We first consider the static logit model. $\beta_0$ is estimated using the conditional likelihood method. For inference we use the two inference methods proposed in Sections (ref) and (ref). In either case, we consider 1000 replications of panels of size $n=5000$ and $T\in \{ 3,5,8 \}$. For $\mathcal{A}_g$, we use a grid of 100 equidistant points between $-5$ and $5$.
The results using the inference method of Section (ref) are based on $\gamma = 0.0001$, and $\mathcal B_{1-\gamma}$ is approximated by a grid of 5000 equidistant points on $\mathcal B_{1-\gamma}$. Outer bounds for this case are obtained by the uniform linear program of Section (ref). Results are presented in Figure (ref). For moderate $T$, which is the main focus of this study, both the bounds and the confidence bands are quite tight. Interestingly, this is despite the fact that the bounds are based on a uniform linear program. In all cases, the confidence bands yield the correct sign for the average effect. The coverage rates of confidence bands are, not surprisingly, conservative. This is expected in partially identified settings and is acceptable given that the bands remain informative.
We next consider the inference approach of Section (ref), the results of which are presented in Figure (ref). Confidence bands are based on $\alpha = \frac{2}{3} \times 0.05$ and $\gamma = \frac{1} {3} \times 0.05$.\footnote{The choice of $\alpha = 2\gamma $ is not crucial and was only imposed to compensate for the fact that the confidence interval ${\rm Conv}(\widehat {\mathcal B} _s)$ is subject to one Bonferroni split, whereas the interval for estimated outer bounds is subject to two.} Both the confidence bands and the outer bounds are based on the linear program defined in (ref). Relative to the inference method of Section (ref), there are two differences: first, the confidence bands are overall visibly closer to the estimated bounds, across all $T$. This is not surprising given that the inference method of Section (ref) is based on the infimum/supremum bands. Second, while the outer bounds improve with $T$, they are not as tight as the bounds produced by the linear program in (ref). This likely results from (ref) incorporating the uncertainty due to $\widehat{\beta}$ in outer bound estimation (as opposed to (ref) which incorporates the same in the inference stage).
We move to the random coefficient static logit example. Figure (ref) presents results based on 1000 replications of panels of size $n=1000$ and $T\in \{ 3,5,10 \}$. We construct $\mathcal{A}_g$ using 50 equidistant grid points between $-5$ and $5$ for $A_{1,i}$, and between $-7$ and $7$ for $A_{2,i}$, leading to 2,500 grid points in total. We note that the average effects and outer bounds are the same as in Section (ref), since no parameter estimation is involved in this setting. The new result is the confidence bands, which are based on Theorem (ref). On average the confidence bands are reasonably close to the outer bounds.
A general observation across the three simulation exercises is that while the confidence bands are generally informative, the coverage rates are conservative. Of course, since we are looking at coverage rates for the true value of the average effects (as opposed to the identified set), it is well-known that coverage will generally be conservative in partially identified settings (e.g., see ImbensManski04 and Stoye21).\footnote{{ImbensManski04 provide a generic method for obtaining narrower confidence bands on sets which can also be used here. In results not reported here we nevertheless observed that the decrease in the width of confidence bands was quite moderate, with little or no change in the actual coverage rates.}}
Finally, we note that the total width of the confidence bands can be decomposed into two components: (i) the gap between the outer bounds and the identified set, which measures the cost of using outer bounds rather than sharp identification, and (ii) the gap between the confidence bands and the outer bounds, which measures the cost of estimation uncertainty. This decomposition helps assess how much of the total uncertainty is due to our methodological choice versus statistical imprecision. We report these results in Section (ref) of the Supplementary Appendix.
We consider an empirical analysis of female labor force participation, using the National Longitudinal Survey of Youth (NLSY) 1979 dataset. Our sample consists of data on women who were married throughout the sample and who were not in active forces or going to school.\footnote{An individual is classified as “in the labor force” if her status was recorded as working, with job not at work or unemployed. Individuals are considered as not in the labor force if their recorded status was keeping house, unable to work or other.} Also, we only include individuals who were observed at all periods under consideration.
First, we consider a random coefficient logit specification\footnote{Labor force participation is often modeled with state dependence hyslop1999state. We focus on static models here for comparability with existing methods; see Section (ref) for simulation evidence on dynamic specifications.}:
where, for individual $i$ and at time $t$, $LFP_{it}$ is the labor force participation indicator whereas $kids3_{it}$ is a binary variable which equals one if the individual has at least one child below the age of three. This is almost identical to the example considered by CFHN13, except that they assume a homogeneous coefficient $\beta$ for all individuals. Our objective is to obtain a confidence interval on the average effect
Our sample period for this analysis covers all even years from 1986 to 1998, which yields data on 929 individuals over seven years. For comparison, we also report the average effects based on the fixed effects logit (FE logit) and probit (FE probit) models, as well as the linear fixed effects model. We note that all these alternatives impose homogeneity of $\beta_i$, and calculate the average effects using estimated $(\alpha_i,\beta)$. Hence, they provide a point-estimate for the average effect. We also use (i) the bias-corrected logit (BC logit) and probit (BC probit) methods, which analytically correct $\widehat{\beta}$ for the incidental parameter bias, following cruz2017bias, and (ii) the split-panel jackknife method (SPJ probit and SPJ logit) of DhaeneJochmans2015, which directly corrects average marginal effects rather than $\widehat{\beta}$. We note that none of these alternative methods are designed for short-$T$ samples where average effects are not necessarily point-identified. For all methods under consideration, we provide the 95% confidence intervals. For the outer bounds this is obtained by using the normal approximation of Theorem (ref).
\afterpage{
}
Results for this first illustration are reported in the top panel of Table (ref). All methods agree that having at least one child younger than three has a negative impact on labor force participation. This is also in line with the results obtained by CFHN13 who consider a shorter sample, covered by our dataset (see their Table III). The confidence intervals for the outer bounds are wider than the rest, but this is normal as it is the only method that allows for heterogeneity of $\beta_i$. Heterogeneity of $\beta_i$ is quite likely, as the effect of having a child younger than three will vary depending on various conditions. For example, families with higher income will have easier (and better) access to child care. Geographical proximity of grandparents (who can, at least from time to time, provide child care) is also likely to have an effect on $\beta_i$. Moreover, the effect of having children younger than three may differ depending on the actual number of children. The wider confidence bands provided by our method reflect all such considerations.\footnote{An alternative approach to accommodating heterogeneity in $\beta_i$ is to use finite discrete mixtures BC10, BC14. However, with $T=7$ periods, such models can identify at most approximately $T/2 \approx 3$ mixture components BC13, which may not fully capture the richness of heterogeneity in labor force participation responses.}
In the second illustration, we consider the static logit specification with a richer set of covariates:
where $educ_{it}$ is the highest completed grade (as of May 1 of the survey year) and $spouseinc_{it}$ is the total income of the spouse from wages and salary in past calendar year. The sample for this exercise covers all even years from 1990 to 1998. We do not include individuals whose spouse had zero income at any point during this period. Average effects for the covariates $kids3_{it}$ and $educ_{it}$ are based on (ref), where we use $(x_1,x_2)=(1,0)$ and $(x_1,x_2)=(educ_{it}+1,educ_{it})$, respectively. Average effects for log spouse income are calculated using $\eqref{marg2}$ with $x_{k,it}=\ln(spouseinc_{it})$. The outer bounds are obtained using the uniform linear program of Section (ref) whereas the inference approach of Section (ref) is used to generate the confidence bands.\footnote{We first obtain the confidence bands across a selection of $\alpha$ and $\gamma$ such that $\alpha+\gamma=0.05$, and then report the shortest confidence interval among these.}
Results are reported in the bottom panel of Table (ref). All methods agree that the average effect of $kids3$ is negative. For $educ$, all confidence bands are ambiguous about the size of the effect. However, for all methods these bands are mostly on the positive side. In addition, estimated average effects and outer bounds all point to a positive effect of $educ$ on labor force participation. We note that although almost all the point estimates for the average effects with respect to $kids3$ and $educ$ are outside the respective outer bounds, they are all covered by the confidence bands for the outer bounds. Finally, for log of spouse income, confidence bands by all alternatives (other than the linear model) are inconclusive about the sign of the average effect, though they mostly lie on the negative side. Interestingly, in this particular case the confidence bands for all methods other than the linear probability model lie partially outside the confidence intervals for the outer bounds. This is not necessarily surprising, given that none of the alternative methods considered here are designed to work in short samples.
In this paper, we have introduced a new method for estimating bounds on average effects in discrete choice panel data models with fixed effects, including two approaches for obtaining asymptotically valid confidence intervals on the average effects. For realistic models and sample sizes, inference based on our outer bounds is easier and more robust than inference based on the sharp identified set. A key strength of our approach is its broad applicability: it is suitable for models with both discrete and continuous covariates, and it can be adapted for a variety of static and dynamic panel models.
We have focused here primarily on the case where the common model parameters $\beta_0$ are point-identified and can be estimated at the parametric rate. In the Supplementary Appendix, we show how our approach extends to cases where the structural parameters are set-identified, with simulation evidence for the static probit model.
Another potential extension is to models with continuous outcomes $Y_{it}$, where the sums over $\mathcal{Y}$ in our linear programs would be replaced by integrals. In principle, such integrals could be approximated by sums over a discretized support, similar to how we approximate the constraints over $\mathcal{A}$ by a grid $\mathcal{A}_g$. We leave the exploration of this extension to future work.
\setstretch{1.39} {5pt plus 0.3ex}
\setcounter{page}{1}
\onehalfspacing