Partial Identification in Nonseparable Binary Response Models with Endogenous Regressors

Partial Identification in Nonseparable Binary Response Models with Endogenous Regressors We are grateful to James Heckman, Marc Henry, Roger Koenker, and to seminar audiences at Columbia University and Michigan State University for helpful feedback. We also thank Martin Weidner and the organizers of the Chamberlain Seminar, and are grateful to Florian Gunsilius, Sukjin Han, Wayne Gao, and Takuya Ura for their questions and feedback, and to Adam Rosen for his thoughtful discussion. Jiaying Gu acknowledges financial support from the Social Sciences and Humanities Research Council of Canada. All errors are our own.

Keywords: Binary Choice, Counterfactual Probabilities, Endogeneity, Hyperplane Arrangement, Linear Programming, Partial Identification

\thispagestyle{empty}

Introduction

This paper considers partial identification of counterfactual parameters in a class of binary response models of the form:

where $Y \in \{0,1\}$ is a binary outcome variable, $X \in \mathcal{X}\subset \mathbb{R}^{d_{x}}$ is vector of (possibly endogenous) covariates with finite support, $U \in \mathcal{U}:=\mathbb{R}^{d_{u}}$ is a vector of latent variables, $\theta \in \Theta \subset \mathbb{R}^{d_{\theta}}$ is a vector of structural parameters, and $\tilde{\varphi}_{1}:\mathcal{X}\times \Theta \to \mathbb{R}^{d_{u}}$ and $\tilde{\varphi}_{2}: \mathcal{X}\times \Theta \to \mathbb{R}$ are known functions. Our approach does not require any parametric distributional assumptions on the latent variables. We then consider counterfactuals that can be expressed as a function $\gamma: \mathcal{X} \to \mathcal{X}$ that reassigns each $x \in \mathcal{X}$ to a new value $\gamma(x) \in \mathcal{X}$. Our main focus is on bounding linear functionals of counterfactual probabilities.

In our setting, nonparametric point-identification of the distribution of latent variables occurs only under restrictive conditions, including strong independence assumptions and large support conditions.\footnote{E.g. ichimura1998maximum.} Control function approaches are often used to address the issue of endogenous regressors, but if endogenous regressors are discrete or the mechanism generating the endogenous regressors is poorly understood, then many of these approaches are not applicable. Partial identification arises as a natural alternative to methods for point-identification in the presence of endogenous and discrete covariates. This paper uses an optimization-based approach to bound counterfactual quantities, which allows researchers to easily construct sharp bounds on counterfactual quantities under a variety of different assumptions by simply altering the constraints in our optimization problems.

In the absence of parametric distributional assumptions, our analysis reveals the importance of a special partition of the latent variable space into response types that have identical responses in all possible counterfactuals. We are not the first to emphasize the importance of response types, and our discussion echoes the insights of balke1994counterfactual and heckman2018unordered, among others. Similar to these works, we show that functional form and monotonicity assumptions amount to assigning zero probability to certain response types. Before bounding counterfactual quantities, it is thus necessary to determine which response type are possible/impossible in a given model. In our particular class of models, the latent variable space admits a partition into cells defined by a collection of hyperplanes, each cell corresponding to a unique response type. Using the cell enumeration algorithm of gu2020nonparametric, we show how to enumerate all response types in a time polynomial in the number of input hyperplanes. After enumerating response types, we then demonstrate how various counterfactual quantities can be bounded by solving a sequence of linear programming problems. We also study the case when the index function is linear in parameters, in which case we show how sharp bounds on counterfactual probabilities can be computed without the need to grid over the entire parameter space. Finally, we present a consistency result, we show how to adapt the inference procedure of cho2021simple to our setting, and we apply our method to study the effects of private health insurance on the decision to seek medical treatment.

A Review of the Relevant Literature

Binary response models with endogenous regressors have been studied extensively. Point identified approaches include linear probability models, maximum likelihood estimation (e.g.\! the bivariate probit), control function approaches, or approaches based on special regressors. All of these approaches have well-documented limitations.\footnote{See lewbel2012comparing for a review. } Nonparametric identification was studied in binary choice and threshold crossing models by matzkin1992nonparametric, and in more general nonseparable models by matzkin2003nonparametric and chernozhukov2005iv, among many others. vytlacil2007dummy studied nonparametric identification of the average treatment effect in a discrete triangular system with a binary endogenous variable under a weak separability assumption in the outcome equation. Since then a number of paper have studied point identification in similar models with discrete endogenous variables (e.g. han2017identification, vuong2017counterfactual, chen2020identification and khan2021informational) and continuous endogenous variables (e.g. imbens2009identification, d2015identification, torgovitsky2015identification). Important precedents to the work presented here from the literature on point identification in random coefficient models include ichimura1998maximum, gautier2013nonparametric and gu2020nonparametric. However, these papers focus almost exclusively on the point-identified case with a linear index function and exogenous covariates with large support.

Many authors have also used partial identification methods to relax the assumptions required for point-identification in these models. In a relevant series of papers, chesher2013instrumental and chesher2014instrumental show how to use random set theory to characterize the identified set of structures in discrete choice models.\footnote{The general formulation of their approach is presented in chesher2017generalized.} Similar to the current paper, they do not provide a model for the endogenous explanatory variables, rendering the discrete choice model incomplete. They construct bounds on structural parameters using a characterization of the sharp set of constraints based on the results of artstein1983distributions.\footnote{See also norberg1992existence and molchanov2017theory Corollary 1.4.11.} Although our paper focuses primarily on computational issues that arise when bounding counterfactual parameters, we present a comparison of our approach with the approach of chesher2013instrumental and chesher2014instrumental in Appendix (ref).

Our approach to identification is closely related to approaches in galichon2011set, laffers2019identification and torgovitsky2019partial, who demonstrate how to construct sharp bounds on various parameters in models with discrete variables by partitioning the latent variable space and discretizing the latent variables. We also use an identification argument based on partitioning the latent variable space, and we demonstrate how to practically compute the relevant partition using results from the literature on computational geometry. In certain cases, we also show how to avoid griding over the entire parameter space when computing the identified set.

There are a number of other relevant papers in the literature on partial identification in discrete choice models. manski2007partial considers counterfactual choice probabilities in a setting with partial identification, and shows how these counterfactual choice probabilities can be bounded using optimization problems. However, the response type approach in our paper is quite different. We also show how to practically incorporate different assumptions, and we allow for endogenous explanatory variables.\footnote{This latter point differentiates our work from chiong2017counterfactual and allen2019identification.} In other related work, tebaldi2019nonparametric study the problem of computing various counterfactual quantities in a nonparametric discrete choice model with an application to consumer choice of health insurance in California. However, they focus on quasi-linear utility functions and use the particular structure of their setting to resolve the issue of endogeneity by conditioning on a set of covariates. Computational considerations are also not the main focus of these papers.

This paper is also related to papers that bound treatment effect parameters. Nonparametric bounds were proposed by manski1990nonparametric, and since then contributions have been made by manski1997monotone, manski2000monotone, heckman2001instrumental, bhattacharya2008treatment, manski2009more, chiburis2010semiparametric, shaikh2011partial, bhattacharya2012treatment and mourifie2015sharp, among many others. However, most papers derive closed-form bounds and provide corresponding proofs of sharpness for a fixed target parameter and under a fixed set of assumptions. Similar to other optimization-based approaches, the advantage of our procedure is its flexibility: the researcher can easily modify the model or target parameter and obtain new bounds without the need to derive closed form expressions, or to provide a corresponding proof of sharpness.\footnote{Examples of optimization-based approaches to bounding treatment effect parameters include balke1997bounds, laffers2019bounding, russell2021sharp, mogstad2018using and gunsilius2020path. } Our optimization-based bounds are shown to be sharp under a variety of different assumptions, including flexible functional form assumptions, and any mix of fixed and random coefficients on either exogenous or endogenous covariates. These assumptions are of direct interest to those familiar with structural models of binary outcome variables, and allow us to obtain many previous results as a special case. We hope our approach may also serve as a useful method of sensitivity analysis for those primarily interested in point identified models.

Finally, our paper makes numerous connections to the literature on computational geometry. Computation of our bounds requires the analysis of a partition of the latent space determined by a finite collection of hyperplanes. This turns out to be a well studied subject in combinatorial geometry, and leads us to consideration of the enumeration algorithm proposed by gu2020nonparametric, who build on the work of rada2018new. Our profiling procedure also makes use of the double-description algorithm proposed by Fukuda.

Paper Outline and Notation

The remainder of the paper proceeds as follows. Section (ref) introduces the main theoretical framework and main assumptions. Section (ref) studies practical implementation of the theoretical framework from Section (ref) and introduces our optimization-based bounding procedure for counterfactual probabilities. Section (ref) then demonstrates how to introduce functional form, independence, and monotonicity assumptions into our bounding procedure. Section (ref) discusses estimation and inference, and Section (ref) applies our methodology to study the impact of health insurance on utilization of health care services. Section (ref) concludes. All proofs can be found in Appendix (ref). Appendix (ref) provides some additional discussion of the results presented in the main text. Appendix (ref) compares our procedure to an approach based on Artstein's inequalities, and Appendix (ref) contains supplementary material for our application. \\ \%, and a comparison of our approach to the approach based on Artstein's inequalities from chesher2013instrumental and chesher2014instrumental is presented in Appendix (ref).\\ \\

Notation: The following notation is relevant for both the main text and the appendices. Given a subset $\mathcal{X}$ of Euclidean space, we use $\mathfrak{B}(\mathcal{X})$ to denote the Borel $\sigma-$algebra on $\mathcal{X}$. For two measurable spaces $(\mathcal{X},\mathfrak{B}(\mathcal{X}))$ and $(\mathcal{X}',\mathfrak{B}(\mathcal{X}'))$, the product $\sigma-$algebra on $\mathcal{X}\times \mathcal{X}'$ is denoted by $\mathfrak{B}(\mathcal{X})\otimes \mathfrak{B}(\mathcal{X}')$. Random variables are denoted using capital letters, and if $X: (\Omega,\mathfrak{A})\to (\mathcal{X},\mathfrak{B}(\mathcal{X}))$ is a random variable defined on the probability space $(\Omega,\mathfrak{A},P)$, then we use $P_{X}$ to denote the probability measure induced on $\mathcal{X}$ by $X$; that is, for any $A \in \mathfrak{B}(\mathcal{X})$, $P_{X}(A) := P(X^{-1}(A))$. Furthermore, we interpret $P_{X\mid X'}(X \in A \mid X' =x')$ as a regular conditional probability measure. Finally, $P_{X\mid X'}$ is used as shorthand for the collection $P_{X\mid X'}:=\{P_{X\mid X'}(\,\cdot\, \mid X'=x') : x' \in \mathcal{X}'\}$. We do not explicitly differentiate between scalars and vectors, or random variables and random vectors. To keep the notation clean, we sometimes omit the transpose when combining column vectors; that is, if $v_{1}$ and $v_{2}$ are two column vectors, rather than write $v=(v_{1}^\top,v_{2}^\top)^\top$ we instead write $v=(v_{1},v_{2})$, where it is understood that $v$ is a column vector unless otherwise specified. The cardinality of a set $\mathcal{X}\subset \mathbb{R}^{d}$ is given by $|\mathcal{X}|$.

General Framework: Theoretical Considerations

Main Assumptions and Definitions

We begin by introducing our main assumptions on the binary response models under consideration.

In Assumption (ref), $U \in \mathcal{U}$ is a vector of latent variables, $\theta \in \Theta$ is a vector of fixed coefficients, and $X\in \mathcal{X}$ is a vector of covariates. From (ref) we restrict the index function to be linear in the latent variables $U \in \mathcal{U}$, although the model in Assumption (ref) still allows for general nonseparability between covariates and latent variables. Importantly, Assumption (ref) imposes that the random vector $X$ has finite support.

In this model the latent variables can also be interpreted as random coefficients. A special case of linearity occurs when the function $\varphi$ is additively separable in a scalar latent variable $U$, which occurs, for instance, when $\tilde{\varphi}_{1}(x,\theta)=1$. A full analysis of this special case using the framework in this paper is taken up in Appendix (ref). Finally, assuming the distribution of $U$ assigns zero probability to sets of the form $\{u \in \mathcal{U} : \varphi(x,u,\theta)=0\}$ allows for a simplification of the cell enumeration algorithm introduced in the next section. This assumption is implied by absolute continuity of the distribution of latent variables with respect to the Lebesgue measure, which is a standard assumption in this literature.

We assume that the researcher's objective throughout is to obtain a sharp set of constraints defining the identified set of latent variable distributions, and to use these constraints to bound various counterfactual quantities, such as counterfactual conditional probabilities.\footnote{ Similar to previous works, we take the selection relation as a primitive relation on which to construct a definition of the identified set. The close connection between the selection relation from random set theory and the concept of observational equivalence from the work in econometrics on identification has been appreciated in beresteanu2011sharp, beresteanu2012partial, chesher2013instrumental, chesher2014instrumental, and chesher2017generalized, among many others. We continue this work here.} Define the set:

chesher2017generalized call this set the $U-$level set; intuitively, it delivers all possible values of the latent variables $u$ consistent with the vector $(y,x,\theta)$ given the binary response model in (ref). A measurable selection from the random set $ \mathcal{U}(Y,X,\theta)$ is a random vector $U: \Omega \to \mathcal{U}$ satisfying $U \in \mathcal{U}(Y,X,\theta)$ a.s.\footnote{A general definition of a selection and a random set is provided in Appendix (ref). In Appendix (ref) we prove that $\mathcal{U}(Y,X,\theta)$ is suitably measurable and thus is a random set under our assumptions (see Lemma (ref)). We also prove the existence of a universally measurable selection (see Lemma (ref)).} Importantly, given a distribution of the observable random vectors $(Y,X)$, a structural function $\varphi$ and a fixed coefficient $\theta \in \Theta$, any two measurable selections $U$ and $U'$ from the random set $\mathcal{U}(Y,X,\theta)$ are observationally equivalent in the sense that both latent variable vectors $U$ and $U'$ are consistent with the observed distribution of $Y$ and $X$ for the vector of parameters $\theta \in \Theta$ through the model (ref). This idea can be used to define the identified set.

This definition depends on the observed distribution of $(Y,X)$ through the almost-sure relation in (ref). Conditioning the latent variable distribution on the vector $(Y,X)$ is carried throughout the paper, and we show in Section (ref) that it allows us to bound counterfactual parameters that may be relevant to policy analysis. Definition (ref) can also be used to define other related identified sets, including identified sets for conditional latent variable distributions of the form $P_{U\mid X}$ or $P_{U}$.

In this paper, we consider counterfactuals characterized by the occurrence of an exogenous intervention that modifies at least one of the explanatory variables.

Assumption (ref) implies that (i) counterfactual response variables indexed by $\gamma \in \Gamma$ exist on the common probability space from Assumption (ref), and (ii) such counterfactual response variables are equal (almost surely) to the values that would arise after an intervention on the system represented by (ref). Here each counterfactual is represented by a function $\gamma: \mathcal{X} \to \mathcal{X}$, which allows us to consider a general class of counterfactuals. Although each function $\gamma$ is seen as a map from $\mathcal{X}$ to itself, this does not prevent consideration of counterfactuals where $\gamma$ selects values of $x$ that have never been observed in the data. Such cases can be accommodated by simply extending the support $\mathcal{X}$ from Assumption (ref) to include any counterfactual pair $x$ of interest.\footnote{In particular, if $x \in \mathcal{X}$ but $\gamma(x)=x' \notin \mathcal{X}$, then redefine $\mathcal{X} \leftarrow \mathcal{X} \cup \{x'\}$. This approach does not affect anything we present in this paper, since we always require any relation to the observed distribution of $(Y,X)$ to hold only almost-surely. }

Note this definition references the identified set $\mathcal{I}_{Y,X}^{*}$ from Definition (ref). This definition can also be used to define other related identified sets, including for counterfactual distributions $P_{Y_{\gamma}\mid Y,X}$ or $P_{Y_{\gamma}\mid X}$, as well as identified sets for average structural functions and average treatment effects.

\setcounter{example}{0}

The following result connects the definitions and assumptions in this section.

Theorem (ref) provides the theoretical link between the identified set of counterfactual conditional distributions, and the identified set for the pair $(P_{U\mid Y,X},\theta)$. While theoretically straightforward, it hides some important practical difficulties. In particular, verifying the existence of a pair $(P_{U\mid Y,X},\theta)$ that satisfies the conditions from Definition (ref) is a nontrivial task. This is at least partly due to the fact that $P_{U\mid Y,X}$ is an infinite dimensional object, even in the case when $X$ has finite support. This infinite dimensional existence problem is exacerbated in practice by the fact that $P_{U\mid Y,X}$ must satisfy a number of constraints to ensure it is consistent with the binary response model through (ref), and to ensure it is a proper conditional probability measure. We consider these practical difficulties in detail in the next section.

General Framework: Practical Considerations

To make progress, define the following vector-valued function:

and for a fixed binary vector $s \in \{0,1\}^{m}$ define the set:

The sets from (ref) partition the space $\mathcal{U}$ into at most $2^{m}$ sets, with each set uniquely associated with a binary vector $s \in \{0,1\}^{m}$. The binary vectors $r(u,\theta)$ represent response types (c.f. balke1994counterfactual and heckman2018unordered).\footnote{The collection of sets defining response types appears to be similar to the “minimal relevant partition” in tebaldi2019nonparametric, as well as the partition described in chesher2014instrumental Appendix B. } In the discrete choice setting, these response types tell us the choices that an individual with type indexed by $(u,\theta)$ would have made had they been assigned an alternate value of $x$. Any two individuals characterized by values of $u$ from the same set $\mathcal{U}(s,\theta)$ make identical choices under every possible counterfactual assignment $\gamma: \mathcal{X} \to \mathcal{X}$, making this a natural grouping of latent types.

After partitioning the space of latent variables using response types, various counterfactual objects of interest can be written as a disjoint union of the sets $\mathcal{U}(s,\theta)$ from (ref) that comprise the partition. For the sake of illustration, consider the binary vectors:

for $j=1,\ldots,m$. Note that each set $S_{j}$ is comprised of all binary vectors that have a $j^{th}$ entry equal to $1$, and thus contain exactly $2^{m-1}$ elements.\footnote{It is useful to note that the sets $\{S_{j}\}_{j=1}^{m}$ are not disjoint; indeed, it is easy to show that $S_{j} \cap S_{k} \neq \emptyset$ and $S_{j} \cap S_{k}^{c} \neq \emptyset$ for every $j\neq k$.} By definition of the sets $\mathcal{U}(s,\theta)$ and $S_{j}$ we have:

Furthermore, for $s' \neq s$ we have $\mathcal{U}(s',\theta) \cap \mathcal{U}(s,\theta) = \emptyset$, so that this union is a disjoint union. Thus, we have the following decomposition:

Such a decomposition holds for any $j=1,\ldots,m$. When the conditioning value $x$ differs from the value $x_{j}$ in the structural function, an application of Theorem (ref) shows that the left hand side of this display represents a counterfactual conditional distribution, illustrating the connection between response types and counterfactual choices.

The next result shows that, in order to rationalize a given collection of counterfactual conditional distribution, for each $\theta$ it is both necessary and sufficient to construct a probability measure on sets of the form $\mathcal{U}(s,\theta)$ from (ref) that satisfies the constraints of Theorem (ref). In the statement of the result we redefine $\gamma:\mathbb{N}\to \mathbb{N}$ to denote the index of the point in $\{x_{1},,\ldots,x_{m}\}$ assigned under counterfactual $\gamma$, and we set $S_{\gamma(j)} := \{ s \in \{0,1\}^{m} : s_{\gamma(j)} =1 \}$, the analog of $S_{j}$ from (ref).

Theorem (ref) reduces our infinite dimensional existence problem to a finite dimensional existence problem. For a fixed $\theta \in \Theta$, instead of verifying there exists a Borel probability measure $P_{U \mid Y, X}$ satisfying (ref) and (ref) (an infinite-dimensional object), Theorem (ref) shows it suffices to verify there exists a finite dimensional probability vector with typical element $P_{U\mid Y,X}\left( \text{int}(\mathcal{U}(s,\theta)) \mid Y=y, X=x\right)$ satisfying (ref) - (ref). Note that this result relies crucially on the finiteness of $\mathcal{X}$. There is a close connection between Theorem (ref) and the bounding approach based on Artstein's inequalities (e.g. chesher2017generalized) and optimal transportation (e.g. galichon2011set).\footnote{In Appendix (ref) we show that Theorem (ref) is equivalent to a characterization based on Artstein's inequalities after conditioning on the value of the endogenous variables. This conditioning allows us to obtain a much smaller number of equality constraints when compared to the full set of unconditional constraints arising from Artstein's inequalities. However, Artstein's inequalities can handle continuous instruments. It is also well known that Artstein's inequalities are equivalent to the existence of a certain zero-cost optimal transport problem (see galichon2009test, galichon2011set, and galichon2016optimal). } Related results have also appeared in laffers2019identification and torgovitsky2019partial. The finite number of linear constraints from Theorem (ref) leads naturally to a linear programming approach to bounds on counterfactual quantities.

Optimization Formulation

For simplicity, we suppose throughout this subsection that our objective is to bound the counterfactual probability:

for some $j \in \{1,\ldots,m\}$. All the results in this section are immediately applicable to the case when we wish to bound some linear function of these counterfactual probabilities. From Theorem (ref) our counterfactual probability can be written as:

where $\gamma(j)$ is the index in $\{1,\ldots,m\}$ assigned to $j$ under counterfactual $\gamma$. Define the parameter:

Furthermore, define:

The vector of parameters $\pi(\theta)$ represents the variable over which we optimize in our result ahead, and has dimension $d_{\pi} = 2m2^m$. Without loss of generality, we suppose that each $(y,x)$ is assigned positive probability by the observed distribution. From conditions (ref) and (ref) in Theorem (ref), we have the constraints:

for $j=1,\ldots,m$. Finally, we require the constraints:

for all $y\in\{0,1\}$ and $j=1,\ldots,m$ and $s \in \{0,1\}^{m}$, and:

for all $y \in \{0,1\}$ and $j=1,\ldots,m$. Note that to impose the constraints from (ref), the researcher must first determine which sets $\mathcal{U}(s,\theta)$ have nonempty interior. We return to this point in the next subsection. We are now ready to state one of the main results.

In one direction, Theorem (ref) implies that any counterfactual conditional probability of the form (ref) belonging to the identified set can be written as:

for some $\theta$ and some vector $\pi(\theta)$ satisfying the constraints (ref), (ref), and (ref). In the opposite direction, the Theorem implies that if for some $\theta$ the vector $\pi(\theta)$ satisfies the constraints (ref), (ref), and (ref) then the conditional probability measure on $\mathcal{U}$ represented by $\pi(\theta)$ can be extended to a (not necessarily unique) Borel probability measure on all of $\mathfrak{B}(\mathcal{U})$ that satisfies the conditions of Theorem (ref). This result can be easily modified to bound any linear function of counterfactual conditional probabilities by simply modifying the objective function in Theorem (ref). We make use of this fact in the application section.

After determining which of the sets $\mathcal{U}(s,\theta)$ have nonempty interior, all the constraints in (ref) and (ref) can be written as linear equality/inequality constraints, so that the optimization problems in (ref) and (ref) are linear programming problems. This is very beneficial, since linear programs can be efficiently solved even in cases with thousands of parameters and constraints.

Interestingly, the proofs for Theorems (ref) and (ref) do not require linearity of the index function in $U$ from Assumption (ref), although this assumption is used starting in the next subsection. This implies that Theorem (ref) can be used to bound counterfactual parameters for models of the form:

without imposing any restrictions on the index function. Without any restrictions, it is always possible to construct a function $\varphi$ such that all regions $\mathcal{U}(s,\theta)$ have nonempty interior, implying that there are $2^m$ response types. In this case, constraint (ref) only imposes that all probabilities are bounded between zero and one, and the rest of Theorem (ref) remains unchanged. We illustrate our procedure using the model (ref) in the application section.

In the general case, elements of $\pi(\theta)$ corresponding to sets $\mathcal{U}(s,\theta)$ with empty interior can be removed from the parameter vector $\pi(\theta)$ without altering the optimal solutions to the linear programs in (ref) and (ref). This allows for further reduction of the dimension of these linear programs. Although the number of sets $\mathcal{U}(s,\theta)$ appear to grow exponentially in $m$, in the subsections ahead we show that under linearity of the index function in $U$ the number of sets $\mathcal{U}(s,\theta)$ that have nonempty interior grows at a rate that is polynomial in $m$, substantially reducing the computational burden.

An Impossibility Result

To this point, the binary response model has been left almost entirely unrestricted. The following corollary to Theorem (ref) states formally the unsurprising result that, when $\varphi$ is completely unrestricted, counterfactual choice probabilities are not identified under our assumptions.

This result is stated as a corollary since it can be proven using Theorem (ref). Intuitively, Theorem (ref) shows that under Assumptions (ref) and (ref) the observed conditional choice probabilities effectively impose no constraints on counterfactual conditional choice probabilities; indeed, this follows by simple inspection of the constraints (ref) and (ref) from the statement of Theorem (ref), as well as the discussion following the statement of Theorem (ref). Without any constraints imposed by the observed conditional choice probabilities, by choosing a sufficiently flexible function $\varphi:\mathcal{X} \times \mathcal{Z} \times \Theta \times \mathcal{U} \to \mathbb{R}$ it is possible to rationalize all counterfactual conditional choice probabilities in the interval $[0,1]$. This result relies on the form of the counterfactual choice probabilities in (ref) and (ref); for other counterfactual parameters---e.g. the average treatment effect (ref)---it is possible to obtain nontrivial bounds under only Assumptions (ref) and (ref).\footnote{It is important to note that this does not occur because Assumptions (ref) or (ref) have imposed any structure, but instead is because the average treatment effect is a weighted average of unidentified counterfactual conditional choice probabilities and identified conditional choice probabilities.}

This impossibility result shows that it is necessary to entertain additional assumptions in order to obtain informative bounds on counterfactual choice probabilities. In the next section, we will explore additional assumptions that can be used to tighten bounds on counterfactual choice probabilities.\footnote{For additional impossibility results in partially identified discrete response models, see the discussion on pages 1396 and 1402 of manski2007partial as well as Corollary 3 in chesher2013instrumental.} \end{comment}

Hyperplane Arrangements and Cell Enumeration

To use Theorem (ref), we must first determine which set $\mathcal{U}(s,\theta)$ have nonempty interior in order to impose the constraints from (ref). This section explains how this can be done in our setting by using linearity of the index function in the latent variables in combination with the cell enumeration algorithm of gu2020nonparametric.

First we demonstrate that linearity of the index function in the latent variables imposes restrictions on the model by limiting the number of sets $\mathcal{U}(s,\theta)$ that can be assigned positive probability. Constraining sets of the form $\mathcal{U}(s,\theta)$ to be assigned zero probability is called eliminating response types. Response types corresponding to sets $\mathcal{U}(s,\theta)$ that survive elimination are called admissible, and response types corresponding to sets $\mathcal{U}(s,\theta)$ that are eliminated are called inadmissible. The following simple example shows how we can eliminate response types under Assumption (ref).

In the general case, it can be shown that when $\varphi$ is restricted to be linear in $U$, there is an upper bound on the number of nonempty sets $\mathcal{U}(s,\theta)$ that grows at a rate that is polynomial in $m$ rather than exponential, which is the case when $\varphi$ is unrestricted.

This result is implied by results in the literature on combinatorial geometry. In particular, linearity of the function $\varphi(\,\cdot\,,u)$ means that for each instance of $(x,\theta)$ the function $\varphi(x,u,\theta)$ defines a hyperplane in $d_{u}-$dimensional space. In the case when the vectors defining these hyperplanes are in general position the upper bound in Proposition (ref) is obtained.\footnote{ A collection of $m$ hyperplanes in $d-$dimensional space are considered to be in general position when any collection of $k$ out of the $m$ hyperplanes intersect in a $d-k$ dimensional space for $1 < k \leq d$, and any collection of $k$ out of $m$ hyperplanes has an empty intersection for $k > d$.} This latter result was initially proven by Buck. Straightforward calculation shows that, for $m>d_{u}+1$, the function $\sum_{j=0}^{d_{u}} \binom{m}{j}$ is bounded above by $(e \cdot m/d_{u})^{d_{u}} = O(m^{d_{u}})$, a polynomial in the number of hyperplanes.

Define the collection:

To impose linearity in the latent variables we must determine which sets $\mathcal{U}(s,\theta)$ have nonempty interior, and then ensure that any distribution of the latent variables assigns zero probability to these sets. To compute the collection $S_{\varphi}(\theta)$ we propose to use the enumeration algorithm of gu2020nonparametric.

When the index function $\varphi$ is linear in $U$, for each fixed $\theta$ and $s \in \{0,1\}^{m}$ the set $\mathcal{U}(s,\theta)$ is a convex polyhedron formed by the intersection of at most $m$ halfspaces whose boundaries are hyperplanes of the form $\{u \in \mathcal{U} : \varphi(x,u,\theta)=0\}$. The hyperplane arrangement algorithm of gu2020nonparametric accepts $m$ hyperplanes as an input, and outputs the binary vectors $s$ corresponding to the sets $\mathcal{U}(s,\theta)$ that have nonempty interior, as well as a point from each of these sets. avis1996reverse were the first to provide an enumeration algorithm that runs in a time proportional to the maximum number of sets with nonempty interior. Improvements to the enumeration algorithm were made by sleumer1999output and rada2018new. The algorithm of gu2020nonparametric was developed for the problem of nonparametric maximum likelihood in a linear index random coefficient model. It runs in a time proportional to $O(m^{d_{u}})$, which is near-optimal when the hyperplanes are in general position, according to Proposition (ref).

To understand the algorithm, note that for each $s \in \{0, 1\}^m$ and fixed $\theta$, we can verify using a linear program whether there exists a point in the interior of $\mathcal{U}(s,\theta)$. Consider the following problem:

where $s_{j}$ is the $j^{th}$ element of our fixed binary vector $s$. If $\varepsilon^{*}$ and $u^{*}$ are the optimal values of the program (ref), then a value $\varepsilon^{*}>0$ indicates that $u^{*}$ is an interior point to the polyhedron $\mathcal{U}(s,\theta)$.

Since the linear program (ref) must be solved for each $s \in \{0, 1\}^m$, checking whether each $\mathcal{U}(s,\theta)$ has nonempty interior requires solving $2^{m}$ linear programs, despite the fact that we know the number of nonempty subsets $\mathcal{U}(s,\theta)$ is polynomial in $m$. To address this issue, the algorithm proposed in gu2020nonparametric builds upon the algorithm in rada2018new. The idea is to add one hyperplane at a time. At step $k$ we start with a collection of $k-1$ hyperplanes from the previous steps, as well as all existing response types found up to step $k-1$. We then introduce a new hyperplane into the arrangement, and determine all newly created response types by solving a collection of linear programs. When a new hyperplane is added the only new cells are those that are created when the existing cells are crossed by the last hyperplane. By efficiently locating those crossed cells at each step, all cells in the full arrangement can be enumerated in polynomial time.

In summary, the hyperplane arrangement algorithm is used as a pre-processing step under Assumption (ref) to determine which sets $\mathcal{U}(s,\theta)$ have nonempty interior in a given application. Eliminating the inadmissible sets $\mathcal{U}(s,\theta)$ can also dramatically reduce the dimension of the parameter vector $\pi(\theta)$ in the bounding optimization problems. In particular, under Assumption (ref) we need only consider a parameter vector $\pi(\theta)$ with typical element $\pi(y,x,s,\theta)$ defined for $s$ corresponding to subsets $\mathcal{U}(s,\theta)$ with nonempty interior. The dimension of the revised parameter vector $\pi(\theta)$ constructed in this way is always upper-bounded by a polynomial in $m$ under Assumption (ref). In the next subsection we show how the assumption of linearity in parameters $\theta \in \Theta$ can be combined with the hyperplane arrangement algorithm to dramatically simplify the bounding procedure suggested by Theorem (ref).

Profiling Under Linearity in the Fixed Coefficients

Constructing bounds on counterfactual probabilities using Theorem (ref) requires evaluating the linear programs (ref) and (ref) at all values of $\theta \in \Theta$ in the parameter space. In practice this procedure is infeasible, and instead the identified set must be constructed using Theorem (ref) by establishing a grid over $\Theta$. The following proposition demonstrates that, theoretically speaking, the researcher need only consider a finite grid.

We call the points in $\Theta'$ the representative points, although these points are generally not unique. If the representative points can be determined by the researcher, Proposition (ref) implies that the union over $\theta \in \Theta$ in (ref) can be replaced with a union over $\theta \in \Theta'$. That is, the linear programs in (ref) and (ref) need only be solved at the representative points.\footnote{See the discussion at the end of Section 3 in torgovitsky2019partial for the same idea. Proposition (ref) also implies that the identified set for counterfactual conditional distributions in Theorem (ref) is always a closed (but possibly disconnected) set.}

In the case when $\varphi$ is also linear in $\theta$, we provide a polynomial-time algorithm for finding a collection of representative points.\footnote{Linearity in $\theta$ is not necessarily restrictive. Since $X$ is discrete, the researcher can construct a saturated specification with an index function $\sum_{x \in \mathcal{X}} \mathbbm{1}\{X=x\}\theta_{x} + U$. Defining $\theta :=(\theta_{x})_{x \in \mathcal{X}}$, this index function is linear in $(\theta,U)$, and so all of our computational results are applicable.} To introduce the approach, recall the set $S_{\varphi}(\theta)$ from (ref). Note that for any two values of $\theta, \theta' \in \Theta$ with $\theta \neq \theta'$, if $S_{\varphi}(\theta) = S_{\varphi}(\theta')$ then the linear programming problems in Theorem (ref) at $\theta$ and $\theta'$ are identical, since they have an identical set of constraints. The points $\theta$ and $\theta'$ are thus equivalent in the sense that we only need to solve the linear programming problems for one of them. Extending this idea, we can define an equivalence class by the set of all $\theta \in \Theta$ delivering the same collection $S_{\varphi}(\theta)$. We then only need to solve the linear programming problems at one value of $\theta$ belonging to each equivalence class. These values of $\theta$ selected from each equivalence class are exactly what we call representative points.

To see how to find the representative points, partition $x = (x_{r},x_{f})$ and suppose the index function is of the form:

where $x_{r}$ has dimension $d_{u}$ and $x_{f}$ has dimension $d_{f}=d_{x} - d_{u}$. Here $x_{r}$ is a subvector of $x$ associated with a random coefficient $u \in \mathcal{U}$, and $x_{f}$ is a subvector of $x$ associated with a fixed coefficient $\theta \in \Theta$. For any binary vector $s \in \{ 0,1\}^{m}$, define:

These sets form a unique partition of the space $\mathcal{U}\times \Theta$ into cells defined by $m$ hyperplanes of the form:

To find the representative points, we project the sets $\mathcal{R}(s)$ onto the parameter space $\Theta$, and intersect the projections of each set $\mathcal{R}(s)$ across $s$. This produces a collection of sets on $\Theta$ corresponding exactly to the equivalence classes discussed above. The main challenge of the procedure is to obtain a tractable characterization of the projected sets.

Define the set:

That is, $S_{\varphi}$ collects all vectors $s \in \{0,1\}^{m}$ corresponding to sets $\mathcal{R}(s)$ with nonempty interior. We begin by determining all vectors $s \in S_{\varphi}$ by running the hyperplane arrangement algorithm of gu2020nonparametric on the $m$ hyperplanes of the form (ref). These $m$ hyperplanes can be combined as $X_{r} u + X_{f} \theta =0$, where $X_{r}$ is an $m \times d_{u}$ matrix and $X_{f}$ is an $m\times d_{f}$ matrix. Now fix any $s \in S_{\varphi}$, and let $D(s) = \text{diag}(2s-1)$ denote the $m\times m$ diagonal matrix with the sign vector $2s-1$ along its main diagonal. Define $X_{r}(s) := D(s) X_{r}$ and $X_{f}(s):= D(s) X_{f}$. Then $\mathcal{R}(s)$ can be rewritten as:

The projection of $\mathcal{R}(s)$ on $\Theta$ is:

The objective is now to write $\Theta(s)$ only in terms of linear inequality constraints in $\theta$; in other words, to “eliminate” the latent variables $u$ from the system of inequalities in (ref).\footnote{It is possible to use Fourier-Motzkin elimination for this purpose, which was explored in a similar context in Section 8.2 of chesher2019generalized. In practice we find that Fourier-Motzkin elimination leads to a prohibitively large number of constraints defining (ref), most of which are redundant. This remains true even when the number of inequalities defining the set (ref) is small.} To this end, consider the set:

kohler1967projections showed that the projected set (ref) can be rewritten as:

Furthermore, the Minkowski-Weyl Theorem also allows us to re-write the set $\mathcal{C}(s)$ as:

where $R(s)$ is some matrix.\footnote{For a general convex set defined by $\Lambda = \{\lambda \in \mathbb{R}^d: A\lambda \leq b\}$, the Minkowski-Weyl Theorem states that every vector $\lambda \in \Lambda$ can be written as $\lambda = \lambda_{1}+\lambda_{2}$, where $\lambda_{1} \in \text{conv} \{ v_1, \dots, v_k\}$ and $\lambda_{2} \in \text{cone}\{v_{k+1}, \dots, v_n\}$. Here $v_1, \dots, v_k$ are called vertices of $\Lambda$ and $v_{k+1}, \dots, v_n$ are the extreme rays of $\Lambda$. In the special case of $b = 0$, where all hyperplanes are through the origin, then $\Lambda$ becomes a polyhedral cone and $k=0$, so that $\Lambda = \text{cone}\{v_1, \dots, v_n\}$. This latter case is what is relevant for us, and the columns of the matrix $R(s)$ are the collections of these extreme rays. } That is, every element belonging to the polyhedral cone $\mathcal{C}(s)$ can be written as a nonnegative linear combination of the columns of $R(s)$. The matrix $R(s)$ is called the generating matrix of the polyhedral cone $\mathcal{C}(s)$, and the problem of finding the minimal generating matrix $R(s)$ is called the extreme ray enumeration problem.\footnote{A minimal generating matrix for $\mathcal{C}(s)$ is a generating matrix with the property that no proper submatrix also generates $\mathcal{C}(s)$ (c.f. Fukuda). Note the minimal generating matrix is unique only up to multiplication by a positive scalar.} After obtaining the matrix $R(s)$, we have the following representation of the projection of $\mathcal{R}(s)$:

The two representations of the cone $\mathcal{C}(s)$ from (ref) and (ref) are called its H-representation and its V-representation, respectively. Converting from one representation of a convex polyhedron to another is called the double description problem in computational geometry. An efficient double description algorithm was proposed by Fukuda, and we use the R implementation in the package Rcdd by Geyer. avis1997good provides a comparison of different algorithms. The projection procedure is then repeated for all $s \in S_{\varphi}$, resulting in the collection $\{ \Theta(s) : s \in S_{\varphi}\}$. To obtain the final representative points, we stack $R(s)^\top X_f(s) \theta = 0$ for all $s \in S_{\varphi}$, and then run the enumeration algorithm of gu2020nonparametric on this final collection of hyperplanes.

This procedure also sheds light on how to construct the identified set for $\theta$. Since our final arrangement involves only hyperplanes through the origin, each representative point is selected from a polyhedral cone. For some of these representative points, the linear programming problems in our bounding procedure from Theorem (ref) may have an empty feasible region. In these case, the representative points---and all points belonging the corresponding cone---lie outside of the identified set $\Theta^{*}$. Thus, the identified set $\Theta^*$ is the union of the polyhedral cones with representative points associated with nonempty feasible regions in our linear programming problems from Theorem (ref).\footnote{This implies that the identified set $\Theta^*$ may not be connected, and for any $\theta \in \Theta^*$, we also have $\lambda \theta \in \Theta^*$ for all $\lambda \geq 0$. An appropriate normalization---for example, fixing $||\theta||=1$---leads to a bounded identified set $\Theta^*$.}

Additional Assumptions

In this section we describe how to impose additional independence and monotonicity assumptions in our framework.

Independence Assumptions

In some cases the researcher may have access to an observed variable that is independent of the latent variables. If such a variable enters as an argument in the index function, it induces variation in the observed conditional probabilities without affecting the distribution of the latent variables. We refer to such variables as exogenous covariates. A similar intuition applies if the variable is independent of the distribution of latent variables, does not enter as an argument in the structural function, but has nontrivial dependence with the variables that do enter the index function. We refer to such variables as instruments. Both exogenous covariates and instruments can be used to produce a smaller identified set for counterfactual parameters. To introduce our independence assumption, we partition $X=(W,Z)$, where $W \in \mathcal{W} \subset \mathbb{R}^{d_{w}}$ and $Z \in \mathcal{Z} \subset \mathbb{R}^{d_{z}}$ and $d_{w}+d_{z}=d_{x}$.

The independence assumption constrains the set of admissible latent variable distributions, and links the conditional distributions of $U \mid Z=z$ across values of $z \in \mathcal{Z}$. Although Assumption (ref) posits full independence between $Z$ and the vector of latent variables $U$, the assumption can be easily modified for the case when a subvector of $Z$, say $Z_{1}$, is conditionally independent of $U$ given some other subvector of $Z$, say $Z_{2}$.

Definition (ref) in Appendix (ref) provides the extension of Definition (ref) to the case when Assumption (ref) also holds, and Corollary (ref) in Appendix (ref) provides the analogous extension of Theorem (ref). To extend the linear programming result of Theorem (ref) we must include additional constraints. Without loss of generality we assume all values of $(y,w,z)$ are assigned positive probability. Then the appropriate constraints on the vector $\pi(\theta)$ are given by:

for $k=1,\ldots,m_{z}-1$, where $m_{z} = |\mathcal{Z}|$. A formal statement of the extension of Theorem (ref) to the case when the constraints (ref) are also imposed is provided by Corollary (ref) in Appendix (ref).

Monotonicity Assumptions

Let $\mathcal{M} \subset \{1,\ldots,m\}\times \{1,\ldots,m\}$ denote any collection of integer tuples $(j,k)$, where $1 \leq j,k \leq m$.

This monotonicity assumption states that, when comparing two points $x_{j}$ and $x_{k}$, the value of the structural function can be ordered by the researcher. For instance, if the order determined by the researcher's monotonicity assumption for the points $x_{j}$ and $x_{k}$ is $\varphi(x_{j},u,\theta) \leq \varphi(x_{k},u,\theta)$, then the researcher automatically rules out response types with $\mathbbm{1}\{\varphi(x_{j},u,\theta) \geq 0\} > \mathbbm{1}\{\varphi(x_{k},u,\theta)\geq 0\}$.

Similar monotonicity assumptions in triangular systems have been extensively explored by heckman2018unordered. In particular, heckman2018unordered explore how choice theory can be used to impose monotonicity assumptions and to eliminate response types, and many of their insights are applicable here.

Let $S_{M}$ be the collection of vectors $s \in \{0,1\}^{m}$ that respect the monotonicity relations from Assumption (ref). Definition (ref) in Appendix (ref) provides the extension of Definition (ref) to the case when Assumption (ref) is also imposed. The extension of Theorem (ref) to the case when Assumption (ref) is imposed is provided by Corollary (ref) in Appendix (ref). To extend the results of Theorem (ref) we must simply include the set of constraints imposed by Assumption (ref) in our optimization problems. These constraints are provided in Corollary (ref), and can be written in terms of the parameter vector $\pi(\theta)$ as:

for all $y\in \{0,1\}$ and $j=1,\ldots,m$ occurring with positive probability. Corollary (ref) in Appendix (ref) then shows the extension of Theorem (ref) to the case when Assumption (ref) is imposed using the constraints (ref).

Consistency, Inference, and Bias Correction

The previous section focused on identification and computation issues, and set up the main optimization procedure for bounding counterfactual quantities. In this section we briefly describe estimation and inference in our setting. Our inference procedure is adapted from the procedure of cho2021simple.

In most applications of our procedure, we expected both the constraints and objective function in the linear programming problems (ref) and (ref) to be data dependent, and our results in this section are designed with this case in mind. We make use of this in the application section. Although motivated by the previous sections, the results in this section do not require any of our previous assumptions to hold. In this sense this section is independent of the previous sections, and may be of interest to researchers working with a similar class of problems.

Let $\mathscr{P}$ denote the set of all probability measures on $\mathcal{Y}\times \mathcal{X}$, and let $\Pi:=\{\pi \in \mathbb{R}^{d_{\pi}} : A \pi \leq b\}$ denote a compact and convex polytope defined by some matrix $A$ and some vector $b$. In our context, $\Pi$ represents the parameter space for the conditional probability vector $\pi$ from Theorem (ref). To introduce our results, we convert each of the equality constraints in the linear programs from Theorem (ref) to two equivalent inequality constraints. Then for a fixed $\theta \in \Theta$, the constraints in the programs (ref) and (ref) can be written as a finite number of moment inequality constraints:

where $\mathcal{J}(\theta)$ is an index set that may depend on $\theta$.\footnote{These moment functions should also include parameter space constraints.} Note this formulation also allows the use of data-dependent constraints in the programs (ref) and (ref), so long as these constraints can be written as a collection of moment inequalities. Furthermore, the objective function in the programs (ref) and (ref) can be written as $\E_{P}[\psi(Y_{i},X_{i},\pi,\theta)]$, which also permits the use of a data-dependent objective function. We require the following assumption.

Part (i) restricts the objective function to be Lipschitz continuous in $\pi$ for each $\theta$. This assumption is satisfied, for instance, when bounding counterfactual probabilities or average treatment effects, as in the application section. Part (ii) requires each moment function to be linear in $\pi$, an assumption which is also satisfied for the model and assumptions in this paper. Assumption (iii) is standard in the moment inequalities literature, and constrains $\mathcal{P}$ to be the collection of data generating processes that satisfy the moment conditions. Part (iv) imposes a uniform integrability requirement on the vector $(Y,X)$ needed for the procedure of cho2021simple, and part (v) imposes a uniform integrability condition on the gradients. These are both easily satisfied with discrete random variables $Y$ and $X$. Part (vi) allows us to replace $\Theta$ with a finite subset $\Theta'$ without impacting the bounding problem. Proposition (ref) shows this is the case in our setting. Part (vii) assumes the existence of a bounded envelope function on a slight expansion of $\Pi$. This is also easily verified, for instance, when bounding counterfactual probabilities or average treatment effects. Finally part (viii) assumes that the sample under consideration is i.i.d. from some $P \in \mathcal{P}$ satisfying the other conditions.

Using these assumptions we will prove a consistency result, and demonstrate how to do inference in our class of problems. For any $c \in \mathbb{R}_{+}$ and any $P \in \mathscr{P}$ let us define:

where:

The value functions $\Psi_{\ell b}(\theta,P,0)$ and $\Psi_{u b}(\theta,P,0)$ are analogous to the value functions $\pi_{\ell b}(y,x,\theta)$ and $\pi_{u b}(y,x,\theta)$ from Theorem (ref). Furthermore, let $\Psi^{*}(P):=\Psi^{*}(P,0)$, and denote the empirical measure as $\mathbb{P}_{n} \in \mathscr{P}$. The following theorem shows that a slight (shrinking) enlargement of the set $\Psi^{*}(\P_{n})$ is a consistent estimator for the set $\Psi^{*}(P)$, where consistency is defined using the Hausdorff metric.\footnote{For two sets $A, B \subset \mathbb{R}^{d}$, the Hausdorff metric is defined as:

}

Proposition (ref) is related to results and discussions found in molchanov1998limit, manski2002inference, and chernozhukov2007estimation. It is also a special case of a more general consistency result presented in Appendix (ref), and suggests a consistent estimator for the identified set from Theorem (ref). In the application section we take $b_{n} = b/\sqrt{\log(n)}$ for some $b>0$, and call $\Psi^{*}(\P_{n},b_{n})$ the “plug-in” estimate of the identified set.

In Section (ref) we also use the inference method of cho2021simple, designed for uniform inference on value functions in stochastic linear programming problems. In general, this inference problem is highly irregular, but cho2021simple show that regularity can be restored by introducing infinitesimal random perturbations to the constraints and objective function. After perturbing the problem, they prove consistency of a simple and fast nonparametric bootstrap procedure to construct a confidence set. However, a modification of the procedure of cho2021simple is needed to fit our setting to allow for the profiling points $\theta \in \Theta'$.

Let $\xi\sim P_{\xi}$ denote the random perturbation vector from the procedure of cho2021simple. By Proposition (ref), there exists a finite set $\Theta' \subset \Theta$ of representative points satisfying:

To apply the procedure of cho2021simple, we take each $\theta \in \Theta'$, and use the procedure of cho2021simple to construct a set $C_{n}(1-\alpha,\theta)$ satisfying:

Here, probability is taken with respect to the product measure $\text{Pr}_{P}\times P_{\xi}$ to account for the random perturbations $\xi \sim P_{\xi}$ introduced to restore regularity. We refer to cho2021simple for additional discussion. We then set:

The following result shows that this confidence set is uniformly valid over $\mathcal{P}$. The proof of the result proceeds by verifying the assumptions of cho2021simple for (ref), and then showing that the modified confidence set $CS_{n}(1-\alpha)$ has the correct coverage.

In our application ahead, we also report bias-corrected estimates of the lower and upper endpoints of the identified set. Convexity of the minimum in combination with Jensen's inequality shows that the sample analog lower bound is biased upward. Similarly, the sample analog upper bound is biased downward. This leads to an identified set that is on average too narrow. In response, chernozhukov2013intersection proposed the use of half-median unbiased estimators. Half-median unbiased estimates $\hat{\Psi}_{\ell b}$ and $\hat{\Psi}_{ub}$ satisfy $\hat{\Psi}_{\ell b} \leq \Psi_{\ell b}(P)$ and $\Psi_{ub}(P)\leq \hat{\Psi}_{u b}$, each holding with probability at least $1/2$. We construct a bias-corrected estimate of $\Psi_{\ell b}(P)$ ($\Psi_{u b}(P)$) using the inference procedure of cho2021simple by setting $\alpha=1/2$ and by taking our estimate to be the endpoint of the $1-\alpha$ lower (upper) confidence set for $\Psi_{\ell b}(P)$ ($\Psi_{u b}(P)$). In our application we report both the plug-in estimates of our bounds based on Theorem (ref) and the half-median unbiased estimates.

Application

In this section we apply our method to study the impact of private health insurance on an individual's decision to visit a doctor. In general, insurance markets are plagued by problems arising from asymmetric information between consumers and insurance providers (c.f. rothschild1978equilibrium). For example, adverse selection occurs in the health insurance market when individuals have more information about their latent health determinants than the providers of health insurance. A robust prediction of the classical theory of asymmetric information is that those who are more likely to purchase insurance are also those who are more likely to experience the insured risk.\footnote{The “insured risk” refers to the event for which insurance was purchased. In our context, it is any event that would typically require a visit to the doctor.} On the other hand, there has been little and mixed empirical evidence of adverse selection in health insurance markets (see cardon2001asymmetric for a discussion).

In this section we compute various counterfactual parameters while remaining agnostic on the exact nature of the latent variables linking health insurance and health care utilization decisions. We take the decision to visit a doctor as our binary outcome variable of interest, and consider the individuals' private health insurance status as an endogenous explanatory variable. This is consistent with the idea that private insurance status may be dependent with individual-specific latent factors---most importantly, unobserved health determinants and attitudes towards risk---that influence an individual's propensity to visit a doctor. We use data from the 2010 wave of the Medical Expenditure Panel Survey (MEPS), which has also been recently analyzed by han2019estimation and acerenza2021testing. We focus on the same sub-sample considered in these papers. In particular, we focus on the month of January 2010, consider only individuals between ages $25$ and $64$, and drop individuals who obtain either federal or state insurance in 2010 and individuals who are self-employed or unemployed. These restrictions leave us with a sample of $7555$ individuals.

In all specifications $W$ is a binary endogenous variable representing an individual's private insurance status, and we consider a binary health status variable ($Z_{1}$) and a binary marital status variable ($Z_{2}$) as regressors.\footnote{The MEPS data includes information on self-reported health status on a scale from $1- 5$, and we consider values less than or equal to $2$ as being “unhealthy.”} Finally, we use the number of employees working for the individual's firm ($Z_{3}$) as an instrument. This variable provides a measure of the size of a firm and has discrete support in the range $[1,500]$, which we further discretize into 11 bins.\footnote{Variable $Z_3$ is supported on the range $[1, 500]$ and is clearly top-coded. We notice that there is bunching of observations at firm sizes in multiples of five, and some regions of the support of $Z_3$ contain very few observations. In order to get reliable estimates of the conditional choice probabilities, we further discretize the firm size into 11 bins. The bins are respectively $[1, 5]$, $(5, 10]$, $(10,20]$, $(20,30]$, $(30,40]$, $(40,50]$, $(50,60]$, $(60,70]$, $(70,100]$, $(100,200]$ and $(200,500]$.} Using firm size as an instrument is consistent with the evidence that larger firms are more likely to provide health insurance benefits, but do not directly influence an individual's decision to visit a doctor.\footnote{From cardon2001asymmetric p.408: “Another observed symptom, consistent with the theoretical predictions, is that the uninsured tend to work for small employers. Large employers can overcome adverse selection by risk pooling.” }

A possible concern with using firm size as an instrument is that risk averse individuals may be more likely to select into a job with a larger firm size. In an attempt to address this issue, we include an alternate independence assumption that assumes the firm size $Z_3$ is conditionally independent of $U$ given $(Z_1, Z_2)$ only when $Z_{3}$ lies within a certain range. The idea is that once we condition on a particular range of firm size, the remaining variation in firm size is independent of $U$ conditional on $(Z_1, Z_2)$. We consider four ranges, given by $(1, 10]$, $(10, 50]$, $(50, 100]$ and $(100,500]$, and impose our conditional independence assumption for each range separately.

The first parameter we consider is the average treatment effect, defined as:

This parameter provides the average causal effect of obtaining health insurance on the decision to visit a doctor. Second, we consider the counterfactual choice probability:

for $y \in \{0,1\}$. We focus on the parameter $\mu_{ccp}(0)$ for simplicity, which represents the counterfactual choice probability of visiting a doctor when given private health insurance for the set of individuals who have no insurance and who have chosen not to visit a doctor, averaged across health and marital status. We construct our bounds under the following set of assumptions:

Note that the general index function takes the form $\varphi(w,z_{1},z_{2},u,\theta)$. When monotonicity is imposed in (A2), we impose:

for each $z_{2} \in \{0,1\}$. This implies that for an unhealthy individual, the propensity to visit a doctor when the person has private insurance is always weakly greater than without insurance, regardless of marital status. Finally we consider three different models for the binary outcome variable $Y$:

Recall that the extension of our procedure to cover model (ref) was discussed briefly at the end of Section (ref). Indeed, under model (ref) the index function $\varphi$ need not be explicitly specified and it may not satisfy the linearity assumption made under Assumption (ref). This makes model (ref) the most flexible. Models (ref) and (ref) impose linearity of $\varphi$ in the latent variables and in the parameters. Here we distinguish two cases. In the first case, (ref) regards $(U_{1},U_{2})$ as the latent variables in the model. Model (ref) is the same as (ref) except that we have replaced the random slope coefficient $U_{1}$ from (ref) with a fixed coefficient. Model (ref) represents the additively separable linear index model that is commonly used in the empirical literature, except for the fact that we do not assume a parametric distribution for $U$ and do not have a model for the endogenous variable $W$.

The identified sets for $\mu_{ate}$ under assumptions (A1) - (A8) and models (ref) - (ref) are reported in Table (ref). For simplicity, we report the convex hull of the estimated identified set for each specification. Table (ref) also reports our modified plug-in estimator (see Section (ref)) as well as half-median unbiased estimators and $90\%$ confidence sets constructed using the modified procedure of cho2021simple. Due to a confluence of factors---including the dimension of the empirical choice probability vector, the large number of constraints, and the sample size---we find that the bootstrap standard errors are small, resulting in half-median unbiased estimates that are only slightly more narrow than the 90% confidence sets.

Unsurprisingly, the plug-in bounds on $\mu_{ate}$ shrink as the strength of our assumptions increase. The most flexible model is (ref) under assumption (A1), in which case the length of the bound on $\mu_{ate}$ is one.\footnote{Note that in a potential outcome framework with a binary treatment and binary outcome (and no other additional assumptions), the worst-case bound on the average treatment effect always has a length of one (c.f. manski1990nonparametric). } The identified set for $\mu_{ate}$ also always overlaps zero for model (ref). Results in Table (ref) suggest that full independence of $Z_3$ (A5 and A6) produces more informative bounds for $\mu_{ate}$ than our alternate conditional independence assumption (A7 and A8). In fact, our alternate conditional independence assumption does not provide much identifying power (compare the results under Assumptions (A3) and (A7)). On the other hand, full independence of $Z_3$ does induce a noticeable narrowing of the identified set for $\mu_{ate}$ (compare the results under Assumptions (A3) and (A5)). The results for this model are a useful benchmark to compare with cases where we impose linearity on the index function.

Next, we see in Table (ref) that the linear models from (ref) and (ref) narrow the bounds relative to the case of the general index function under some of the assumptions. Unsurprisingly, the smallest interval for $\mu_{ate}$ for model (ref) is obtained under Assumption (A6), in which case the sign of $\mu_{ate}$ is identified. For models (ref) and (ref) we make use of our method for profiling $\theta$, as described in Section (ref). In model (ref) we must profile on $\theta \in \mathbb{R}^2$ and there are 8 representative points. Interestingly, we find that under Assumptions (A1) - (A4) and (A7) - (A8), the identified set for $\theta$ is the entire Euclidean space $\mathbb{R}^2$. This illustrates that non-trivial bounds on $\mu_{ate}$ are possible even when the structural parameters are unidentified. Figure (ref) in Appendix (ref) shows the intervals computed using the linear programs of the form (ref) and (ref) for each representative point of $\theta$ under our various assumptions.

In the second linear model (ref), all coefficients are fixed. Thus, we need to profile on a parameter vector $\theta\in \mathbb{R}^3$. Our profiling procedure from Section (ref) returns $96$ representative points, each associated with a polyhedral cone in $\mathbb{R}^3$. Under Assumptions (A1) and (A2), the identified set for $\theta$ is $\mathbb{R}^3$, while for all other assumptions (A3) - (A8) we get a more informative identified set for $\theta$. In Figure (ref) in Appendix (ref) we also show the intervals computed using the linear programs of the form (ref) and (ref) for each representative point of $\theta$ under our various assumptions. Interestingly, the ATE bounds under (A1) for model (M2) and (M3) are the same as those under model (M1), which suggests that the functional form restrictions only become informative when combined with other modelling assumptions. The sign of $\mu_{ate}$ is identified for model (ref) under assumptions (A4) - (A8). The narrowest bounds for $\mu_{ate}$ under model (ref) are obtained under Assumption (A6), where the $90\%$ confidence interval is $[0.12,0.29]$. For comparison, using the same data but a slightly different model, acerenza2021testing also bound the ATE and obtain a $95\%$ confidence interval of $[0.02,0.29]$.\footnote{In addition to using a different model, acerenza2021testing also use the number of employees (without discretization) as their instrument, and they use inference procedure of chernozhukov2013intersection combined with a sample splitting procedure, which is valid under a very different set of assumptions than those presented in the current paper. }

Next we consider the counterfactual choice probability $\mu_{ccp}(0)$. Table (ref) reports the convex hull of the estimated identified set for $\mu_{ccp}(0)$ under various model specifications and under various assumptions. Similar to the bounds for $\mu_{ate}$, the half-median unbiased estimates are only slightly more narrow than the 90% confidence sets. We also see that the bounds on counterfactual choice probabilities tend to be wide and uninformative for most assumptions. Note that under Assumption (A1) to (A4) we always obtain the interval $[0,1]$ for the estimated identified set. The narrowest bounds are found in model (ref) under Assumptions (A5) and (A6). These bounds allow us to conclude that the probability an individual visits a doctor when provided private health insurance, given that they have no private health insurance and did not visit a doctor increases noticeably, with a magnitude in the interval $[0.37, 0.55]$ and a 90% confidence set of $[0.22, 0.62]$.

For the sake of comparison, we estimate the following bivariate probit:

where $(Z_1, Z_2, Z_3)$ are assumed to be independent from $(\varepsilon_1, \varepsilon_2)$, which are bivariate normal with mean zero, unit variance and correlation $\rho$. This model was estimated with our data using maximum likelihood, and $\mu_{ate}$ was estimated as $0.16$ with a bootstrapped confidence interval of $[0.11, 0.20]$. This value for $\mu_{ate}$ is not in the plug-in bounds under (M3) and (A5) and (A6), but it does lie within all of the 90% confidence sets in Table (ref), and seems to suggest strong evidence of a positive causal effect of health insurance on the decision to visit the doctor.\footnote{han2019estimation obtain a similar result in a model allowing for $\varepsilon_{1}$ and $\varepsilon_{2}$ to have unrestricted marginals, and a flexible dependence structure. However, they consider a different model from us, and the average treatment effect in han2019estimation is different from ours; we consider the average treatment effect averaged over all values of $(w,z)$, while they report the average treatment effect at the average value of their conditioning variables. They also report the average treatment effect at various quantiles of their conditioning variables.} However, the bivariate probit model is highly parameterized, and the results from Table (ref) suggest that under weaker assumptions the sign of $\mu_{ate}$ may not be identified.\footnote{In fact, acerenza2021testing reject the assumptions of the bivariate probit model using the same data but with a slightly different specification.}

The previous literature studying the effects of health insurance on the utilization of health care services is full of mixed results, and Table (ref) suggests that highly parameterized models may give highly significant, but possibly misleading results relative to models that make weaker assumptions.

Conclusion

This paper considers (partial) identification of a variety of counterfactual parameters in binary response models with possibly endogenous regressors. Importantly, our class of models allows for nonseparability of the index function in latent variables, and does not require any parametric distributional assumptions. Our specific partition of the latent variable space is key to our procedure, and we show how to enumerate the sets in this partition using results from the literature on computational geometry and hyperplane arrangements. In doing so, we provide a feasible method of constructing bounds on counterfactual quantities under a variety of different assumptions with multi-dimensional and nonseparable latent variables. We also thoroughly study the special case when the index function is linear in parameters, and show how to compute exact (i.e. not approximate) sharp bounds on counterfactual quantities. We also show how to adapt a recent inference procedure to the setting in this paper in order to construct confidence sets and bias-corrected estimates of the identified set. Finally, we show how to impose independence and monotonicity assumptions, and we present an application of our method to study the effects of private health insurance on the utilization of health care services.

The consideration of multinomial choice models, triangular systems, or general simultaneous discrete choice models (e.g. games, network formation, or models of social interactions) are all natural future extensions of the framework presented here which we intend to pursue. In addition, this paper emphasizes computational issues that arise in models that are partially identified. We believe exploring applications of state-of-the-art algorithms in computer science to problems in econometrics---as we have attempted here---is a fruitful avenue of future research.