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.
113,062 characters · 15 sections · 74 citation commands
On Quantile Treatment Effects, Rank Similarity, and Variation of Instrumental Variables
{ Rewriting in progress: please excuse the debris!}
This paper develops a nonparametric framework to identify and estimate distributional treatment effects under nonseparable endogeneity. We begin by revisiting the widely adopted rank similarity (RS) assumption and characterizing it by the relationship it imposes between observed and counterfactual potential outcome distributions. The characterization highlights the restrictiveness of RS, motivating a weaker identifying condition. Under this alternative, we construct identifying bounds on the distributional treatment effects of interest through a linear semi-infinite programming (SILP) formulation. Our identification strategy also clarifies how richer exogenous instrument variation, such as multi-valued or multiple instruments, can further tighten these bounds. Finally, exploiting the SILP's saddle-point structure and Karush-Kuhn-Tucker (KKT) conditions, we establish large-sample properties for the empirical SILP: consistency and asymptotic distribution results for the estimated bounds and associated solutions.
In the treatment effect literature, RS (or in its simplest form Rank Invariance) is widely used and plays a central role in identifying treatment effects in the presence of self-selection into treatment; see, e.g., foresi1995conditional, heckman1997making, chesher2003identification,Che05, chernozhukov2005iv, VY07, jun2011tighter, SV11, d2015identification, torgovitsky2015identification, vuong2017counterfactual, han2021identification, among many others. Specifically, RS requires that, within any given compliance group, the population ranks of the treated and untreated potential outcomes have the same distribution. The plausibility of this assumption is often questionable in applications and lacks firm theoretical justification (see, e.g., maasoumi2019gender, and also ChernozhukovHansen2013ARE for a survey). In response, several papers propose tests for the validity of RS; see frandsen2018testing, dong2018testing, and kim2022testing.
In this paper, we characterize RS by the model restrictions it imposes on the relationship between observed and counterfactual outcome distributions. Under a mild additional condition, we show that RS is equivalent to a two-sided preservation of first-order stochastic dominance (FOSD) across treatment statuses: for any two groups formed as linear mixtures of compliance types, if one group's treated-outcome distribution is FOSD-dominated by the other's, then the same dominance holds for their untreated counterfactual-outcome distributions; and vice versa. This connection between RS and FOSD preservation across potential-outcome distributions appears to have received limited attention in the literature. By establishing this equivalence, we bridge the structural IV formulation of chernozhukov2005iv with the potential-outcomes framework of rubin1974estimating.
Furthermore, we propose a relaxation of RS by introducing a weaker condition based on one-way preservation of FOSD. This new condition can be interpreted as requiring that the potential outcome without treatment is Rank Noisier (RN) than the treated potential outcome. Our proposed RN notion is related to, but distinct from, the Noisier concept developed in pomatto2020stochastic. We further show that this condition arises naturally in a broad class of structural models derived from primitive economic assumptions and motivated by empirical applications. Under this weaker condition, we derive bounds on distributional treatment effects via an SILP formulation.
Nonparametric identification of treatment effects with essential endogeneity has long been a challenging goal. For instance, manski1990nonparametric, manski1997monotone, MP00, among others, construct sharp bounds on the average treatment effect (ATE) under assumptions on the direction of treatment effects and on selection, while allowing instruments to be invalid in specific ways. Even with valid instruments, however, bounds on the ATE are typically wide and uninformative for precise policy prediction. The local ATE (LATE) imbens1994identification and local QTE abadie2002instrumental have been popular alternatives when we impose a monotonicity assumption on selection to treatment. However, the local group for which the treatment effect is identified may not be the group of policy interest. Therefore, the extrapolation of the local parameters is crucial for policy analysis (e.g., treatment allocation), in which case the identification challenge still remains mogstad2018using,han2020sharp.
Building on our identification strategy, we further solve empirical SILP problems to construct estimators for our bounds and establish their large-sample properties. We reformulate the SILP as a saddle-point problem ChristensenConnault2023. Following, e.g., shapiro1991asymptotic,ekeland1999convex, we show that the optimal-value map is Hadamard directionally differentiable with respect to the function-valued SILP coefficients to be estimated from the data. This establishes consistency and an asymptotic distribution for the bound estimators; See also GoffMbakop2025 for related results in finite-dimensional linear programming problems.
In our SILP, the optimal-value functional is, in general, only Hadamard directionally (rather than fully) differentiable with respect to the SILP coefficients. As emphasized by fang2019inference, this typically leads to non-Gaussian limiting distributions and can invalidate standard bootstrap procedures for inference. To obtain the full differentiability, we impose a local stability condition requiring uniqueness of the primal-dual optimal pair bonnans2000perturbation. Under this condition, the set of binding constraints forms a smooth manifold near the optimum, ensuring that both the primal and dual solutions vary smoothly under small perturbations of the SILP coefficients.
In addition, we exploit the SILP's saddle-point structure and KKT conditions to establish large sample properties for the optimal solution. Without assuming a full rank condition on the set of active constraints, we restructure the problem as a two-step nested optimization: an inner SILP that solves for a subvector of the solution while holding the other components fixed, and an outer problem that minimizes a convex criterion given by the inner problem's value function. Applying milgrom2002envelope's generalized envelope theorem, we obtain first- and second-order conditions for the outer problem, which we use to derive the asymptotic distribution for its solution. Moreover, exploiting the KKT conditions for the inner SILP, we obtain the limiting distribution of the remaining components. A key insight to this derivation is to combine complementary slackness with primal feasibility, providing us first-order conditions from the slackness equalities, since the slack attains its local minimum (or maximum) value at the active constraints. Our analysis contributes to the growing literature on inference for optimization-based estimators that confronts nondifferentiability and nonstandard asymptotics; see, e.g., andrews1999estimation,chernozhukov2007estimation,fang2019inference,hsieh2022inference,ChristensenConnault2023,fang2023inference,GoffMbakop2025, among many others.
The paper is organized as follows. Section (ref) introduces the treatment effect model, along with the necessary notation and assumptions. To motivate our approach, we begin by characterizing the RS assumption and then introduce a new condition that relaxes RS. We then derive bounds for the distributional treatment effects of interest under this weaker condition. Section (ref) discusses the computation of these bounds by solving an SILP problem, and Section (ref) establishes the asymptotic properties of the estimated bounds derived from the empirical SILP problem. Section (ref) presents numerical studies to illustrate our method. All proofs are provided in Appendix (ref).
This section presents a treatment effect model that accounts for individual self-selection into treatment, while allowing treatment effects to vary across individuals with identical observed characteristics (i.e. covariates). To fix ideas, and also in line with both the econometrics and statistics literature, we adopt Rubin's potential outcome framework heckman1986alternative,imbens1994identification, augmented with structural assumptions that address treatment selection. We then motivate our key identification assumption and use it to derive bounds on the distributional treatment effects of interest.
Let $Y_1$ denote the counterfactual outcome under treatment and $Y_0$ the counterfactual outcome without treatment. For simplicity, we assume both $Y_1$ and $Y_0$ are continuously distributed. The observed outcome $Y \in \mathcal{Y} \subseteq \mathbb{R}$ is defined as $Y \equiv Y_D$, where $D \in \{0, 1\}$ is the observed treatment indicator, representing an individual's treatment decision in response to instrumental variables (IVs) $Z$ and other covariates $X \in\mathcal X\subseteq \mathbb{R}^{d_X}$. We assume that $Z$ is either a vector of binary IVs or a multi-valued IV, taking $L$ distinct values, i.e., $Z \in \mathcal{Z} \equiv \{z_1, \dots, z_L\}$. Multi-valued or multiple IVs are common in both observational studies (e.g., multiple natural experiments affecting the same policy) and experimental studies (e.g., randomized control trials with multiple treatment arms implemented simultaneously or sequentially).\footnote{See mogstad2021causal for a recent survey.}
Following chernozhukov2005iv, each latent outcome $Y_d$ is expressed in terms of its quantile function. Specifically, for $d \in \{0,1\}$, \[ Y_d = q(d, X, U_d), \] where $q(d, x, \cdot)$ is the {\it quantile function} of $Y_d$, given $X=x$. Because $Y_d$ is continuously distributed, $U_d \sim U(0,1)$ and $q(d,x,\cdot)$ is strictly monotone. Note that $U_d$ represents the rank of $Y_d$ within its population cumulative distribution, i.e., $U_d=F_{Y_d|X}(Y_d|X)$. Moreover, the treatment selection $D$ is given by \[ D = h(Z, X, \eta), \] where function $h$ describes treatment selection by the individual, and $\eta\in\mathcal T\subseteq \mathbb R^{d_\eta}$ is a random vector capturing unobserved heterogeneity. As noted by imbens1994identification, individuals can be classified into distinct compliance types based on how their treatment choices would respond to each possible value of the IVs, as determined by $\eta$ (with covariates $X$ fixed).
In our analysis, we focus on evaluating the treatment effects on the distributional features of potential outcomes $Y_d$ for the treated (or untreated) populations. For $d\in\{0,1\}$ and $x\in\mathcal{X}$, the {\it quantile treatment effects} (QTE) at $\tau\in(0,1)$ for the group with treatment status $D=d$ and covariates $X=x$ is define as:
for $\tau\in(0,1)$.\footnote{$Q_{A| B}(\tau| b)\equiv \inf\{a\in\mathbb{R}: F_{A| B}(a| b)\ge \tau\}$ for (generic) random variable $A$ and random element $B$.} Similarly, the {\it (conditional) average treatment effects} (ATE) for the population with the treatment status $D=d$ and $X=x$ is given by:
Note that the unconditional QTE/ATE can be derived when these parameters are identified across all $d\in\{0,1\}$ and $x\in\mathcal{X}$. These parameters are of significant interest to researchers and policymakers for assessing the impact of interventions in policy evaluations. Because our model is fully nonparametric, we fix $X=x$ throughout this section.
Throughout, we maintain the assumption that the IVs are valid, i.e., the instrument $Z$ must be (strongly) relevant to individuals' treatment selection, and satisfy the exclusion restriction.
Assumption (ref) has been widely assumed in the literature; see e.g. imbens1994identification. Applying $Y_d$'s quantile expression, the first half of Assumption (ref) can be equivalently rewritten as $Z\perp (U_{d},\eta)\, | X$. In Assumption (ref), note that covariates $X$ are not required to be exogenous with respect to either the potential outcomes or the IVs.
To identify treatment effects under nonseparable endogeneity, a prevalent approach in the literature is to restrict heterogeneity in treatment effects via RS (or, in its simplest form, rank invariance).
By definition, RS ensures that the rank distribution of $Y_d$ remains unchanged across treatment status $d\in\{0,1\}$, regardless of the compliance type. As noted by chernozhukov2005iv, the strongest yet simplest form of the RS condition is Rank Invariance, which requires that $U_0=U_1$ almost surely.
Interestingly, we notice that RS is inherently linked to the preservation of first-order stochastic dominance (FOSD) between potential outcome distributions to their counterfactual counterparts, a connection that appears to have received limited attention in the existing literature.
Condition (ref) describes a biconditional relationship for linear mixtures of potential outcome distributions and their counterfactual counterparts, particularly ensuring the preservation of FOSD between potential outcome distributions. To see this, let \( W(\cdot) \) and \( \tilde{W}(\cdot) \) be arbitrary two probability distribution functions defined over \( \mathcal{T} \). Then, Condition (ref) implies that: the following FOSD relationship \[ \int_\mathcal T F_{Y_{1}|X,\eta}(\cdot|x,t)dW(t)\leq \int_\mathcal T F_{Y_{1}|X,\eta}(\cdot|x,t)d\tilde W(t) \] holds if and only if \[ \int_\mathcal T F_{Y_{0}|X,\eta}(\cdot|x,t)dW(t)\leq \int_\mathcal T F_{Y_{0}|X,\eta}(\cdot|x,t)d\tilde W(t). \]This result can be derived by defining \( G(\cdot) = W(\cdot) - \tilde{W}(\cdot) \) and setting \( c = 0 \) in Condition (ref).
The following lemma shows that under an additional “regularity” condition, RS can be equivalently characterized by Condition (ref).
\proof see Appendix (ref).\qed
Lemma (ref) shows that Condition (ref) is a key component of RS. To see this, note that RS can be equivalently represented as follows: For any absolutely continuous function \(G:\mathcal T\to\mathbb R\), \[ \int_{\mathcal T} F_{U_0| X,\eta}(\cdot| x,t)\,dG(t) = \int_{\mathcal T} F_{U_1| X,\eta}(\cdot| x,t)\,dG(t).\footnote{A simple sketch of the proof: If RS holds, integrating both sides with respect to $dG$ gives this representation. On the other hand, restricting $G$ to point-mass functions recovers the standard RS formula.} \] Thus, under Condition (ref), Lemma (ref) implies that the RS equivalence extends from the class \(\{W_u(\cdot): u\in(0,1)\}\) to all absolutely continuous functions \(G\).
In this subsection, we propose two alternative assumptions to relax RS. The first one corresponds to the “only if” part of Condition (ref), while the other corresponds to its “if” part. It is worth pointing out that we also drop the regularity condition in Lemma (ref) for our identification analysis. Grounded in economic theory, the proposed two assumptions have distinct empirical applicability; depending on the context, researchers may find that one, or both, applies.
Under Assumption (ref), we can further derive explicit model restrictions implied by Condition (ref). For any $(x, z) \in \mathcal{X} \times \mathcal{Z}$ and $d \in \{0,1\}$, note that \[ \Pr(Y_d \leq \cdot, D = 1 | Z = z, X = x) = \int_{\mathcal{T}} \Pr(Y_d \leq \cdot, D = 1 | \eta = t, Z = z, X = x) \, dF_{\eta | X,Z}(t | x, z). \] Since $D$ is binary and satisfies $D = h(Z, X, \eta)$, we have \[ \Pr(Y_d \leq \cdot, D = 1 | \eta = t, Z = z, X = x) = \Pr(Y_d \leq \cdot | \eta = t, Z = z, X = x)\times h(z,x,t), \]where $h(z,x,t)$ takes the binary value $0$ or $1$. By Assumption (ref), it follows that \[ \Pr(Y_d \leq \cdot, D = 1 | Z = z, X = x) = \int_{\mathcal{T}} F_{Y_d | X,\eta}(\cdot | x, t) \, dG(t |x, z), \] where $G(t |x, z) \equiv \int_{-\infty}^{t} h(z, x, u) \, dF_{\eta | X}(u | x)$. Thus, Condition (ref) implies the following model restriction: For any $\gamma_0\in\mathbb R$ and $ (\gamma_{11}, \dots, \gamma_{1L}) \in \mathbb{R}^{L}$, if
then
As will be discussed later, the above condition will be crucial for deriving bounds on distributional treatment effects.
Similarly, we define the converse of Condition (ref), denoted as Condition (ref).
Clearly, combining Conditions (ref) and (ref) together are equivalent to Condition (ref). It is of interest to understand when Condition (ref), Condition (ref), or both would hold in practice. To illustrate, we now derive e.g. Condition (ref) from primitive conditions on the joint distribution of $(U_0, U_1)$.
Following vuong2017counterfactual's counterfactual mapping approach, consider the mapping from \( F_{U_1| X, \eta}(\cdot | x, t) \) to \( F_{U_0 | X, \eta}(\cdot | x, t) \) for each \( t \in \mathcal{T} \). By the law of iterated expectations, we have \[ F_{U_0 | X, \eta}(\cdot\, | x, t) = \int_0^1 F_{U_0 | U_1, X, \eta}(\cdot \,| u_1, x, t) \, dF_{U_1 | X, \eta}(u_1 | x, t). \] This indicates that there exists a linear operator that maps the distribution \( F_{U_1| X, \eta}(\cdot \,| x, t) \) to \( F_{U_0 | X, \eta}(\cdot\, | x, t) \), with the conditional distribution \( F_{U_0 | U_1, X, \eta}(\cdot\, | \cdot, x, t) \) serving as the kernel function. Suppose in addition that \( U_0 \perp \eta \, | U_1, X \). Then, for any absolutely continuous function $G: \mathcal{T} \rightarrow \mathbb{R}$, \[ \int_\mathcal T F_{U_0 | X, \eta}(\cdot\, | x, t) \, d G(t)= \int_0^1 F_{U_0 | U_1, X}(\cdot\, | u_1, x) \, d\left\{\int_\mathcal TF_{U_1 | X, \eta}(u_1 | x, t)\, dG(t)\right\}, \] which is a Fredholm integral equation of the first kind. In this literature, significant attention has been given to understanding how the function \( \int_\mathcal T F_{U_0 | X, \eta}(\cdot | x, t) \, d G(t) \) can inform us about the behavior of its inverse function \( \int_{\mathcal T}F_{U_1 | X\eta}(\cdot |x, t)\, dG(t) \) under certain conditions on the kernel function $F_{U_0 | U_1, X}(\cdot |\cdot, x)$; See e.g. buchholz2021semiparametric.
In Assumption (ref), the first part implies that the rank \( U_1 \) captures all the information related to the compliance type. The second part introduces a weaker notion of positive dependence, i.e., positive regression dependence, than commonly used alternatives such as positive affiliation or decreasing inverse hazard rate Castro2007. In the following discussion, we will present structural examples in which Assumption (ref) is satisfied.
Let \( \mathcal{L}^1_+([0,1]) \) denote the set of all non-negative Lebesgue integrable functions on the interval \([0,1]\), i.e., \[ \mathcal{L}^1_+([0,1]) = \left\{ g: [0,1] \to \mathbb{R}_+ \,\Big|\, \int_0^1 g(t)\, dt < \infty \right\}. \] Let further \[ \mathcal{W}(x) \equiv \left\{ -\frac{\partial F_{U_0 | U_1, X}(u_0 \,| u_1, x)}{\partial u_1} : u_0 \in [0,1] \right\} \] be the collection of functions in $u_1\in[0,1]$, derived from $F_{U_0 | U_1, X}(u_0 \,| u_1, x)$.
In Lemma (ref), the first part shows that Condition (ref) follows directly from Assumption (ref). This suggests tha Condition (ref) can be justified on economic grounds. Similarly, one can consider the counterpart to Assumption (ref) by switching the roles of \( U_0 \) and \( U_1 \), so that Condition (ref) holds under an alternative economic interpretation. The second part of the lemma is a topological condition. Specifically, Condition (ii) requires that the collection \( \mathcal{W}(x) \), derived from the kernel function \( F_{U_0 | U_1, X}(\cdot \,| \cdot, x) \), is sufficiently rich to approximate any function in \( \mathcal{L}^1_+([0,1]) \) through linear combinations. This condition can be satisfied in some settings, e.g., \( F_{U_0 | U_1, X}(\cdot| u_1, x) = \mathbf{1}(\cdot \leq u_1) \) for any $u_1\in[0,1]$. In general, however, it is difficult to motivate for this condition or to verify it empirically.
Given Condition (ref) or Condition (ref), we are now prepared to construct bounds for distributional treatment effects in policy analysis. In the following discussion, we maintain Condition (ref) (respectively, Condition (ref)) and consider the quantile treatment effect \( QTE_{\tau}(1,x) \) (respectively, \( QTE_{\tau}(0,x) \)) at quantile \( \tau \in (0,1) \). By definition, note that \[ QTE_{\tau}(1, x) = Q_{Y_1 | D, X}(\tau | 1, x) - Q_{Y_0 | D, X}(\tau | 1, x), \] where \( Q_{Y_1 | D, X}(\tau | 1, x) \) is trivially identified by the “observed” conditional quantile \( Q_{Y | D, X}(\tau | 1, x) \). Thus, identifying \( QTE_{\tau}(1, x) \) reduces to identifying the counterfactual quantile \( Q_{Y_0 | D, X}(\cdot\, | 1, x) \), which requires to invert the conditional distribution \( F_{Y_0 | D, X}(\cdot \,| 1, x) \). For simplicity, we focus on constructing the lower bound of \( Q_{Y_0 | D, X}(\tau | 1, x) \) in the following discussion. It is straightforward to extend our approach analogously to the upper bound.
To proceed, we apply Condition (ref) - (ref) to a particular class of weights \( (\gamma_0, \gamma_{11}, \ldots, \gamma_{1L}) \in \mathbb{R}^{L+1} \) satisfying \( \sum_{\ell=1}^L \gamma_{1\ell} = 0 \), and then obtain the following result.
Because there may exist multiple (indeed, infinitely many) vectors \( (\gamma_0, \gamma_{11}, \ldots, \gamma_{1L}) \) satisfying \( \sum_{\ell=1}^L \gamma_\ell = 0 \) and eq. (ref), we aim to further tighten the bounds. To this end, we first introduce some notation to simplify the exposition. For each \( (z, x) \in \mathcal{Z} \times \mathcal{X} \), let \( p(z, x) = \Pr(D = 1 | Z = z, X = x) \) be the propensity score. Moreover, for \( d \in \{0,1\} \), define \[ \Delta_{dz}(\cdot \, | x) \equiv \Pr(Y \leq \cdot, D = d \, | Z = z, X = x) - \Pr(Y \leq \cdot, D = d \, | Z = z_L, X = x), \] supported on \( [\underline{y}, \overline{y}] \). In the above definition, we treat $z_L$ as the reference point of $Z$ on the support $\mathcal Z$. Next, define the \((L-1)\)-dimensional profile vector of functions \[ \Delta_d(\cdot \, | x) \equiv \big( \Delta_{dz_1}(\cdot \, | x), \ldots, \Delta_{dz_{L-1}}(\cdot \, | x) \big) \in \mathbb{R}^{L-1}, \] which captures the variation in \( \Delta_{dz}(\cdot | x) \) induced by the instrument. Finally, we denote \( \gamma_1 \equiv (\gamma_{11}, \ldots, \gamma_{1,L-1})' \in \mathbb{R}^{L-1} \) as a column vector. Note that by letting \( \gamma_{1L} = -\sum_{\ell=1}^{L-1} \gamma_{1\ell} \), the constraint \( \sum_{\ell=1}^L \gamma_{1\ell} = 0 \) is automatically satisfied.
The proof of Corollary (ref) is straightforward and thus omitted. Together, Theorem (ref) and Corollary (ref) highlight the identifying power of multi-valued IVs. In particular, in the above LP problem for the upper bound \( F^{UB}_{Y_0 | D, X}(\cdot | 1, x)\), its feasible region depends critically on the variations of the IVs.
To see this, consider the linear subspace of a Hilbert space (e.g., \( L^2([\underline{y}, \, \overline{y}]) \)) spanned by the set of functions \( \{1\} \cup \{ \Delta_{dz_\ell}(\cdot | x) : \ell \leq L \} \). Suppose further that each \( \Delta_{dz_\ell}(\cdot | x) \) belongs to the logistic family of distribution functions, and let \( L \to \infty \). Then the resulting linear subspace becomes dense in the space of distribution functions on \( [\underline{y}, \overline{y}] \), allowing for arbitrarily accurate approximation of $F_{Y | D, X}(y | 1, x) $; See, e.g., hornik1989multilayer. In this sense, \( F_{Y_0 | D, X}(\cdot | 1, x) \) becomes point identified in the limit by \( F^{UB}_{Y_0 | D, X}(\cdot | 1, x) \). On the other hand, when \( L \) is small, even an extra variation of $\mathcal Z$ expands the dimensionality of the linear subspace spanned by \( \{1\} \cup \{ \Delta_{dz_\ell}(\cdot | x) : \ell \leq L \} \), potentially tightening the bounds on \( F_{Y_0 | D, X}(\cdot | 1, x) \) substantially. In Section (ref), we investigate identification power of the IVs' variation by using Monte Carlo studies.
Finally, we construct a lower bound for $QTE_{\tau}(1,x)$ from \( F^{UB}_{Y_0|D,X}(\cdot\, | 1, x) \). Using the worst case bounds for the conditional quantile manski1994selection,blundell2007changes, we have
where $Q_{Y_{0}|D,X}^{LB}(\tau|1,x)$ is the $\tau$-th quantile of $F_{Y_{0}|D,X}^{UB}(\cdot\, |1,x)$. Because $F_{Y_{0}|D,X}^{UB}(\cdot\, |1,x)$ is point-wisely constructed, it may not be monotone over its support. Next, we define $Q^{LB}_{Y_{0}|DX}(\tau|1,x)$ by inverting the upper bound function $F_{Y_{0}|D,X}^{UB}( \cdot\, |1,x)$ as follows: \[ Q_{Y_{0}|D,X}^{LB}(\tau|1,x)= \inf_{ y\in\mathcal Y}\, \left\{ y: F_{Y_{0}|D,X}^{UB}( y\, |1,x)\geq \tau\right\}. \] Furthermore, we construct an upper bound for $QTE_{\tau}(1, x)$ by \[ QTE^{UB}_{\tau}(1, x) =Q_{Y_1 | D, X}(\tau | 1, x) - Q^{LB}_{Y_0 | D, X}(\tau | 1, x). \]
In this subsection, we first present an illustrative example to clarify the intuition behind our identification strategy. Next, we provide empirical justifications for the proposed key conditions, i.e., Conditions (ref) or (ref), by deriving them from primitive structural assumptions.
We begin with a setting involving multiple-valued IVs within imbens1994identification's LATE framework. Throughout, we fix \( X = x \). Assume that \( h(z_{\ell}, x, \cdot) \leq h(z_{\ell+1}, x, \cdot) \) for \( \ell = 1, \ldots, L \), and let \( D_z \equiv h(z, x, \eta) \) denote the potential treatment status when \( Z = z \) and \( X = x \). This setup leads to the following monotone selection condition:
Under this condition, the observed treated group \( \{D = 1\} \) comprises both always-takers (AT) and compliers (C). Define the always-takers (AT) as \( \{ \eta \in \mathcal{T} : D_1 = \cdots = D_L = 1 \} \), i.e., individuals who always receive treatment regardless of the IV value. For \( 1 \leq \ell \leq L - 1 \), define the complier group induced by shifting the instrument from \( z_\ell \) to \( z_{\ell+1} \) as: \[ (z_\ell, z_{\ell+1})_C \equiv \left\{ \eta \in \mathcal{T} : D_{z_j} = 0 \text{ if and only if } j \leq \ell \right\}. \] For instance, \( (z_1,z_2)_C \) denotes eager compliers, and \( (z_{L-1}, z_L)_C \) denotes reluctant compliers, following the terminology of mogstad2021causal. Within this LATE framework, the following lemma restates Theorem (ref) under the monotone selection assumption.
In Lemma (ref), the inequalities directly follow Conditions (ref)-(ref) under the monotone selection assumption, so the proof is omitted. According to Lemma (ref), if \( \sum_{k=1}^L \gamma_{1k} = 0 \), then Condition (ref) yields testable model implications, since \( \Pr[Y_d \leq \cdot \,, \eta \in (z_\ell, z_{\ell+1})_C] \) is identified for both \( d = 0 \) and \( d = 1 \); See imbens1997estimating. An exception arises when \( L = 2 \), as the condition holds trivially. When \( \sum_{k=1}^L \gamma_{1k} \neq 0 \), the model yields bounds on the counterfactual distribution \( \Pr(Y_0 \leq \cdot \,; \eta \in \text{AT}) \), and increasing \( L \) would tighten these bounds. Such extrapolation is central to Theorem (ref).
Next, we examine Condition (ref) under primitive structural assumptions. By Lemma (ref), Condition (ref) holds under Assumption (ref), which has an intuitive interpretation: Rank $U_0$ is {\it noisier} than $U_1$.
Our RN concept extends the “noisier” concept introduced by pomatto2020stochastic, where a random variable \(Z'\) is considered noisier than \(Z\) if \(Z' = Z + W\) and \(W\) is independent of \(Z\). Hence, if $X + Z\prec_{FOSD} \tilde X + Z$ for some $Z$ that is independent of $X$ and $\tilde X$, then $X +Z' \prec_{FOSD}\tilde X + Z'$ for any independent $Z'$ that is {\it noisier} than $Z$. In contrast, our RN concept supports a different type of preservation of stochastic dominance. Namely, for absolutely continuous function $G$ and $\tilde G$, if \[ \int_\mathcal T F_{Y_{1}|X,\eta}(\cdot\,| x,t)dG(t)\prec_{FOSD} \int_\mathcal T F_{Y_{1}|X,\eta}(\cdot\,| x,t)d\tilde G(t), \] then \[ \int_\mathcal T F_{Y_{0}|X,\eta}(\cdot\,| x,t)dG(t)\prec_{FOSD} \int_\mathcal T F_{Y_{0}|X,\eta}(\cdot\,| x,t)d\tilde G(t). \] Note that the reverse does not necessary hold for the same intuition provided in pomatto2020stochastic.
We now illustrate that the RN condition, namely, that \( Y_0 \) is rank-noisier than \( Y_1 \) or vice versa, can be justified within a class of empirical structural models.
In Example (ref), the assumption that \( \xi_0 \) is noisier than \( \xi_1 \) implies that \( \xi_0 = \xi_1 + \zeta \) for some independent noise \( \zeta \). It follows that \( V_0 = V + \xi_0 = V_1 + \zeta \). Thus, conditional on \( X \), the rank of \( V_0 \) is more dispersed than the rank of \( V_1 \), consistent with the RN condition. Importantly, this specific structure for the idiosyncratic shocks is not essential. For instance, one could instead impose a nonlinear form as follows: \[ \xi_0 = \max\{\phi(\xi_1), \zeta\}, \] where \( \phi \) is a strictly increasing function and \( \zeta \) is again independent.
In both Examples (ref) and (ref), Condition (ref) holds under the assumption that \( U_0 \) is rank-noisier than \( U_1 \). In this case, Theorem (ref) and Corollary (ref) yield bounds on \( QTE_{\tau}(1, x) \), the quantile treatment effect for those who receive the treatment. The next example, however, illustrates the converse case.
Example (ref) provides justification for Condition (ref), under which bounds on \( QTE_{\tau}(0, x) \), i.e., the quantile treatment effect for individuals who abstain from treatment, can be derived through a similar argument. Assuming either Condition (ref) or (ref), our approach yields partial identification of distributional treatment effects for the treated or the untreated, respectively.
In practice, it is often of interest to policymakers to determine which treatment parameter is most relevant for policy evaluation. Suppose the policymaker is primarily concerned with risk-averse individuals, who tend to prefer treatment options associated with lower uncertainty. For such a policymaker, she would be interested in evaluating the policy that would offer a form of “insurance”, i.e., either in the literal sense (e.g., health insurance) or through interventions that mitigate risk (e.g., vaccination, regulations on high-risk treatments). In this sense, our procedure offers a statistical tool for bounding the treatment effects for individuals with \( D = d \), when rank \( U_{d} \) is less noisier than \( U_{1-d} \).
In this section, we employ optimization methods to systematically calculate the bounds introduced in Corollary (ref). For clarity, we focus on the upper bound under the assumption that the IVs $Z$ are discrete, i.e., $Z \in \mathcal{Z} \equiv \{z_1, \dots, z_L\}$. Extending the analysis to the case of continuous $Z$ is more challenging and left for future work. Since our model specification is fully nonparametric, we condition on $X = x$ throughout this section, regardless of whether $X$ is continuously or discretely distributed.
To begin with, we introduce some notation. For expositional simplicity, we assume that the domain $\mathcal{Y}$ is a finite interval on $\mathbb{R}$, i.e., $\mathcal{Y} = [\underline{y},\, \overline{y}]$. For $d \in \{0,1\}$ and $z \in \mathcal{Z}$, recall that \[ \Delta_{dz}(\cdot| x) =\Pr(Y \leq \cdot, D = d | Z = z, X = x) - \Pr(Y \leq \cdot, D = d | Z = z_L, X = x) \] and \[ \Delta_{d}(\cdot | x)=\big( \Delta_{dz_1}(\cdot| x),\cdots,\Delta_{dz_{L-1}}(\cdot| x)\big) \] supported on $[\underline y,\, \overline y]$. By construction, $\Delta_{dz}(-\infty|x) = 0$ and $\Delta_{dz}(+\infty|z) = (-1)^{d+1}[p(z, x) - p(z_L, x)]$. Moreover, $\Delta_{dz}(\cdot|x) $ is bounded and differentiable.
We now consider the LP problem introduced in Corollary (ref) for the upper bound of \( F_{Y_0 | D, X}(y_0|1, x) \), where \( y_0 \in [\underline{y}, \overline{y}] \). Let \( \mathbb{S}(x) \) denote the feasible region of the LP, defined as \[ \mathbb{S}(x) = \left\{ (\gamma_0,\gamma'_1) \in \mathbb{R}^L : F_{Y | D, X}(y | 1, x) \leq \gamma_0+\gamma_1' \Delta_1(y | x) \quad \forall y \in [\underline{y}, \,\overline{y}] \right\}. \] By definition, \( \mathbb{S}(x) \) is constructed by a continuum of inequality constraints, with \( F_{Y| D X}(\cdot | 1, x) \) and \( \Delta_1(\cdot | x) \) as the intercept and slope functions, respectively. Moreover, it is straightforward that \( \mathbb{S}(x) \) is convex and unbounded. Thus, the upper-bound LP problem can be rewritten as: \[ F^{UB}_{Y_0 | D, X}(y_0 | 1, x) = \min_{(\gamma_0,\gamma'_1) \in \mathbb{S}(x)} \quad \gamma_0-\gamma_1' \Delta_0(y_0 | x). \] In this LP, any feasible move in the direction \( \big(-1,\,\Delta_{0}(y_0| x)\big)\) decreases the objective value.
In the empirical LP problem introduced later, it is essential that the corresponding population LP satisfies certain regularity conditions such that the optimal value is stable with respect to small perturbations in these coefficients.
Under Assumption (ref), the conditional distribution $F_{Y|DZX}(\cdot | 1, z, x)$ is continuous and differentiable a.e.. This assumption ensures that the functions $F_{Y|DX}(\cdot|1,x)$ and $\Delta_{1}(\cdot|x)$, which determine the LP constraints, vary continuously over the compact support $[\underline{y}, \overline{y}]$.
\noindentRemarks. Assumption (ref) strengthens the usual “no improving direction” (which only requires the infimum to be $\ge 0$) by imposing a strict margin $\varepsilon_0>0$. Geometrically, $(1,-\Delta_0'(y_0\,|\,x))$ strictly separates from the recession cone: along any feasible recession direction, the objective increases at a uniformly positive rate. Together with Slater's condition (introduced later), this rules out unbounded rays in the objective direction and ensures finiteness and robustness of the optimal value, and a perturbation-stable, bounded solution set. The larger the margin $\varepsilon_0$, the more robust the problem is to small perturbations.
With Assumptions (ref) and (ref), we now examine the well-posedness of the SILP problem and the stability of \( F^{UB}_{Y_0 | D, X}(y_0 | 1, x) \). As a first step, we verify that Slater's condition holds. This follows from the fact that, for any value \( \gamma_0 > 1 \), the vector \( (\gamma_0, 0, \dots, 0) \in \mathbb{R}^L \) lies in \( \mathbb{S}(x) \) as an interior point. By the Strong Duality Theorem, we obtain the following result. Its proof is standard in the SILP literature goberna1998linear,hettich1993 and is therefore omitted.
In Lemma (ref), strong duality (under Slater) justifies the min-max Lagrangian representation (i.e., a saddle point) of the upper bound used below and guarantees dual attainment for the ensuing sensitivity analysis. It should also be noted that the optimal solution $\gamma^*$ may not be unique. Throughout, we denote the set of optimal solutions by $\Gamma^*(x,y_0)$. Similarly, let $\lambda^*$ and $\Lambda^*(x,y_0)$ be an optimal solution to the dual problem and the collection of dual solutions, respectively. As is standard in the LP literature, one can use the interior-point algorithm to solve the above large-scale LP dual problem (i.e. high-dimensional $\lambda$).
In practice, when the solution set is unbounded, “small” but random perturbations in the LP's parameters might significantly affect the optimal objective value. To deal with such a sensitivity issue, we introduce an additional regularity condition. Let $\tau>0$ and consider the following modified LP problem, denoted as $LP(\tau)$:
where $\mathbb S(x,\tau)$ is defined as follows: \[ \mathbb S(x,\tau)=\big\{\gamma\in \mathbb S(x): \|\gamma\|^2_2\leq \tau \big\}. \] By definition, $\mathbb S(x,\tau)$ is a convex and compact subset of $\mathbb R^L$.
Note that the original LP can be written as $LP(+\infty)$. By Lemma (ref), there exists a finite optimal solution $\gamma^\ast$. Hence there is a threshold $\bar\tau<\infty$ such that, for all $\tau\ge \bar\tau$, the optimal value of $LP(\tau)$ remains the same as that of the original problem. The reason is that for some optimal solution of the unregularized LP the added regularity constraint is slack, so that solution remains feasible (and optimal) in $LP(\tau)$ and the objective value is unchanged.
Given our preceding discussion, the proof of Lemma (ref) is straightforward and therefore omitted. Throughout our following discussion, we maintain the assumption that $\tau \geq \bar{\tau}_x$.
The basic idea behind our estimation procedure is straightforward. In the \( LP(\tau) \) problem, the coefficients \( \Delta_{0z_\ell}(y_0 | x) \), \( \Delta_{1z_\ell}(\cdot | x) \), and \( F_{Y|DX}(\cdot | 1, x) \) are unknown but can be estimated nonparametrically from the data. This motivates the following two-step plug-in estimator. In the first step, we construct nonparametric estimates of these coefficients. In the second step, we solve the plug-in LP problem, denoted \( \widehat{LP}(\tau) \), under its dual formulation to obtain the estimated upper bound \( \hat{F}^{UB,\tau}_{Y_0|DX}(y_0 |1, x) \).
To clarify the idea, we consider an i.i.d. random sample \( \{(Y_i, D_i, X_i, Z_i) : i \leq n\} \) of size $n$, drawn from the population distribution of \( (Y, D, X, Z) \). For each \( (y, d, z, x) \in \mathcal{Y} \times \{0,1\} \times \mathcal{Z} \times \mathcal{X} \), we estimate the conditional distributions \( F_{Y|DX}(y | d, x) \) and \( F_{Y|DZX}(y | d, z, x) \) using kernel-smoothed empirical CDFs:
where \( K \) is a kernel function with compact support on $\mathbb R^{d_X}$, and \( h_0, h_1, h'_0, h'_1 \) are bandwidths associated with each conditioning group. For simplicity, we assume \( x \in \mathcal{X} \) is an interior point to avoid boundary issues in the above kernel estimates. However, note that the boundaries on \( \mathcal{Y} \) do not affect the consistency of the kernel estimates.
Furthermore, for each \( (y, d, z, x) \in \mathcal{Y} \times \{0,1\} \times \mathcal{Z} \times \mathcal{X} \), we estimate ${\Delta}_{dz}(y | x) $ by
in which $\hat p(z,x)$ is a kernel estimator of $p(z,x)$, i.e., \[ \hat{p}(z, x) = \frac{\sum_{i=1}^n \mathbbm{1}(D_i = 1, Z_i = z)\, K\left( \frac{x - X_i}{h^\dagger} \right)}{\sum_{i=1}^n \mathbbm{1}(Z_i = z)\, K\left( \frac{x - X_i}{h^\dagger} \right)}. \]where $h^\dag$ is also a bandwidth. Finally, let \[ \hat{\Delta}_d(y | x) = \big(\hat{\Delta}_{dz_1}(y | x), \ldots, \hat{\Delta}_{dz_{L-1}}(y | x) \big). \]
It is worth noting that the proposed estimators \( \hat{F}_{Y|DX}(\cdot | d, x) \) and \( \hat{\Delta}_{d}(\cdot | x) \) are not smooth over the support \( [\underline{y}, \overline{y}] \). In fact, they are step functions, as is typical for empirical distribution functions. As a result, \( \widehat{LP}(\tau) \) reduces to a finite-dimensional linear program, whose dimensionality increases with the sample size \( n \). Alternatively, one can construct smoothed estimators of the conditional distributions \( F_{Y|DX}(y | d, x) \) and \( F_{Y|DZX}(y | d, z, x) \) as follows:
where \( \Phi \) denotes the standard normal CDF and \( b_n > 0 \) with \( b_n \to 0 \) as \( n \to \infty \). With this modification, the estimated functions \( \tilde{F}_{Y|DX}(\cdot | d, x) \) and \( \tilde{F}_{Y|DZX}(\cdot | d, z, x) \) become differentiable. Interestingly, introducing such smoothness into the LP coefficients has little impact on the derivation of the asymptotic properties for $\hat{F}^{UB,\tau}_{Y_0|DX}(y_0|1,x)$, as will be discussed below.
Furthermore, let \( \hat{\gamma}^* \) denote an optimal solution to \( \widehat{LP}(\tau) \), and let \( \hat{F}^{UB,\tau}_{Y_0|DX}(y_0 | 1, x) \) be the corresponding optimal value, i.e.,
Again, the optimal solution \( \hat{\gamma}^* \) may not be unique. In this case, let $\hat \Gamma^*(x,y_0)$ be the set of optimal solutions to the empirical LP problem \( \widehat{LP}(\tau) \).
We now introduce some conditions to the above nonparametric estimators. For notation simplicity, let $\xi( x, y_0) \equiv \big( \Delta_0(y_0 | x),\, \Delta_1(\cdot | x),\, F_{Y | D, X}(\cdot | 1, x) \big)$ be the triple of vector and functional inputs to the SILP. For any compact subset \( A \subset \mathbb{R}^d \), define \[ \ell^\infty(A) \equiv \left\{ g : A \to \mathbb{R} \;\text{such that}\; \sup_{a \in A} |g(a)| < \infty \right\} \] as the space of bounded, real-valued functions on \( A \), equipped with the supremum norm. Moreover, let $\mathbb S_{\xi}\equiv \mathbb{R}^{L-1} \times \ell^\infty([\underline{y},\, \overline{y}])^{L-1} \times \ell^\infty([\underline{y},\, \overline{y}])$ be the support of $\xi(x,y_0)$. Let further $\hat \xi( x, y_0) \equiv \big( \hat\Delta_0(y_0 | x),\, \hat\Delta_1(\cdot | x),\, \hat F_{Y | D, X}(\cdot | 1, x) \big)\in\mathbb S_\xi$ denote the estimator of $\xi(x,y_0)$.
While Assumption (ref) is high level, it can be derived from standard primitive assumptions on, e.g., kernel functions and bandwidths in the nonparametric estimation literature; see, e.g., paganullah1999.
Assumption (ref) imposes a functional central limit theorem on the first-stage estimators. Since weak convergence in $\ell^\infty$ entails stochastic equicontinuity, we obtain, for any $\epsilon_n \downarrow 0$,
where $\hat{\mathbb G}_{1}(\cdot|x) = n^{\kappa}\big[\hat \Delta_1(\cdot | x) - \Delta_1(\cdot | x)\big]$ and $\hat{\mathbb G}_F(\cdot|x) = n^{\kappa}\big[\hat F_{Y| D,X}(\cdot | 1,x) - F_{Y| D,X}(\cdot |1,x)\big]$. These properties are standard in the empirical process literature and can be established under primitive conditions such as smoothness of the underlying distributions and regularity of kernel or series estimators; see, e.g., vandervaart1996weak and hardle2004nonparametric.
To illustrate (i) how having multiple instrument values sharpens our bounds and (ii) our two-step SILP estimation procedure, we conduct a focused Monte Carlo exercise. We consider the following DGP that satisfies Conditions (ref) but violates (ref). We generate a random sample $\{(Y_i,D_i,Z_i): i\leq n\}$ of $(Y,D,Z)$ generated as follows:
with $(\pi_0,\pi_1)=(0.2,0.5)$. For $d=0,1$, let further \[ U_d = U + \xi_d,\quad \ \xi_0=\xi_1+\nu, \] with (i) $(\nu,\,\xi_1) \bot\, (U,\eta)$; (ii) $\nu\bot\, \xi_1$; (iii) both $\nu$ and $\xi_1$ conform to standard normal distributions; (iv) $(U,\eta)$ conform to joint normal with mean $(0,0)'$ and variance $\Sigma=[1,0.8; 0.8,1]$. Moreover, let $Z\bot \, (\nu, \xi_1,U, \eta)$ and $Z \sim \text{Bin}(L-1,p)/(L-1) \;\in [0,1]$ for $L\in \{2,3,4\}$ and $p=0.5$. Here, $Z$ is normalized so that the endpoints of the support are invariant regardless of the value of $L$. This is intended to understand the role of the number of values Z takes while fixing the role of instrument strength fixed.
Our object of interest is the distribution of the potential untreated outcome among the treated, $F_{Y_0 | D}(\cdot \;| 1)$, and we compute identified upper and lower bounds as defined in Corollary (ref). Figure (ref) displays the true $F_{Y_0 | D}(\cdot \;| 1)$ together with the estimated lower and upper bounds for $L \in \{2,3,4,5\}$. All curves are obtained from simulations with a large sample size ($n=10{,}000{,}000$). The tuning parameter is fixed at $\tau=100$, but the regularization constraint $\|\gamma\| \le \tau$ is slack at the optimum and thus does not affect the SILP solutions (with a few exceptions due to sampling error). The upper and lower bounds are computed pointwise for each $y_0 \in [-6,6]$. As $L$ increases, i.e., as the IV offers richer support, the bounds tighten.
Fix $L=2$. We now examine the finite sample performance of the proposed SILP estimator. Specifically, we set $N=1000,2000$ and $4000$. For each $N$, we run $R=200$ replications to estimate the mean of the upper/lower bound and then construct 95% CI. To ensure Assumption (ref) hold for a sufficient large $\varepsilon_0$, we focus on $y_0\in[-0.520,\,1.688]$, where the two endpoints correspond to the 25th and 75th quantiles of $F_{Y_0| D}(\cdot\,| 1)$, respectively. For $y_0\notin[-0.520,\,1.688]$, note that the objective direction is close to the boundary of the dual cone of the feasible region's recession directions. Accordingly, rather than solving the SILP problems, we take the set of optimal solutions computed on $y_0\in[-0.520,\,1.688]$ and optimize the objective function over that set rather than over the SILP's feasible region. The resulting minimum and maximum objective values are reported as the upper and lower bounds outside the highlighted interval $[-0.520,\,1.688]$, respectively.
Figures (ref) and (ref) display the means and pointwise 95% CIs of the SILP bounds based on $R=200$ replications. Within the highlighted interval $[-0.520,\,1.688]$, we find: (i) the true bound lies within the confidence band; (ii) the mean of the lower bound estimator perfectly matches the true lower bound, while the mean of the upper bound estimator converges to the true upper bound as the sample size increases\footnote{For the upper bound problem, Assumption (ref) holds in a weak sense for small values of $y_0$.}; and (iii) all the dispersions of CIs decrease with sample size.
Section (ref) introduces a nonparametric estimator \( \hat{F}^{UB,\tau}_{Y_0 \mid D, X}(\cdot | 1, x) \) for the upper bound \( {F}^{UB,\tau}_{Y_0 \mid D, X}(\cdot | 1, x) \) by solving an empirical SILP. In this section, we study the asymptotic properties of this estimator. In addition, we also consider the large sample behavior of the associated empirical optimizer $\hat\gamma^*$ when the solution is unique. Throughout, our analysis is pointwise in the covariate value, i.e., fixing \( X = x \). Extending the theory to process convergence over the support of \(X\) is left for future work.
The asymptotic analysis of empirical SILP estimators has been studied by establishing Hadamard directional differentiability of the value map with respect to small perturbations of the SILP coefficients; see, e.g., shapiro1991asymptotic and bonnans2000perturbation. Under Slater's condition and compactness (in the weak-* topology) of the primal and dual optimal sets, $\Gamma^*(x,y_0)\subseteq\mathbb{R}^L$ and $\Lambda^*(x,y_0)\subseteq\Lambda$, they show that, for convex programs (including SILPs), the value function is Hadamard directionally differentiable, and that its derivative admits a min-max representation over the associated optimal sets.\footnote{Note that $\Lambda=\mathcal P([\underline y,\overline y])$ is defined as the set of Borel probability measures on the compact interval $[\underline y,\overline y]$, endowed with the topology of weak convergence (i.e., $\lambda_n\Rightarrow\lambda$ iff $\int f\,d\lambda_n\to\int f\,d\lambda$ for all $f\in\mathcal C([\underline y,\overline y])$). Because $[\underline y,\overline y]$ is a compact metric space, $\Lambda$ is compact in this topology (by Prokhorov's theorem). Note that the dual problem imposes linear moment constraints, which are weakly closed, hence the dual-feasible set is a closed subset of $\Lambda$ and therefore compact.} Following their approach, we treat our upper-bound estimator as a min-max object and begin by establishing its consistency.
For fixed $(x,y_0)$, define the Lagrangian \[ \mathcal{L}\big(\xi(x,y_0);\gamma,\lambda\big) \equiv \gamma_0 - \gamma_1' \Delta_0(y_0| x) + \int_{\underline y}^{\overline y}\!\big[F_{Y| D,X}(y| 1,x) - \gamma_0 - \gamma_1' \Delta_1(y| x)\big]\,d\lambda(y), \] with $\gamma\in\mathcal B_\tau\equiv\{\gamma\in\mathbb R^L:\|\gamma\|\le\tau\}$ and $\lambda\in\Lambda$. Define the value functional $\phi:\mathbb S_\xi\to\mathbb R$ by \[ \phi(\xi)\;\equiv\;\min_{\gamma\in\mathcal B_\tau}\ \max_{\lambda\in\Lambda}\ \mathcal L(\xi;\gamma,\lambda). \] By Lemma (ref) (strong duality) and Sion's minimax theorem, since $\mathcal B_\tau$ is compact, $\Lambda$ is weakly compact billingsley1999convergence, and $\mathcal L$ is jointly continuous, convex in $\gamma$, and concave in $\lambda$, sion1958general's minimax theorem gives \[ F^{UB,\tau}_{Y_0| D,X}(y_0| 1,x) =\min_{\gamma\in\mathcal B_\tau}\max_{\lambda\in\Lambda}\mathcal L\big(\xi(x,y_0);\gamma,\lambda\big) =\max_{\lambda\in\Lambda}\min_{\gamma\in\mathcal B_\tau}\mathcal L\big(\xi(x,y_0);\gamma,\lambda\big). \] See, e.g., rockafellar1970convex. It implies that any $(\gamma^*,\lambda^*)\in\Gamma^*\times\Lambda^*$ is a saddle point satisfying \[ \mathcal L(\xi(x,y_0);\gamma^*,\lambda)\ \le\ \mathcal L(\xi(x,y_0);\gamma^*,\lambda^*)\ \le\ \mathcal L(\xi(x,y_0);\gamma,\lambda^*)\quad \forall(\gamma,\lambda)\in\mathcal B_\tau\times\Lambda, \] which we use to establish continuity of $\phi$ (hence consistency) and, under additional regularity, its Hadamard directional derivative.
We next derive the limiting distribution of $\hat{F}^{UB,\tau}_{Y_0|D,X}(y_0\,|\,1,x)$. By the generalized envelope theorem of milgrom2002envelope, the value functional $\phi$ is Hadamard directionally differentiable at $\xi(x,y_0)$. The derivative takes the form of a support-function-type expression over the primal-dual optimal sets $(\Gamma^*,\Lambda^*)$ and is therefore, in general, only piecewise linear in the perturbation direction.
Under Assumption (ref), note that the empirical process $n^{\kappa}[\hat\xi(x,y_0)-\xi(x,y_0)]$ converges in distribution to $\mathbb G_\xi(\cdot|x,y_0)=(\mathbb G_0(x,y_0),\mathbb G_1(\cdot|x),\mathbb G_F(\cdot|x))$, a tight Gaussian element in the space $\mathbb S_{\xi}$. For each $(\gamma,\lambda)\in\mathcal B_\tau\times\Lambda$, define \[ \mathbb Z(\gamma,\lambda|x,y_0)\equiv -\gamma'_1\mathbb G_0(x,y_0) \;+\;\int_{\underline y}^{\overline y}\!\Big[\mathbb G_F(y| x)-\gamma_1' \mathbb G_1(y| x)\Big]\,d\lambda(y). \] Because $(\gamma,\lambda)\mapsto \mathbb Z(\gamma,\lambda| x,y_0)$ is a continuous linear functional of $\mathbb G_\xi(\cdot|x,y_0)$ and the index set $\mathcal B_\tau\times\Lambda$ is compact, the process $\{\mathbb Z(\gamma,\lambda| x,y_0):(\gamma,\lambda)\in\mathcal B_\tau\times\Lambda\}$ is a tight Gaussian process.
If the optimal set $\Gamma^*$ or $\Lambda^*$ is not a singleton, the limiting distribution is a min-max functional of a Gaussian process, reflecting the piecewise-linear geometry of the value functional $\phi$. In particular, $\phi$ is only Hadamard directionally differentiable at $\xi(x,y_0)$, its derivative map is generally nonlinear, rather than fully (linearly) Hadamard differentiable. As emphasized by fang2019inference, this lack of full differentiability can render the standard bootstrap invalid. In such cases one typically employs nonstandard methods for inference (e.g., the derivative-based bootstrap in fang2019inference and the numerical delta method in hong2018numerical).
In this subsection, we establish asymptotic normality of the SILP's optimal value and, under uniqueness and nondegeneracy, of the optimal solution. Specifically, we impose additional regularity ensuring local stability of the primal and dual solutions, so that the solution mapping is locally unique and Lipschitz in $\xi(x,y_0)$. Under these additional conditions, the value functional $\phi$ is fully Hadamard differentiable at $\xi(x,y_0)$ (i.e., it admits a continuous linear derivative), in contrast to the merely directional differentiability that holds without such stability.
In the LP sensitivity literature, local stability means that under small perturbations of the coefficients $\xi(x,y_0)$, the set of binding constraints and their multipliers remain stable, and the optimal primal and dual solutions vary smoothly (e.g., locally Lipschitz). A standard sufficient condition is the Linear Independence Constraint Qualification (LICQ) bonnans2000perturbation, which, in our SILP with optimaization vector $\gamma\in\mathbb R^L$, requires exactly $L$ linearly independent active constraints at the optimum $\gamma^*$. This can be overly restrictive and lacks clear economic motivation in our application. We therefore impose a weaker regularity condition that still guarantees local stability while allowing the number of active constraints to be equal or less than $L$.
In Assumption (ref), linear independence of the active constraints, together with strictly positive dual multipliers (a form of strict complementarity), ensures uniqueness of the dual solution. To see this, note that the first-order KKT conditions gives us\footnote{Note that the constraint $\|\gamma\|^2\leq \tau$ is inactive.} \[ \Delta_{0}(y_0| x)+\sum_{k=1}^K \lambda_k^*\,\Delta_{1}(y^*_{1k}| x)=0. \] Because the \(K\) active constraints (i.e. linear equations with slopes \(\{[\,1;\,\Delta_1(y^*_{1k}| x)\,]\}_{k=1}^K\)) are linearly independent, the above linear system admits a unique solution for the strictly positive dual variables \(\{\lambda_k^*\}_{k=1}^K\). W.l.o.g., we assume that the first \( K \) components of the coefficient vectors in the active constraints are linearly independent, i.e., \[ Rank
=K. \]
Under Assumption (ref), the asymptotic normality of the estimated upper bound follows directly from Theorem (ref). Let \((\gamma^*,\lambda^*)\) be the unique optimal primal-dual solution pair. Then the limiting distribution in Theorem (ref) is normal, i.e., \[ n^{\kappa}\Big[\hat{F}^{UB,\tau}_{Y_0| D,X}(y_0| 1,x)-F^{UB,\tau}_{Y_0| D,X}(y_0| 1,x)\Big] \ \overset{d}{\rightarrow}\ \mathbb{Z}(\gamma^*,\lambda^*\,|x,y_0). \] This shows that the way first-stage estimation error propagates into sampling variation of the bound estimator is governed by the optimal primal vector \(\gamma_1^*\) and the active-constraint multipliers \(\{\lambda_k^*\}_{k=1}^K\).
Next, we move to the asymptotic properties for the estimator $\hat\gamma^*$ as a solution to the empirical SILP problem. Recall that the optimal solution $\gamma^*$ identifies the mixture of complier subgroups: their potential treated outcome $Y_1$ FOSD the probability distribution of treated outcome.
With a unique primal-dual solution under Assumption (ref), we establish consistency of \(\hat\gamma^*\) by adapting the Argmin Theorem vandervaart1998asymptotic to the saddle-point (min-max) structure: For each \(\gamma\in\mathcal B_\tau\), denote
For notation simplicity, let $\xi=\xi(x,y_0)$ and $\hat \xi=\hat \xi(x,y_0)$ in the following discussion. By Assumption (ref), we obtain
where all the three inequalities use the saddle point structure. By uniqueness of \(\gamma^*\) and compactness of \(\mathcal B_\tau\), the standard Argmin Theorem gives \(\hat\gamma\overset{p}{\to}\gamma^*\). By a similar argument, we can also show that $\hat\lambda^*\overset{p}{\rightarrow}\lambda^*$ in the space $\Lambda$ (which is compact in the weak-* topology).
Furthermore, we consider the limiting distribution of $\hat\gamma^*$. Because there are \(K\) active constraints (\(K\le L\)) at the optimum, this suggests that there might be less binding constraints at $\hat\gamma^*$ than the SILP's dimensionality $L$. To deal with this issue, we restructure the SILP problem via a two-step nested optimization. First, we reparametrize the optimization vector \(\gamma\in\mathbb R^L\) as \(\gamma=(\theta_1,\theta_2)\in\mathbb R^K\times\mathbb R^{L-K}\), where \( \theta_1=(\gamma_0,\gamma_{11},\dots,\gamma_{1,K-1})'\in\mathbb R^K\) and \(\theta_2=(\gamma_{1K},\dots,\gamma_{1,L-1})'\in\mathbb R^{L-K}\). For each $y\in[\underline y, \overline y]$ and $d\in\{0,1\}$, let further $\Psi_{d1}(y| x)=\big[(-1)^{d+1},\ \Delta_{d1}(y| x),\dots,\Delta_{d,K-1}(y| x)\big]' \in\mathbb R^{K}$, and $\Psi_{d2}(y| x)=\big[\Delta_{dK}(y| x),\dots,\Delta_{d,L-1}(y| x)\big]' \in\mathbb R^{L-K}$. Moroever, for each fixed \(\theta_2\), we define an inner SILP problem, denoted as $\widetilde {LP}(\theta_2,\tau)$, for solving the \(K\)-dimensional parameter \(\theta_1\) as follows: \[
\] where \({\mathcal B}_{1\tau}=\{\theta_1\in\mathbb R^{K}:\|\theta_1\|^2\leq \tau\}\). In addition, let \(\theta_1^\dagger(\theta_2)\) be the solution to the inner problem. Under Assumption (ref), $\theta_1^\dagger(\theta_2)$ is well defined for $\theta_2$ belongs to a neighborhood of $\theta_2^*$, i.e., there exists a unique solution to the above inner SILP problem, when $\theta_2\in\mathcal N_\epsilon(\theta^*_2)$ for some $\epsilon>0$. By definition, $Q(\cdot|x,y_0)$ is convex.\footnote{To see this, let $\rho\in[0,1]$, and $\theta_2,\tilde \theta_2\in \tilde B_\tau$. Note that $\theta_{1,\rho}\equiv \rho\theta_1^\dag(\theta_2)+(1-\rho)\theta^\dag(\tilde \theta_2)$ belongs to the feasible region of $\widetilde {LP}(\rho\theta_2+(1-\rho)\tilde \theta_2,\tau)$, and \[ \rho Q(\theta_2|x,y_0)+(1-\rho) Q(\tilde \theta_2|x,y_0)= \theta'_{1,\rho}\Psi_{01}(y_0|x)\geq Q(\rho\theta_2+(1-\rho)\tilde \theta_2|x,y_0). \]}
In the outer step, we minimize the criterion \(Q(\cdot| x,y_0)\) over the compact set \( B_{2\tau}=\{\theta_2\in\mathbb R^{L-K}:\ \|\theta_2\|^2\le \tau\}\), i.e. \[ \min_{\theta_2\in B_{2\tau}}\ Q(\theta_2| x,y_0). \] By Assumption (ref), the unique solution \(\gamma^*=(\theta_1^*,\theta_2^*)\) to the population SILP is recovered by the two-step procedure: \(\theta_2^*\) solves the outer problem and \(\theta_1^*=\theta_1^\dag(\theta_2^*)\) solves the inner SILP. Appendix B establishes local stability of the inner SILP at its optimum, covering \(\theta_1^\dag(\theta_2)\), the active set, and the associated multipliers for \(\theta_2\) in a neighborhood \(\mathcal N_\varepsilon(\theta_2^*)\). In addition, we also provide an expression for the first and second order derivative of $Q(\cdot|x,y_0)$ by applying milgrom2002envelope's generalized envelope theorem and exploiting the KKT conditions.
Under local stability (unique optimizer, constant active set, smooth multipliers), the generalized envelope theorem of milgrom2002envelope gives a Hadamard directional derivative for the inner value, which in our setting is linear in \(\theta_2\) within a neighborhood of \(\theta_2^*\). Therefore, \(Q(\cdot| x,y_0)\) is continuously differentiable at \(\theta_2^*\). In particular, since \(\theta_2^*\) uniquely minimizes \(Q(\cdot| x,y_0)\), we have $\frac{\partial Q(\theta_2^*| x,y_0)}{\partial \theta_2}=0$, and the Hessian matrix \(\frac{\partial^2 Q(\theta_2^*| x,y_0)}{\partial \theta_2\partial \theta_2'}\) is positive semi-definite. In Appendix B, we derive the expressions for $\frac{\partial \hat Q(\theta_2|x,y_0)}{\partial \theta_2}$ and $\frac{\partial Q(\theta_2^*| x,y_0)}{\partial \theta_2\partial \theta_2'}$ by applying milgrom2002envelope's generalized envelope theorem.
Using Taylor expansion of the outer FOC gives us the stated representation for \(n^{\kappa}(\hat\theta_2-\theta_2^*)\). The proof is straightforward and therefore omitted.
Next, we derive the asymptotic distribution of $\hat\theta_1^*$ by exploiting the KKT conditions from the inner SILP. In particular, we characterize how the estimated binding points $\{\hat y^*_{1k}\}$ and the constraint residuals (slacks) respond to perturbations in the estimated coefficients $\hat\xi(x,y_0)$.
For $(\theta_1,\theta_2)\in\mathbb R^{K}\times\mathbb R^{L-K}$ and $y_1\in[\underline y,\overline y]$, define the residual functions, which is associated with the constraints, as follows: \[
\] Note that $\hat r$ is the plug-in estimator of $r$. At the primal optima, feasibility implies $r(\cdot\,|x,\theta_1^*,\theta^*_2)\leq 0$ and $\hat r(\cdot\,|x,\hat\theta_1^*,\hat\theta^*_2)\leq 0$; at binding points, we have $r(y_1^*\,|\,x,\theta_1^*,\theta^*_2)=0$ and $\hat r(\hat y_1^*\,|x,\hat\theta_1^*,\hat \theta^*_2)=0$.
By definition, it follows immediately that
These inequalities highlight the distance between the estimated and population residual functions at the active constraints. Moreover, because for any $y_1 \in [\underline y, \overline y]$, \[ \hat r(y_1 | x, \hat \theta_1,\hat \theta_2) = \hat r(y_1 | x, \theta^*_1, \theta^*_2) - (\hat\theta_{1}-\theta^*_{1})' \hat{\Psi}_{11}(y_1 | x) - (\hat\theta_2-\theta^*_2)' \hat \Psi_{12}(y_1 | x), \] therefore, substituting this expression into the inequalities above gives us:
By the equicontinuity assumed in Assumption (ref), we obtain the following result.
The proof is omitted, as it follows directly from the previous discussion.
For notational convenience, let $y_1^* \equiv (y_{11}^*,\ldots,y_{1K}^*)'\in [\underline y, \, \overline y]^K$ and define \[ {r}\big( y^*_{1} | x, \theta^*_1, \theta^*_2 \big)=[{r}\big( y^*_{11} | x, \theta^*_1, \theta^*_2 \big) ;\cdots; {r}\big( y^*_{1K} | x, \theta^*_1, \theta^*_2 \big)] \]as a $K$-dimensional column vector. Let further
By definition, $\mathbb A_{11}(y_1^*|x)\in\mathbb R^{K\times K}$ and $\mathbb A_{12}(y_1^*|x)\in\mathbb R^{(L-K)\times K}$ denote the Jacobian matrices of the active constraints w.r.t. $\theta_1$ and $\theta_2$, respectively. Note that we assume $\mathbb A_{11}(y_1^*|x)$ has full rank, hence is invertible.
TBA