Bounds on inequality with incomplete data

Bounds on inequality with incomplete data

\spacing{1}

\spacing{1.4}

Introduction

It is often important to study inequality in settings where the variable of interest is only partially observed. These include studies when researchers draw upon historical data tables in previously published sources, use data from organizations which group individual records for privacy reasons, or draw upon statistical surveys which allow respondents to report interval responses. In such environments, an inequality index (such as the Gini coefficient, among others) will be set identified in general. While point identification could be obtained under strong assumptions - such as imposing a functional form for the distribution of the variable of interest, or assuming homogeneity within data groups - we instead focus on characterizing and estimating this identified set. Specifically, we develop a unified framework that (i) delivers sharp, nonparametric bounds for a broad class of inequality indices under several observational regimes; (ii) provides computationally tractable algorithms to obtain those bounds; and (iii) gives a general inference procedure for them. These allow us to gain more credible understanding of the level of inequality in these settings, without placing unrealistic assumptions on the data.

Viewed abstractly, we study a setting in which the inequality functional $\mathcal{J}(\mu_0)$ is set identified. Here $Y$ is a latent scalar outcome with distribution $\mu_0$, $\mathcal{J}$ is an inequality functional, and the observed data reveal only a coarsened version of $Y$: for each unit, we observe an interval $[\underline{a}_i,\overline{a}_i]$ such that $Y_i \in [\underline{a}_i,\overline{a}_i]$. In addition, we allow for auxiliary information about the distribution of $Y$, such as overall or subgroup means, which can be encoded as linear (in)equality restrictions on the feasible distributions. The identified set of interest is the set of values taken by $\mathcal{J}(\mu)$, where $\mu$ ranges over all distributions of $Y$ consistent with the observed interval information and any auxiliary linear restrictions. Our goal is to characterize and estimate the sharp lower and upper bounds on this set through optimization of $\mathcal{J}$. Our procedures are simple to implement with off-the-shelf solvers, scale well with sample size, and can easily incorporate additional linear information when available. The bounds are sharp in the sense that every value within the bound interval can be attained by some distribution consistent with the observed information and auxiliary restrictions.

Section (ref) formalizes the setting we study and focuses on two observational regimes common in empirical work, defined by the nature of the available information. Scenario 1 covers the case where there are a collection of non-overlapping intervals within which each observation's value is known to lie. These intervals may be predetermined (e.g. survey brackets) or constructed ex post, either exogenously or on the basis of percentile ranges of the responses. Scenario 2, by contrast, features individual-specific intervals that may overlap across individuals, as is common when sensitive variables are elicited through unfolding bracket designs. We allow intervals to degenerate to point observations in both scenarios; the framework therefore accommodates data sets that combine exact and interval observations\footnote{This occurs for example in our empirical application of Scenario 2, where many respondents report exact wealth values while others report intervals.}. The key distinction is that Scenario 1 imposes a non-overlapping structure on the intervals, while Scenario 2 permits arbitrary overlap, leading to a different feasible set of latent distributions and, correspondingly, different sharp bounds for $\mathcal{J}(\mu_0)$.

Our key insight is that many inequality indices, such as the Gini coefficient, the 90/50 quantile ratio, and others (possibly with mild modifications), can be expressed as linear-fractional functions of the sorted vector of the variable of interest. This representation is the key computational observation: it turns an infinite-dimensional optimization over feasible distributions into a finite-dimensional program with a globally optimal solution. Section (ref) develops this idea in the context of Scenario 1 and shows that the problem of finding sharp bounds can be cast as a linear program in terms of the original variables (e.g., income or wealth), using well-known results in optimization theory developed by charnescooper. Additional linear inequality or equality constraints can be incorporated naturally, so the framework is organized around empirically distinct observation schemes rather than around special cases of auxiliary information. We also illustrate that, in the special case with no additional constraints, the Gini coefficient (which we consider as the most prominent index in the family of linear-fractional indices after sorting) can be rewritten as a quadratic-to-linear objective in terms of the unknown mass shares assigned to the interval boundaries. In this form, the problem can be solved using results from the fractional programming literature, such as dinkelbach67. We propose algorithms to implement this approach.

Section (ref) also provides solution-form results for linear-fractional indices (after sorting) and for general Schur-convex inequality indices within Scenario 1. These results identify conditions under which the optimizers of the inequality indices of interest are supported on either a fixed number of points or on a number of points that increases sufficiently slowly with the sample size, which is later used in the asymptotic analysis. For linear-fractional indices, a broad range of additional linear constraints can be incorporated. For Schur-convex indices, the admissible auxiliary restrictions are more structured, but still include the main forms of auxiliary information emphasized in the paper.

Section (ref) focuses on Scenario 2. We derive bounds on inequality indices in the presence of individual-specific interval data. We first characterize the forms of the optimizers for both linear-fractional and general Schur-convex measures, and then use these results to develop generic computational approaches for this scenario. These approaches rely on identifying mass shares allocated to a set of points known ex ante. The section also presents an alternative, direct method formulated in terms of the original observations for the Gini coefficient, and the same construction extends to other linear-fractional measures. We state the Scenario 2 results for the baseline case without additional distributional constraints because this is the empirically relevant case in our application; adding linear constraints amounts to augmenting the same finite-dimensional formulation. The informational environment of Scenario 2 is somewhat similar to that studied by ManskiTamer2002. However, their analysis focuses on partial identification of conditional mean functions in regression models, whereas our objective is to obtain sharp bounds for global distributional measures, which requires a different geometric, computational, and inferential treatment.

Section (ref) develops statistical inference by viewing the bound endpoints as value functions of a constrained optimization problem over probability measures. Under the regularity conditions stated there, this value function map is directionally differentiable, which yields $\sqrt{n}$ asymptotics for the estimated lower and upper bounds and bootstrap validity for inference. Although much of our computational discussion focuses on linear-fractional indices, the asymptotic argument also applies to several other inequality measures of interest, including the Theil, mean-log deviation, and Atkinson indices, when the conditions in Section (ref) are satisfied. In our formulation, the feasible set is characterized through collections of piecewise-affine linear functionals of the latent distribution, possibly involving estimated nuisance objects such as cutoffs or point-data shares. This abstraction allows the grouped-data and interval-data environments studied in the paper, as well as closely related variants, to be analyzed under a common framework.

Section (ref) presents a brief application for each scenario. Firstly, we illustrate the methodology by looking at wealth inequality using data from the English Longitudinal Study of Ageing (ELSA) which is characterised by having complex interval responses. This application naturally fits our Scenario 2. Our second application uses historical U.S.\ income distribution tables to construct time-series bounds for the Gini coefficient, selected quantile ratios, and the Hoover index under grouped information, with additionally available linear information taken into account where possible. Taken together, the results demonstrate the value of a flexible, transparent toolkit for studying inequality when incomes and wealth are only partially observed. As well as providing credible bounds for inequality in modern-day liquid savings wealth and the historical trends in US income inequality from the 1930s to the 1960s, these applications also demonstrate how the size of the resulting bounds depends on the quantity and complexity of the underlying interval data and the granularity of source material in the form of distributional information or additional data tables.

Related Literature. Within the inequality-measurement literature, related work studies partial identification of inequality or spread measures under coarse data. Stoye2010 provides elegant characterizations of the identification regions for the expectation and several spread parameters, including the Gini coefficient, when the underlying distribution is observed only through interval-censored or otherwise coarsened data. Our setting differs in several important respects. First, we allow for overlapping interval information of the type considered in Scenario 2, which Stoye2010 does not address. Second, we incorporate auxiliary linear restrictions arising from additional data sources. Third, Stoye2010 focuses exclusively on characterizing identification regions and does not provide computational or inferential methods for the grouped-data environment he considers, which is most closely related to our Scenario 1A. By contrast, we develop tractable algorithms and establish $\sqrt{n}$ inference for the bound endpoints via a directional delta method with bootstrap validity, making the bounds implementable in the settings we consider.

Other existing literature has studied selected special cases of Scenario 1, leading to valuable insights for bounding inequality measures. Many contributions, however, examine one particular configuration of auxiliary information at a time -- e.g., known overall mean, known subgroup means, or knowledge of specific points on the Lorenz curve, and are typically not designed to extend beyond the specific case they treat. Much of this literature constructs bounds by interpolating the Lorenz curve from a small set of known points: gastwirth1972 uses subgroup means and produces non-sharp bounds, while mehran1975 uses known Lorenz curve points to attain sharp bounds. murray1978 takes a different route, giving a computational approach to finding Gini bounds under a known overall mean.

These papers illustrate the usefulness of bounding approaches in a variety of grouped-data environments, but each is tailored to a specific piece of auxiliary information. By contrast, our framework treats auxiliary information generically through linear (in)equality constraints and remains applicable when subgroup means or Lorenz-curve points are unavailable. Even within our Scenario 1, this substantially enlarges the set of feasible data environments.

To our knowledge, existing work does not provide sharp bounds and formal inference for settings such as Scenario 2, where individuals are observed through overlapping intervals. Most existing methods are developed for environments in which groups partition the support of income or wealth and are therefore mutually exclusive and ex-ante ordered. These approaches do not directly extend to the informational environments considered here. When intervals overlap across individuals, the feasible set of latent distributions has a different geometric structure, and techniques based on Lorenz-curve interpolation or group-level aggregation are no longer directly applicable. As a result, this setting raises distinct conceptual and computational challenges.

Our paper also contributes to the sparse literature on statistical inference for inequality bounds. MCDONALD1981 provides inference in the setting of gastwirth1972 using their closed-form Gini bounds, while gastwirth_et_al_estimation derives joint asymptotic distributions for Gini bounds and extends the analysis to other measures. DedduwakumaraPrendergast2019 proposes a parametric bootstrap when subgroup means are available. Our approach remains applicable in all these cases and beyond, without requiring closed-form bounds or restrictive informational assumptions.\footnote{A complementary literature imposes parametric structure on the income distribution; see Jorda_at_al2021 for a survey. Our approach is distribution-free and can be used as a sensitivity benchmark for parametric extrapolations.}

Related methodological work includes Cowell1991, who derives solution forms for strictly Schur-convex inequality measures when either the overall mean or the subgroup means are known in our Scenario 1. In contrast, we allow for general (not necessarily strictly) Schur-convex measures and characterize solution forms under a broader class of linear (in)equality restrictions in Scenario 1. We also obtain solution-form results for general Schur-convex measures in Scenario 2.

Finally, our setting naturally connects to ecological inference and data combination, where individual-level distributions are learned from coarse micro data combined with aggregate information from other sources (see CrossManski2002). We contribute by providing a unified value-function approach that yields sharp, finite-sample-implementable bounds for widely used inequality indices and by developing bootstrap-valid inference for the set-identified parameters. Recent econometric work on partial identification under combined data includes Pacini2019, DH2024. We differ from this work in both objective and method by focusing on inequality indices and exploiting structural properties (such as Schur-convexity or linear-fractional representations).

Data Scenarios and Inequality Indices

We study settings where a variable of interest $Y$, such as income, wealth, consumption or wages, is not observed exactly but only through interval restrictions. Formally, for each unit $i=1,\ldots,n$ there is an unobserved realization $y_i$, and we observe an interval $\mathcal{I}_i=[\underline{a}_i,\overline{a}_i]$ known to contain it, i.e.\ $y_i\in\mathcal{I}_i$, where $\underline{a}_i \leq \overline{a}_i$ are (possibly estimated) finite interval limits. The collection $\mathcal{I}:=\bigl(\mathcal{I}_1,\ldots,\mathcal{I}_n\bigr) $ summarizes the observable information, while $\mathbf{y}=(y_1,\ldots,y_n)\top$ denotes the unobserved realization vector. We focus on two scenarios that cover many applications.

A rich set of auxiliary information can be incorporated either by refining the group structure or by adding linear constraints. For instance, if the sample median $Q_{0.5}(\mathbf{y})$ is known and lies in group $\mathcal{G}_{d_0}$, we can split $\mathcal{G}_{d_0}$ into $[\underline{a}_{d_0},\,Q_{0.5}(\mathbf{y})]$ and $\big(Q_{0.5}(\mathbf{y}),\,\overline{a}_{d_0}\big]$ and update the group counts accordingly. Thus, the knowledge of percentiles can be incorporated through the refinement of intervals $\mathcal{G}_d$ rather than in the form of additional constraints to take into account in optimization.

The additional linear-restriction setup covers many cases of interest, including known subgroup means or income-share restrictions. For example, knowledge of the overall sample mean $\widehat{\mu}$ can be imposed through $\bar{y}=\widehat{\mu}$, while knowledge that the mean in group $d$ equals $\widehat{\mu}_d$ adds the restriction $\bar{y}_d=\widehat{\mu}_d$. Likewise, if a point $\big(\sum_{d=1}^{h}\widehat{s}_d,\,\ell_h\big)$ on the Lorenz curve is known for some $h$, where $\ell_h := \bar{y}^{-1}\!\sum_{j=1}^{h}\widehat{s}_j \bar{y}_j$, then one can impose the linear restriction $\sum_{j=1}^{h}\widehat{s}_j \bar{y}_j - \ell_h\,\bar{y}=0$. While prior work typically treats such cases in isolation and rarely develops inference, our framework accommodates a broad class of auxiliary information that can be expressed through the linear equality/inequality restrictions studied in this paper. It also allows such information to be combined across multiple data sources within a unified computational and inferential framework.

A fixed number of non-overlapping groups featured in Scenario 1 is common in historical or privacy-binned data such as the historical income distribution tables reproduced in Figures A.1 and A.2 in the online supplement. We study both the simple case with no auxiliary information (1A), in which only ${\mathbf{1}(y_i \in \mathcal{G}_d), i=1, \ldots, n, d=1, \ldots, D}$ are observed, and the richer case (1B), where interval indicators are supplemented with linear constraints. We allow for the possibility that $\underline{a}_{d} = \overline{a}_d$, which leads to a setting with mixed point and interval data and non overlapping groups. In many applications, however, the data consist purely of intervals.

By contrast, Scenario 2 allows overlapping intervals, as in interval-response surveys that follow decision trees.

Intervals in this scenario are individual-specific and many of them will potentially be different. Like in Scenario 1, this allows for $\underline{a}_i= \overline{a}_i$ for some $i$, ultimately leading to scenarios of mixed point and interval data. An example of this is explored in our application using data from English Longitudinal Study of Ageing (ELSA, see Section (ref)). Auxiliary linear information can be added here by the same logic as in the passage from Scenario 1A to Scenario 1B, but we keep Scenario 2 in its baseline form because this is the case commonly encountered in interval-response survey data.

In both scenarios we characterize the identified set via a finite-dimensional optimization problem, which yields both fast computation and a useful statistical theory.

A key ingredient is that many common inequality indices can be written, after sorting, as linear-fractional functions of the outcome vector: \[ G_n(\mathbf{y})=\frac{r_1(n)^\top \mathbf{y}}{r_2(n)^\top \mathbf{y}}, \] where $r_1(n),r_2(n)\in\mathbb{R}^n$ are known vectors and $\mathbf{y}$ is understood to be sorted so that

To ensure these indices are well-defined, we assume throughout that denominators are strictly positive, i.e.\ $r_2(n)^\top \mathbf{y}>0$ for all feasible $\mathbf{y}$. Our two leading examples are the following:

Gini coefficient. Under (ref), this takes the form of $G_n(\mathbf{y}) = \frac{1}{n^2 \bar{y}} \sum_{i=1}^n (2i - n - 1)y_i$, where $\bar{y} = \frac{1}{n} \sum_{i=1}^n y_i$. This fits the above equation with $r_1(n) = (1-n, 3-n, \ldots, n-1)^\top$ and $r_2(n) = n \iota_n$, where for any $k\in\mathbb{N}$, $\iota_k$ denotes the $k$-vector of ones.

Quantile ratio. The sample quantile ratio $y_{\lceil \tau_2 n \rceil} / y_{\lceil \tau_1 n \rceil}$ for quantile indices $\tau_1, \tau_2$ (e.g., $\tau_1 = 0.5, \tau_2 = 0.9$) fits our setup by taking $r_1(n)=e_{\lceil \tau_2 n \rceil}$ and $r_2(n)=e_{\lceil \tau_1 n \rceil}$, where $e_j$ denotes the $j$th canonical basis vector in $\mathbb{R}^n$.

Other linear-fractional statistics include weighted or generalized Gini indices, the top-$p$ income share, the Palma ratio, percentile ratios, and the Bonferroni index. While our main computational results apply to the linear-fractional family, our characterization of the identified set extends to any continuous and Schur-convex inequality function. Beyond the linear-fractional examples, this includes the Generalized Entropy family (including the Theil and mean-log-deviation indices), as well as the Atkinson, Eltető–Frigyes, Kolm, and Zenga indices, and the Herfindahl–Hirschman index when applied to income shares. In certain special cases, such as the Hoover index, we can still leverage linear programming ideas to obtain fast computation even though the index itself is not linear-fractional (discussed for Hoover in more detail before Proposition (ref)). We require some mild assumptions on the support of data in order to work with certain indices: for example, log-based indices such as the Theil index require all \(y_i\) to be bounded away from \(0\), while quantile-based indices such as the 90-50 ratio require a separate regularity condition: a Lipschitz assumption on the relevant quantiles, together with a local density bound away from zero around those quantiles. All details are given in relevant sections and in the appendix.

Computation in Scenario 1

Because the feasible set is compact and connected, and $G_n(\cdot)$ is continuous, the image of the feasible set under $G_n$ is a compact connected subset of $\mathbb R$, hence a (possibly degenerate) closed interval. Therefore the identified set is fully characterized by its upper and lower bounds, obtained by minimizing and maximizing $G_n(\mathbf y)$ over the feasible set.

Charnes-Cooper transformation for optimizing linear-fractional inequality indices

To compute sharp nonparametric bounds for a linear-fractional $G_n(\mathbf{y})$, we solve

subject to (ref),

and letting $n_\ell$ denote the number of observations in group $\ell$:

Optimization of (ref) subject to (ref), (ref) and (ref) constitutes a linear-fractional problem, and it is well known how to convert it to a linear program by means of the charnescooper transformation. To illustrate how this can be accomplished, let us rewrite the constraints shaping the feasible set in the matrix form. The ordering (ref) can be written as $E_n \, \mathbf{y} \leq 0,$ where $E_n$ is the $(n-1)\times n$ matrix of first differences whose elements $(r,r)$, $r=1, \ldots, n-1$, are 1, elements $(r,r+1)$, $r=1, \ldots, n-1$, are -1, and all the other elements are 0. The constraints $y_i \in [\underline{a}_d, \overline{a}_d]$ can be written as $I_{n}\mathbf{y} \le b_{U,n}$ and $-I_{n}\mathbf{y} \le -b_{L,n}$, where \[ b_{U,n} := (\overline{a}_1 \iota_{n_1}^\top, \ldots, \overline{a}_D \iota_{n_D}^\top)^\top, \qquad b_{L,n} := (\underline{a}_1 \iota_{n_1}^\top, \ldots, \underline{a}_D \iota_{n_D}^\top)^\top, \] and $I_n$ is the identity matrix of size $n$. Overall, if we denote \[ H_n := [E_n^\top,\ I_n,\ -I_n,\ C_n^\top]^\top, \qquad b_n := (0_{(n-1)\times 1}^\top,\ b_{U,n}^\top,\ -b_{L,n}^\top,\ f_n^\top)^\top, \] then the constraint set can be written as $H_n \mathbf{y} \le b_n$.

Charnes-Cooper transformation. Our linear-fractional problem can be solved as the following linear program: $$\max_{z,t} \, r_1(n)^\top \, z \qquad \text{subject to } \qquad H_n z -b_n t \leq 0, \quad r_2(n)^\top z=1, \quad t > 0. $$ This reformulation of the linear-fractional program is obtained by the well-known Charnes-Cooper transformation with $z=\frac{1}{r_2(n)^\top\mathbf{y}} \cdot \mathbf{y}$, $t = \frac{1}{r_2(n)^\top\mathbf{y}}.$ The optimal solution $(z^*,t^*)$ for $(z,t)$ yields the solution of the original problem as $\mathbf{y}^*=\frac{z^*}{t^*}$. This reformulation is attractive because it reduces the problem to linear programming, and also naturally allows one to incorporate a variety of additional information through the constraints $C_n \mathbf{y} \leq f_n$. If $n$ is very large (say, in millions), one can utilize the population-size invariance property of an inequality metric to scale down the computational problem to a more feasible one, at the cost of a small approximation error.

Solution form for linear-fractional inequality indices

The form of the optimal solution $\mathbf{y}^*$ depends on the additional linear constraints $C_n\mathbf{y} = (\leq) f_n$. Propositions (ref) and (ref) characterize the solution. Proposition (ref) treats Scenario 1A, first for general linear-fractional measures and then for strictly Schur-convex ones.

The result of Proposition (ref) allows reformulation of the optimization problem over a $D$-dimensional parameter $\hat{\mathbf{p}} = (\hat{p}_1, \ldots, \hat{p}_D)^\top$, where $\hat{p}_d \in \{0, 1/n_d, \ldots, 1\}$ represents the proportion of units in group $d$ assigned to $\underline{a}_d$, with $1 - \hat{p}_d$ assigned to $\overline{a}_d$. The exact reformulation will depend on the vectors $r_1(n)$ and $r_2(n)$ in the linear-fractional definition. In the case of the Gini coefficient, we can reformulate the objective function as

where $\hat{s}_d = n_d / n$, $\hat{\mathbf{s}}=(\hat{s}_1, \ldots, \hat{s}_D)^\top$, $A(\hat{\mathbf{s}})$ is a $2D \times 2D$ symmetric matrix of differences between interval boundaries, and $b(\hat{\mathbf{s}})$ is a $2D \times 1$ vector of weighted boundaries:

This function needs to be optimized over $\mathcal{P}_n \equiv \times_{d=1}^D \{0,1/n_d,2/n_d,\ldots, 1-1/n_d,1\}$.

An analogous $G(\hat{\mathbf{p}}, \hat{\mathbf{s}})$ form can be constructed for any linear-fractional $G_n(\mathbf{y})=r_1(n)^\top\mathbf{y}/r_2(n)^\top\mathbf{y}$ for any linear-fractional index with the strict Schur-convex property; non-Schur indices (such as quantile ratios) are handled separately below. For the lower bound of the strict Schur-convex one, Proposition (ref) implies that a minimizer $\hat{\mathbf{p}}^*_{\min}$ takes the form $(0, \ldots, 0, 1, \ldots, 1)$, with a switch from 0 to 1 at some group $d_0$. The sharp lower bound is computed by evaluating $G(\hat{\mathbf{p}}, \hat{\mathbf{s}})$ at $D-1$ vectors $(0_{1 \times m}, \iota_{D-m}^\top)$, $m = 1, \ldots, D-1$, and selecting the minimum. For the upper bound, a maximizer $\hat{\mathbf{p}}^*_{\max}$ takes the form $(1, \ldots, 1, \hat{p}_{d_0}, 0, \ldots, 0)$ with $d_0\in \{1,\ldots,D\}$ and $\hat{p}_{d_0}$ potentially any in $[0,1]$.\footnote{The atomic nature of the argmax and argmin may suggest that imposing smoothness on the c.d.f. of the latent variable $Y$ could sharpen the bounds. However, smoothness alone is unlikely to help, since the step-function c.d.f.s corresponding to the argmax and argmin can be approximated arbitrarily closely by smooth c.d.f.s. Restrictions on the upper or lower bounds of the p.d.f. of $Y$, if it exists, could potentially tighten the bounds.}

The observation that our sharp inequality bounds can be taken as dependent only on a finite-dimensional parameter $\hat{\mathbf{p}}$\footnote{This is when $D$ is fixed. Our statistical theory will allow $D$ to grow with sample size.} will be a cornerstone of our statistical theory; analogous finite-dimensional reductions appear in the more general settings we consider.

We now outline a general approach to maximizing or minimizing the quadratic-to-linear objective in ((ref)). While minimization is unnecessary in Scenario 1A by the results of Proposition (ref), developing this framework prepares us for Scenario 2.

Consider the maximization of ((ref)). Following dinkelbach67, we introduce a family of subproblems indexed by $\lambda$: \[f_{\max,n}(\lambda)= \max_{\hat{\mathbf{p }} \in \mathcal{P}} \left(\frac{1}{2} (\hat{\mathbf{p}}^\top, \iota_D^\top-\hat{\mathbf{p}}^\top ) \; A(\hat{\mathbf{s}}) \; (\hat{\mathbf{p}}^\top, \iota_D^\top-\hat{\mathbf{p}}^\top )^\top - \lambda (\hat{\mathbf{p}}^\top, \iota_D^\top-\hat{\mathbf{p}}^\top) \; b(\hat{\mathbf{s}})\right)\] where $\mathcal{P} := [0,1]^D$, and we let $f_{\min,n}$ define the analogous $\min$ problem. By dinkelbach67, the function above is continuous, strictly decreasing on any interval where $B(\hat{\mathbf p}) := (\hat{\mathbf{p}}^\top, \iota_D^\top-\hat{\mathbf{p}}^\top) \; b(\hat{\mathbf{s}})>0$, and it is the supremum of affine functions of $\lambda$, hence convex. Under the sign conditions $f_{\max,n}(0)\ge 0$ and $f_{\max,n}(1)\le 0$ (which hold here since $G\in[0,1]$), there is a unique zero $\lambda^\star=\max_{\hat{\mathbf p}\in\mathcal{P}} G(\hat{\mathbf p},\hat{\mathbf s})$.

We solve for $\lambda^\star$ using bisection $[0,1]$: at each step evaluate $f_{\max,n}$ at the midpoint and update the bracket according to the sign; this is given in Algorithm (ref), which iteratively locates the solution with arbitrary tolerance $\varepsilon$ (e.g., $\varepsilon=10^{-6}$). This differs from the classical dinkelbach67 procedure (and also procedures in subsequent literature), and achieves geometric convergence since $|\lambda_{i+1}-\lambda_i| = 2^{-(i+1)}$, and $\lambda^*$ is always contained within whichever of the intervals $[\lambda_i, \lambda_i + 2^{-(i+1)}]$ and $[\lambda_i - 2^{-(i+1)}, \lambda_i]$ are selected in the $i$-th iteration of the algorithm. Note that the optimization in $f_{\max,n}$ is over the full set $\mathcal{P}$, not the grid $\mathcal{P}_n$. This greatly speeds computation, at the cost of an asymptotically negligible $O(1/n)$ approximation error as $\mathcal{P}_n$ becomes dense (see Lemma (ref)). Optimizing over $\mathcal{P}$ also allows one to work solely with sample shares in each interval, as the sample size is no longer needed to construct $\mathcal{P}_n$.

The lower bound $\min_{\hat{\mathbf{p}} \in \mathcal{P}} G(\hat{\mathbf{p}}, \hat{\mathbf{s}})$ can be computed by adapting Algorithm (ref): replace $f_{\max,n}(\widetilde{\lambda})$ with $f_{\min,n}(\widetilde{\lambda})$, keep the same sign update rule as in the maximization case (increase $\lambda$ if $f_{\min,n}(\widetilde{\lambda})>0$ and decrease $\lambda$ if $f_{\min,n}(\widetilde{\lambda})<0$), and use the stopping condition $0 \ge f_{\min,n}(\widetilde{\lambda}) \ge -\varepsilon$. This is redundant in Scenario 1A (Proposition (ref)) but may be needed more generally. Let us now discuss the computation of other linear-fractional inequality indices.

\vskip 0.05in

Quantile ratio. This inequality index is not strictly Schur-convex. In Scenario 1A, sharp bounds for the sample quantile ratio are immediate: the upper bound is $\overline{a}_{d_2}/\underline{a}_{d_1}$ and the lower bound is $\max\{1,\underline{a}_{d_2}/\overline{a}_{d_1}\}$, where $\sum_{d=1}^{d_j}\hat{s}_d\ge \tau_j$ but $\sum_{d=1}^{d_j-1}\hat{s}_d<\tau_j$ for $j=1,2$.\footnote{We interpret the upper bound as $+\infty$ if $\underline{a}_{d_1}=0$.} Hence, a maximizer (minimizer) $\mathbf{y}^*$ can be taken at the interval boundaries. Since only $y^*_{\lceil\tau_2 n\rceil}$ and $y^*_{\lceil\tau_1 n\rceil}$ matter for the objective, the optimal solution is not unique.

\vskip 0.05in

Hoover index. While the Hoover index is not linear-fractional, we can leverage linear programming ideas to ensure fast computation. For each $k=0,\ldots,n$, consider a linear-fractional objective $\frac{1}{2n\bar{y}}\left(\sum_{i=1}^k(\bar{y}-y_i) + \sum_{i=k+1}^n(y_i-\bar{y}) \right)$ subject to (ref), (ref) (in Scenario 1B also subject to (ref)) and additional constraints $y_k \leq \bar{y}$, $y_{k+1} \geq \bar{y}$ for $k<n$. For some such $k$ the additional constraints may result in an empty feasible set. But there will be $k$ for which it is non-empty. For all $k$ with non-empty feasible set optimize it using the methods outlined above. Among solutions for different $k$, choose the one that gives the maximum value of the objective and the one that gives the minimum.

Proposition (ref) characterizes the solution in Scenario 1B. With a fixed number of constraints, the optimal solution is supported on finitely many values (possibly sample-dependent) and the proposition bounds how the number of distinct values can grow as constraints are added. This result is important for our asymptotic theory.

General results for Schur-convex indices

Although our computational emphasis is on indices that are linear-fractional after sorting, we also present solution-form results for a general Schur-convex inequality index $G_n(\mathbf{y})$. This identifies the structure of sharp optimizers beyond the linear-fractional class while keeping the computational development focused on the measures used in our applications.

First, note that Theorem (ref) allows constraints on just some (not necessarily all) elements of group $d$. Second, note that conditions on matrix $\widetilde{C}_{n,d}$ allow for constraints on several subgroups within group $d$. Suppose $d=1$. Then we could have constraints (equality as well inequalities from above or/and below) e.g. on the sum of $y_1+...+y_{[n_1/2]-1}$, on the median value $y_{[n_1/2]}$ in that group and on the sum $y_{[n_1/2]+1}+\ldots+y_{n_1}$. A special case of the condition given in Theorem (ref) is when $\widetilde{C}_n$ is invariant to a permutation of columns corresponding to elements within each group $d$: that is, $\widetilde{C}_n$ is invariant to permutation of columns $\sum_{k=1}^{d-1} n_k+1, \ldots, \sum_{k=1}^{d} n_k$. From the practical point of view, the condition on matrix $\widetilde{C}_n$ in Theorem (ref) is consistent with many empirically relevant settings, as most real-world constraints will treat some (or all) variables within a given interval group in a symmetric way. All the examples of constraints we consider (such as subgroup means, income ratios, points on the Lorenz curve, overall means) satisfy this condition. This same symmetry is what makes the later measure formulation natural: once the objective and constraints depend only on blockwise masses or moments, permutations within a group are irrelevant, so the explicit ordering constraints used in the computational derivations become bookkeeping devices rather than substantive restrictions.

An important conclusion from Theorem (ref) is that asymptotically we can look at optimizers of $G_n(y)$ as those composed from a finite number of values, both in a sample and asymptotically (as long as the structures of $C^{(1)}_n$ and $C^{(2)}_n$ remain stable as not to lead to the increase of $k_d+o_d$ with $n$).

Computation in Scenario 2

Scenario 2 allows mixed point observations and individual-specific intervals that may overlap. Unlike Scenario 1, the data generally cannot be summarized by counts in a fixed set of non-overlapping brackets and therefore cannot be fully ordered ex ante (only partially ordered). As a result, the linear-fractional structure of the inequality measures of interest cannot be exploited or transformed into linear programs, since these measures typically admit known linear-fractional representations only under a full ordering.

We start by splitting individuals into two subsets: those that have point data on $y_i$ and those with genuinely interval data. Define the Point Set as \(P = \{i : \underline{a}_i = \overline{a}_i \}\), and the Interval Set as \(Q = \{i : \underline{a}_i < \overline{a}_i\}\). Let $\mathcal{B} := \{\underline{a}_i : i\in Q\}\cup \{\overline{a}_i : i\in Q\}$ and $\mathcal{U} := \mathcal{B}\cup \{\underline{a}_i : i\in P\}$. Thus $\mathcal{B}$ records boundary values coming only from the genuinely interval observations in $Q$, while $\mathcal{U}$ augments $\mathcal{B}$ with the realized point observations from $P$. Conceptually, point data are just degenerate intervals; we separate them only for computation, because their locations are already fixed whereas only observations in $Q$ generate unknown allocations.

Unlike Section (ref), we first present solution-form results, and only after that discuss computational aspects. In our solution structure results we first consider general Schur-convex inequality indices in Theorem (ref), and then Schur-convex indices with a linear-fractional representation (with a fixed ordering of $y_i$ in the sample) in Theorem (ref).

The main insight of Theorem (ref) is that at an extremum almost all interval observations can be placed at the interval boundary values in $\mathcal{B}$. This identifies the structure of sharp solutions in Scenario 2 beyond the linear-fractional class. Theorem (ref) refines the result of Theorem (ref) for linear-fractional indices and shows that at most one interior placement (for the max) and at most one interior value overall (for the min) - as described in Theorem (ref) - have to be in $\mathcal{U}$.

Theorem (ref) turns the structure of solutions into implementable programs as it allows us to look for a solution in the form of numbers $N_u$ allocated to each point $u \in \mathcal{U}$ when we consider linear-fractional inequality indices. Equivalently, point data are just degenerate intervals whose allocation is already known, so they enter as known masses at values in $\mathcal{U}$, whereas only observations in $Q$ generate unknown allocations and, hence, the ones that need to be constrained.

\paragraph*{First step: Ordering} Without loss of generality we can suppose that all unique elements of $\mathcal{B}$ are arranged in the increasing order \({b}_1 < {b}_2 < \cdots < {b}_{K}\), and all unique elements in $\mathcal{U}$ are ordered in an increasing way. Let $u(d)$ index the position of $b_d$ in ordered $\mathcal{U}$.

\paragraph*{Second step: Exhaustive inequality constraints} on $N_{u}$ and their various sums, where $N_u$ is the number of points from $Q$ found at particular $u \in \mathcal{U}$.

Take any consecutive block $[b_d, b_{d+k}]=\bigcup_{j=0}^{k-1}[b_{d+j},b_{d+j+1}] $ ($k=1,\dots,K-1$; $d=1,\dots,K-k$), The lower bound counts intervals $i \in Q$ fully contained in the block: \[ \sum_{i \in Q} \mathbf{1}(\mathcal{I}_i \subseteq [b_d,b_{d+k}]) \leq \sum_{u=u(d)}^{u(d+k)} N_u. \] The upper bound counts intervals $i \in Q$ that overlap the block (even partially): \[ \sum_{u=u(d)}^{u(d+k)} N_u \leq \sum_{i \in Q} \mathbf{1}(\mathcal{I}_i \cap [b_d,b_{d+k}] \neq \emptyset). \] The case $k=K-1$ yields the equality $\sum_{u=u(1)}^{u(K)} N_u = |Q|$.

Constraints for non-consecutive unions are implied by those for consecutive unions (lower bounds additive; upper bounds subadditive with gap adjustments), adding no new information. Indeed, for $d_1+1<d_2$, for the lower bound the constraint $$\sum_{u \in [b_{d_1}, b_{d_1+1}] \cup [b_{d_2}, b_{d_2+1}]} N_u \geq \sum_{i \in Q} \mathbf{1}(\mathcal{I}_i \subseteq [b_{d_1}, b_{d_1+1}] \cup [b_{d_2}, b_{d_2+1}])$$ is the same as \[ \sum_{u \in [b_{d_1}, b_{d_1+1}]} N_u +\sum_{u \in [b_{d_2}, b_{d_2+1}]} N_u \geq \sum_{i \in Q} \mathbf{1}(\mathcal{I}_i \subseteq [b_{d_1}, b_{d_1+1}]) +\sum_{i \in Q} \mathbf{1}(\mathcal{I}_i \subseteq [b_{d_2}, b_{d_2+1}]), \] with the latter implied by the existing lower bounds on consecutive unions. For the upper bound constraint for $d_1+1<d_2$,

Since $\sum_{i \in Q} \mathbf{1}(\mathcal{I}_i \subseteq (b_{d_1+1}, b_{d_2})) \leq \sum_{u \in (b_{d_1+1},b_{d_2})} N_u,$ the desired inequality from above would be implied by

or equivalently, $\sum_{u \in [b_{d_1}, b_{d_2+1}]} N_u \leq \sum_{i \in Q} \mathbf{1}(\mathcal{I}_i \cap [b_{d_1}, b_{d_2+1}] \neq \emptyset)$, which is the upper bound inequality for a consecutive union $[b_{d_1}, b_{d_2+1}]$ and which is accounted for already.

\paragraph*{Third step: Optimization} Since it is difficult to optimize subject to integer restrictions on all $N_u$, it makes sense to rewrite the constraints obtained in the second step in terms of shares in the overall sample by dividing both the left-hand and the right-hand sides by $n$. Then for any $k=1,..,K-1$ and any $d=1,...,K-k$, letting $\widehat{\phi}_u := \frac{N_u}{n}$ implies that

We now focus on Gini in the class of linear-fractional inequality indices. Let $\widehat{\psi}_u := \frac{1}{n}\sum_{i \in P} \mathbf{1}(\underline a_i = u)$. These are known shares as they measure the frequency of encountering a point from $\mathcal{U}$ in the point sample ${P}$. Some of them will naturally be zero as not all the points in $\mathcal{U}$ are found in $P$. Then the overall share of every point $u \in \mathcal{U}$ in the sample is the sum of a known $\widehat{\psi}_u$ and unknown $\widehat{\phi}_u$, with the latter being subject to constraints in ((ref)).

One could use various equivalent forms of Gini in terms of $\widehat{\psi}_u$ and $\widehat{\phi}_u$. We write the elements of $\mathcal{U}$ in increasing order as $\mathcal{U}=\{u_1<\cdots<u_{|\mathcal{U}|}\}$. For any $u\in\mathcal{U}$, let $\widehat{\phi}_u$ and $\widehat{\psi}_u$ denote the (unknown) interval-share and (known) point-share at $u$, respectively. Consider then the following form.

This is a quadratic-to-linear-fraction in $\widehat{\phi}_{u_1}, ..., \widehat{\phi}_{|\mathcal{U}|}$ that needs to be optimized in $\widehat{\phi}$ subject to ((ref)). Its optimization is the quadratic-to-linear-fractional program. Let $\mathcal{S}_{\widehat{\phi}}$ denote all $\widehat{\phi}$ that satisfy ((ref)).

In practice the optimization of $G_{2,n}(\widehat{\phi})$ can follow dinkelbach67 via the parametric quadratic programs $$ f_{2; max}(\lambda)= \max_{{\widehat{\phi}} \in \mathcal{S}_{\widehat{\phi}}} \left(\frac{1}{2} \sum_{i =1}^{ |\mathcal{U}|} \sum_{j =1}^{ |\mathcal{U}|} (\widehat{\psi}_{u_i}+\widehat{\phi}_{u_i})(\widehat{\psi}_{u_j}+\widehat{\phi}_{u_j})|u_i-u_j|- \lambda\sum_{i=1}^{|\mathcal{U}|} (\widehat{\psi}_{u_i}+\widehat{\phi}_{u_i})u_i \right)$$ (and $f_{2;\min}(\lambda)$ for the minimum). The unique root $\lambda^*$ where $f_{2;\max}(\lambda^*)=0$\footnote{$f_{2;max}(\lambda)$ is continuous, strictly decreasing and is convex in $\lambda \in [0,1]$.} gives the sharp upper bound; Algorithm (ref) (analogous to Scenario 1, using bisection on $\lambda$ with global optimization at each step) converges to $\lambda^*$. The lower bound is obtained analogously by finding the unique root $\lambda^*_{min}$ where $f_{2;\min}(\lambda^*_{min})=0$, with the same bisection sign rule as in the maximization case: increase $\lambda$ when $f_{2;\min}(\widetilde{\lambda})>0$ and decrease $\lambda$ when $f_{2;\min}(\widetilde{\lambda})<0$.

Like in Scenario 1, in case of linear-fractional measures there may be a way to create computational procedures in terms of the original $\mathbf{y}$ rather than using the process with unknown and known shares. Exact procedures may depend on the exact form of the share. If we take Gini index again as an example, then we can establish the result in Proposition (ref) below that would lead to alternative computational procedures.

Proposition (ref)(a) establishes that the minimizer compresses unknown values $y_{i}$ for $i \in Q$ toward a common threshold $u_0 \in \mathcal{U}$, pushing bounded intervals to their nearer extreme and aligning overlapping ones at $u_0$. The proposition implies an efficient minimization strategy based on a search over candidate thresholds in $ \mathcal{U}$. For each candidate $u_0 \in \mathcal{U}$ construct the assignment for $y_i$, $i \in Q$, according to the pattern in part (a) of the proposition. A single pass over $|\mathcal{U}|\leq 2|Q| + |P|$ candidates, computing the Gini index each time, is guaranteed by the proposition to return the global minimum and a corresponding minimizer $\mathbf{y}^*_{\min}$.

Proposition (ref)(b) shows that a maximizer stretches the distribution by pushing uncertain values away from a common threshold $u_0 \in \mathcal{U}$. Unlike minimization, it does not uniquely prescribe the boundary assignment for the (typically small) subset $Q(u_0) = \{i \in Q : u_0 \in \mathcal{I}_i\}$. When $|Q(u_0)|$ is small, exhaustive enumeration of the $2^{|Q(u_0)|}$ possible boundary assignments for these intervals is feasible and exact. For larger $|Q(u_0)|$, fast steepest-ascent greedy-flip heuristics (or beam search) often recover the global optimum in practice. Additional dominance rules can further prune the search: if $i,i' \in Q(u_0)$ satisfy $\underline{a}_i = \underline{a}_{i'}$ and $\overline{a}_i < \overline{a}_{i'}$, then assigning the narrower interval $i$ to its upper bound $\overline{a}_i$ precludes assigning the wider interval $i'$ to its lower bound $\underline{a}_{i'}$, and symmetrically for the dual case. Crucially, the proposition dramatically simplifies computation even within proportion-based optimization frameworks discussed earlier as all uncertain values lie on the finite boundary set $\mathcal{B}$, eliminating the need for continuous proportions over $\mathcal{U} \setminus \mathcal{B}$, which leads to a potentially large reduction in the number of unknowns.

Asymptotics of sharp bounds

This section establishes asymptotic theory and bootstrap validity for the vector of sharp bound estimators. Our key observation is that each endpoint can be written as the optimal value of a constrained optimization problem over probability measures. Letting the constraint right-hand sides be \(c\), the nuisance parameters indexing the constraint maps (such as group boundaries) be \(\theta\), and letting \(\pi\) and \(q\) denote the proportion and distribution of point observations, we collect these into a vector \(\eta=(c,\theta,\pi,q)\). We show below that the population bound endpoints can be written as $V_\infty(\eta)$ for a function $V_\infty: \mathbb{H} \rightarrow \mathbb{R}^2$. Similarly, the estimators computed in Sections (ref)--(ref) can be written as $\hat V^{(n)}(\hat\eta)$, where $\hat V^{(n)}$ approximates $V_\infty$ and $\hat\eta\to\eta_0$.

After formalizing this representation, the argument proceeds in three steps. First, solution-form results imply that each estimator can be represented by measures supported on a small number of points. Second, this representation yields the expansion \[ \hat V^{(n)}(\hat\eta) - V_\infty(\eta_0) = \big[ V_\infty(\hat\eta) - V_\infty(\eta_0) \big] + o_p(n^{-1/2}). \] Finally, we show that $V_\infty$ is Hadamard directionally differentiable and apply the directional delta method of FangSantos2019 together with a functional CLT for $\sqrt{n}(\hat\eta - \eta_0)$. The presence of only directional differentiability in general means that inference relies on an $m$-out-of-$n$ bootstrap, though the standard nonparametric bootstrap is also valid in special cases.

Throughout, we let circumflexes denote sample analogues and the subscript $0$ population counterparts. Let $\mathcal{Y}=[\underline y,\overline y]\subset\mathbb{R}$ denote the support of $Y$, and let $\mathcal M$ be the set of Borel probability measures on $\mathcal{Y}$, equipped with the Wasserstein-1 ($W_1$) metric. For $y\in\mathcal{Y}$, let $\delta_y$ denote the Dirac measure at $y$, and let $\mathcal J:\mathcal M\to\mathbb{R}$ denote the inequality index.

Writing the Problem Using Measures

Population problem: We first write the population bounds as an optimization problem over probability measures. Let \(\mu_0 \in \mathcal M\) denote the latent distribution of the variable of interest. When both grouped and point observations are present, write \(\mu_0 = (1-\pi_0)\mu_0' + \pi_0 q_0\), where \(q_0\) denotes the distribution of point data and \(\pi_0\) their population share.\footnote{In Scenario 1, degenerate brackets can be absorbed into \(\mu_0'\), so that \(\pi_0 = 0\).} As described in Sections (ref)--(ref), the data imply linear restrictions on \(\mu_0'\). For example, the statement that \(10\%\) of the data lie in a group \(\mathcal G(\theta)\) is encoded as \(\int \mathbf 1\{y\in\mathcal G(\theta)\}\,d\mu(y)=0.1\), while a mean restriction is encoded as \(\int y\,d\mu(y)=50{,}000\). This covers the restrictions used in the paper, including group probabilities, subgroup means, and overall means.

Our estimands are the minimum and maximum of \(\mathcal J\) subject to these restrictions. Because the number of restrictions may increase with \(n\), we allow the population problem to contain countably many of them. Letting \(H,G:\mathcal M\times \ell^\infty(\mathbb N)\to \ell^\infty(\mathbb N)\) collect the equality and inequality restrictions, respectively, we suppose that for each \(j\in\mathbb N\), \(\mu\mapsto H(\mu;\theta)_j\) and \(\mu\mapsto G(\mu;\theta)_j\) are linear-functionals generated by the piecewise-affine integrands described below. This includes moment restrictions such as \(\mu\mapsto\int y\,d\mu(y)\) and indicator restrictions such as \(\mu\mapsto\int \mathbf 1\{y\in\mathcal G(\theta)\}\,d\mu(y)\). We therefore define the feasible set as \[ \mathcal C_\infty(c, \theta) := \Big\{ \mu\in\mathcal M:\ H(\mu;\theta)=c^1,\ \ G(\mu;\theta)\le c^2 \Big\}, \] where \(c=(c^1,c^2)\in \ell^\infty(\mathbb N)\times \ell^\infty(\mathbb N)\) collects the right-hand sides of the equality and inequality restrictions. In Scenario 1A, only equality restrictions are imposed and these are encoded through \(H\). Scenario 1B adds further equality and/or inequality restrictions through \(H\) or \(G\), while Scenario 2 uses inequality restrictions only. In the population, the only information available about \(\mu_0'\) is that \(H(\mu_0';\theta_0)=c_0^1\) and \(G(\mu_0';\theta_0)\le c_0^2\) componentwise. As such, our information about $\mu_0'$ is summarized by the requirement that \(\mu_0' \in \mathcal C_\infty(c_0,\theta_0)\).

This representation lets all sampling uncertainty enter through \(\eta\). Define \(F(\mu;\pi,q):=\mathcal J((1-\pi)\mu+\pi q)\). Let \(\mathcal F_0:=\{x\mapsto \mathbf 1(x\le t): t\in\mathbb R\}\), and when convenient identify \(q\) with its cumulative distribution function viewed as an element of \(\ell^\infty(\mathcal F_0)\). Write \[ \mathbb H:=\big(\ell^\infty(\mathbb N)\times \ell^\infty(\mathbb N)\big)\times \ell^\infty(\mathbb N)\times \mathbb R \times \ell^\infty(\mathcal F_0), \] equip \(\mathbb H\) with the product sup norm, and for \(c=(c^1,c^2)\) write \(\|c\|_\infty:=\max\{\|c^1\|_\infty,\|c^2\|_\infty\}\). When \(\theta\) has finitely many components, we view it as an element of \(\ell^\infty(\mathbb N)\) by padding with zeros. Then letting \[ V_{\infty}^{\sup}(\eta):=\sup_{\mu\in\mathcal{C}_\infty(c,\theta)} F(\mu;\pi,q), \qquad V_{\infty}^{\inf}(\eta):=\inf_{\mu\in\mathcal{C}_\infty(c,\theta)} F(\mu;\pi,q), \] we define our estimand by \(V_\infty(\eta_0)\), where \(V_\infty(\eta) := \big( V_{\infty}^{\inf}(\eta),\, V_{\infty}^{\sup}(\eta) \big)^{\top} \in \mathbb{R}^2\). Thus \(V_\infty(\eta_0)\) is the vector of lower and upper bounds on the inequality index implied by the population restrictions and the point-data component.

Sample problem: The same representation applies at the sample level. Conditional on the observed coarsening pattern, the sample problem gives sharp bounds for the inequality index of the realized sample, and under sampling it serves as a plug-in estimator of the population bounds. The original formulations in Sections (ref)--(ref) optimize over vectors \(y\in\mathcal{Y}^n\), but because both \(\mathcal J\) and the restrictions depend on \(y\) only through its empirical measure, the same problem can be written as optimization over empirical measures \(\mu_y:=n^{-1}\sum_{i=1}^n\delta_{y_i}\).\footnote{In the measure formulation there is no need to impose an ordering on \(y\).}

In addition, \(\eta_0\) is unobserved. The data yield an estimator \(\hat\eta:=(\hat c,\hat\theta,\hat\pi,\hat q)\), where \(\hat c\) collects the estimated right-hand sides of the restrictions, \(\hat\theta\) collects any additional estimated nuisance quantities entering \(H(\cdot;\theta)\) and \(G(\cdot;\theta)\) (for example empirical cutoffs, quantiles, or smooth transforms thereof), and \(\hat\pi\) and \(\hat q\) are the empirical share and distribution of point data. The distinction between \(c\) and \(\theta\) is therefore only that \(c\) enters on the right-hand side, whereas \(\theta\) enters through the restriction maps. Let \(\mathcal C_J(c,\theta)\) denote the feasible set in the \(J\)-th sample problem, obtained by imposing all equality restrictions together with the inequality restrictions included at stage \(J\); write \(J_n\) for the choice used at sample size \(n\). For integers \(a,b\), let \(\mathcal C_J^{(a)}(c,\theta)\subseteq \mathcal C_J(c,\theta)\) be the subset supported on at most \(a\) points, and let \(\mathcal C_J^{(a,b)}(c,\theta)\subseteq \mathcal C_J^{(a)}(c,\theta)\) further require each mass to be an integer multiple of \(1/b\). The exact indexing convention for the \(J\)-th problem is given in the Appendix.

As a result, the sample restrictions on \(y\) are equivalent to \(\mu_y \in \mathcal{C}_{J_n}^{(n,n)}(\hat c, \hat \theta)\). We therefore define \[ \hat V^{(n)}(\eta):=

, \] so that the computed estimator is \(\hat V^{(n)}(\hat\eta)+\varepsilon_n\), where \(\varepsilon_n\in\mathbb R^2\) denotes numerical optimization error and satisfies \(\|\varepsilon_n\|_2=o_p(n^{-1/2})\). Our theory studies the convergence of \(\hat V^{(n)}(\hat\eta)\) to \(V_\infty(\eta_0)\) as \(n\to\infty\), allowing \(J_n\to\infty\) as additional inequality restrictions are included.

Regularity conditions

Each regularity condition below controls one part of the difference between \(\hat V^{(n)}(\hat\eta)\) and \(V_\infty(\eta_0)\). For integers \(a,b,J\), let \(V_J^{a,b}(\eta)\) denote the vector of lower and upper endpoints obtained by replacing \(\mathcal C_J(c,\theta)\) with \(\mathcal C_J^{(a,b)}(c,\theta)\). When \(\pi_0>0\) we take \(\mathbb D_0:=\mathbb H\). In the case \(\pi_0=0\), fix \(q_0:=q^\dagger\) for some \(q^\dagger\in\ell^\infty(\mathcal F_0)\), and define \(\hat q:=q^\dagger\) on \(\{\hat\pi=0\}\), since \(F(\mu;0,q)\) does not depend on \(q\), and let \[ \mathbb D_0:=\big(\ell^\infty(\mathbb N)\times \ell^\infty(\mathbb N)\big)\times \ell^\infty(\mathbb N)\times \mathbb R\times \{0\}\subset\mathbb H. \] Assumption (ref) formalizes the idea that, for purposes of $\sqrt n$ asymptotics, the optimization problem defining each bound can be effectively reduced to measures with a small number of support points.

Assumption (ref) requires that restricting attention to measures supported on at most \(k_n\) points changes the vector of endpoints by at most \(o(n^{-1/2})\), uniformly over local perturbations of \(\eta\). In the settings considered here, this follows directly from the characterization results in earlier sections. Recalling that \(D\) is the number of groups and \(q_1+q_2\) is the number of constraints, one may take \(k_n\le 2D\) in Scenario 1A, \(k_n\le 2D+q_1+q_2\) for linear-fractional indices in Scenario 1B, and \(k_n\le |\mathcal B|+1\) in Scenario 2; see Propositions (ref), (ref), and Theorems (ref)--(ref). Aside from numerical optimization error, the approximation error is then exactly zero. The substantive content of Assumption (ref) is that \(k_n\) is tied to the complexity of the restriction set rather than to the ambient dimension \(N\). This holds in all of the above settings, since \(k_n\) is bounded by features of the restriction set.

Assumption (ref) places restrictions on the constraints allowed in the problem. Write \(\theta=(\tau,\gamma)\), where \(\tau\) collects the cutoff values at which a restriction may change form and \(\gamma\) the remaining nuisance quantities. After relabeling coordinates if needed, each cutoff is either a fixed endpoint or a coordinate of \(\tau\). For each \(n\), let \(t_{0,n}(\tau)<\cdots<t_{s_n,n}(\tau)\) denote the relevant cutoffs, and write \[ I_{1,n}(\theta):=[t_{0,n}(\tau),t_{1,n}(\tau)], \qquad I_{d,n}(\theta):=(t_{d-1,n}(\tau),t_{d,n}(\tau)],\quad d=2,\ldots,s_n. \] Let \(\mathcal L_{\mathrm{eq},n}\) and \(\mathcal L_{\mathrm{ineq},n}\) index the equality and inequality restrictions, let \(f_{u,n}(\cdot;\theta)\) denote the integrand for restriction \(u\), and let \(\mathcal L_{\mathrm{step,eq},n}\subseteq\mathcal L_{\mathrm{eq},n}\) collect those equality restrictions whose integrands are constant on each partition interval. This covers grouped-share and overlap restrictions, mean restrictions, and ratio or Lorenz restrictions; when groups are defined by quantiles of the latent distribution, the estimated quantiles enter through \(\tau\). In Scenario 2, \(\mathcal L_{\mathrm{step,eq},n}=\varnothing\).

Assumption (ref) says that, after partitioning \(\mathcal{Y}\) at the relevant cutoffs, each restriction depends on \(\mu\) within an interval only through the probability assigned to that interval and its first moment. This is the feature used in the proofs to replace \(\mu\) within an interval by a discrete approximation while preserving the restrictions. We take \(\tau\) to be the cutoffs themselves, so the local interval bounds used later are affine in \(\tau\). Part (iii) ensures that the probabilities of the subintervals determined by the relevant cutoffs can be recovered from cumulative-probability equalities, possibly after adding equalities implied by the original system; in Scenario 1 these are implied by grouped shares, while in Scenario 2 the condition is vacuous. As such, these constraints do not need to be separately added into the sample or population problem, beyond the constraints we already consider in each of our scenarios. The assumption therefore covers the grouped-share, overlap, mean, ratio, and Lorenz restrictions used in Sections (ref)--(ref).

Assumption (ref) is a regularity condition on the objective \(F\). Part (i) is global. For part (ii), fix \(n\), one of the two endpoints, and an optimizer \(\mu^\star\) of the corresponding \(k_n\)-point problem at \(\eta_0\). On a compact neighborhood of \(\mu^\star\), write nearby feasible \(k\)-point measures as \(\mu_x=\sum_{j=1}^k p_j\delta_{z_j}\), with the support points remaining in their current partition cells, and define \(\Phi_k(x;\pi,q):=F(\mu_x;\pi,q)\).

Assumption (ref)(i) is a \(W_1\)-Lipschitz condition for the objective. Part (ii) requires local smoothness of the finite-dimensional \(k_n\)-point problems around each relevant optimizer, together with a first-order expansion in \((\pi,q)\) that is uniform on compact sets of directions. The Appendix verifies these conditions for the objective classes used in the paper: smooth functions of finitely many moments (including the mean log deviation, the Theil index, \(\mathrm{GE}_\alpha\), the Atkinson class, the Kolm class, and linear moment functionals, with positivity restrictions where needed), quantile ratios such as the \(90/50\) ratio under local regularity at the relevant quantiles, and the Gini and Hoover indices under local regularity of \(F_q\) and a mean bounded away from zero.

Assumption (ref) is our key condition on the sampling process. We require that \(\hat\eta\) satisfy a functional central limit theorem as an element of \(\mathbb H\), and that the equality restrictions in \(\mathcal L_{\mathrm{step,eq},n}\) have right-hand sides on the empirical grid \(n^{-1}\mathbb Z\). We impose the same conditions on the bootstrap analogue \(\eta_m^*:=(c_m^*,\theta_m^*,\pi_m^*,q_m^*)\), computed from the \(m\)-out-of-\(n\) resample described below.

A convenient sufficient condition is that the data are coarsened versions of underlying i.i.d.\ random variables \(W_1,\dots,W_n\). Let \(P_0\) denote the distribution of \(W_i\). In Scenario 1, \(W_i:=Y_i\). In Scenario 2, \(W_i:=(D_i,Y_i,L_i,U_i)\), where \(D_i\) indicates whether observation \(i\) is a point or interval observation and \((L_i,U_i)\) are the lower and upper bounds on \(Y_i\). Let \(\mathcal W\) be a class of measurable functions of \(W_i\), define \(\nu_0:=(\int g\,dP_0)_{g\in\mathcal W}\) and \(\hat\nu:=(n^{-1}\sum_{i=1}^n g(W_i))_{g\in\mathcal W}\), and suppose \(\eta_0=\Psi(\nu_0)\) and \(\hat\eta=\Psi(\hat\nu)\) for a Hadamard differentiable map \(\Psi:\ell^\infty(\mathcal W)\to\mathbb H\). If \(\mathcal W\) is \(P_0\)-Donsker, then \(\sqrt n(\hat\eta-\eta_0)\Rightarrow Z\) for a tight, mean-zero Gaussian \(Z\); whenever the relevant statistics can be recomputed from resampled sampling units, the same conditions yield the bootstrap limits for \(\eta_m^*\). If \(\pi_0=0\), fixing \(q\) as above implies that the \(q\)-coordinate has no first-order effect, so \(Z\in\mathbb D_0\) almost surely. This setup covers empirical shares, empirical moments, empirical distribution functions, ratios of empirical averages, and empirical quantiles under the usual local positive-density condition at the relevant quantiles.

In the designs studied here, these empirical-process conditions are routine. In Scenario 1, the relevant coordinates of \(\eta\) are built from threshold indicators and bounded moment functions such as \(\mathbf 1\{Y_i\le t\}\) and \(\mathbf 1\{Y_i\in I\}Y_i\). In Scenario 2, they are built from point-data terms such as \(D_i\) and \(D_i\mathbf 1\{Y_i\le t\}\), and interval-data terms such as \(\mathbf 1\{D_i=0,L_i\ge a,U_i\le b\}\) and \(\mathbf 1\{D_i=0,L_i\le b,U_i\ge a\}\). Indicator classes indexed by thresholds or rectangles, and bounded products of such indicators, are standard Donsker classes. The grid conditions are also natural: in Scenario 1, \(\hat c_u\) and \(c_{m,u}^*\) are empirical group shares and hence exact multiples of \(1/n\) and \(1/m\), while in Scenario 2, \(\mathcal L_{\mathrm{step,eq},n}=\varnothing\).

Assumption (ref) controls the error from replacing the full inequality system by the \(J_n\)-constraint problem. The equality restrictions are eventually fixed; what may grow with \(n\) is the set of included inequalities. The assumption says that, uniformly over \((c,\theta)\) near \((c_0,\theta_0)\), any measure feasible for the \(J_n\)-problem violates the omitted inequalities by at most \(o(n^{-1/2})\).

Assumption (ref) says that omitted inequalities are relaxed by at most \(\kappa_{J_n}=o(n^{-1/2})\) uniformly. Because the inequalities included in the \(J_n\)-problem are imposed exactly, the display concerns only omitted inequalities. If the full inequality system is finite, then all inequalities are eventually included and \(\kappa_{J_n}=0\) from some \(n\) onward; this already covers both applications in Section (ref). More generally, the condition holds whenever the included inequalities uniformly approximate the omitted ones.

Under Assumption (ref), each endpoint problem reduces to a finite-dimensional optimization problem over the masses and support points of a \(k_n\)-point distribution, with \(k_n\) possibly increasing with \(n\). We now impose a regularity condition that delivers a constraint qualification for these reduced \(k_n\)-point problems; the formal version is given in Appendix (ref) as Assumption (ref).

Assumption (ref) has two roles. The first part of the assumption is a uniform Slater condition, which provides a uniformly feasible benchmark measure - both for the full inequality system and on a fixed support. We use this in the approximation arguments and to obtain a strict feasible direction for the local finite-dimensional problem. The rank part of Assumption (ref) is used, after fixing each support point's partition cell, to verify a standard constraint qualification (Robinson's condition), to bound the associated Lagrange multipliers uniformly, and to make those multipliers unique. In Scenario 2 there are no non-step equality restrictions and no step-equality restrictions beyond the simplex constraint. In our applications grouped-share restrictions pin down block masses, mean and Lorenz-type restrictions are controlled by moving a small fixed number of interior atoms, and overlap restrictions enter only through the locally binding inequalities.

Taken together, Assumptions (ref)--(ref) hold for the constraint types and data structures in our applications, and for the objective classes verified in the Appendix.

Main Results

Because the assumptions are local, take neighborhoods around \(\eta_0\) small enough that Assumptions (ref)--(ref) hold simultaneously, with Assumption (ref) applied on the corresponding \((c,\theta)\)-neighborhood. Let \(\mathbb D_0\subseteq\mathbb H\) denote the tangent set associated with \(Z\). Our main result is that \(V_\infty\) is Hadamard directionally differentiable.

Proposition (ref) follows by proving that an envelope theorem applies to the finite-dimensional \(k_n\)-point problem, together with approximation results that transfer differentiability to \(V_\infty\). Proposition (ref) then shows that the errors arising from the finite-support reduction, the discretization of the masses, truncation of the constraint system, and numerical optimization are all \(o_p(n^{-1/2})\) in \(\mathbb R^2\). Consequently, \[ \hat V^{(n)}(\hat\eta)-V_\infty(\eta_0) = \bigl[V_\infty(\hat\eta)-V_\infty(\eta_0)\bigr] + o_p(n^{-1/2}), \] so the directional delta method yields the following result.

The limiting distribution is Gaussian when \(V'_{\eta_0}\) is linear, but need not be otherwise. In the latter case the standard nonparametric bootstrap may fail, so we use the \(m\)-out-of-\(n\) bootstrap of Shao1994. Let \(\eta_m^*\) be the analogue of \(\hat\eta\) recomputed on an \(m\)-out-of-\(n\) resample drawn with replacement from the original sampling units. Our theory requires only that \(m\to\infty\) and \(m/n\to0\); if \(V'_{\eta_0}\) is linear, the standard nonparametric bootstrap (\(m=n\)) is also valid.

On each bootstrap sample we solve the same optimization problem as in the original sample, writing \(\hat V_m^*:=V_{J_n}^{m,m}(\eta_m^*)\). Since \(c_{m,u}^*\in m^{-1}\mathbb Z\) conditionally almost surely for every equality restriction \(u\in\mathcal L_{\mathrm{step,eq},n}\) whose integrand is constant on each partition interval, the bootstrap problem is well defined. As in the original sample, any numerical optimization error can be absorbed into an \(o_{P^*}(m^{-1/2})\) remainder and is suppressed in notation.

Implementation is straightforward in Scenarios 1A and 2: one resamples the sampling units and recomputes \((\hat c,\hat\theta,\hat\pi,\hat q)\). This works because our observed data in those scenarios can be viewed as empirical averages over the sampling units. In some Scenario 1B applications, however, the auxiliary aggregates cannot be reconstructed from a nonparametric resample without additional information such as microdata, replicate estimates, or a model for the covariance structure of the reported statistics. A practical alternative is to bootstrap the weaker problem that omits those auxiliary constraints and to report the sharper point estimate separately; this yields inference for an outer identified interval rather than for the sharper endpoints. As such, critical values computed from \(\sqrt{m}\{\hat V_m^*-\hat V^{(n)}(\hat\eta)\}\) can then be used to construct asymptotically valid simultaneous inference for the lower and upper endpoints in each of these scenarios.

Example Applications

This section illustrates the framework with two applications, one for each observational scenario. The first uses the English Longitudinal Study of Ageing (ELSA), where many households report point values while others provide respondent-specific intervals generated by unfolding brackets (Scenario 2). The second uses published U.S.\ income distribution tables from the mid-twentieth century, which report frequencies in non-overlapping income brackets and, in some years, additional aggregates such as subgroup means or selected quantiles (Scenario 1). In both applications we report conventional point estimates based on common missing-data treatments alongside our sharp identified intervals. This comparison highlights when imputation yields deceptively precise conclusions and when auxiliary linear information materially tightens identification.

Household level wealth data

Household wealth surveys routinely face item nonresponse: respondents may be unwilling to report exact amounts or may not know them precisely. To reduce missingness while limiting respondent burden, many wealth surveys use unfolding brackets: if a respondent does not provide a point value, they are routed into a short sequence of Yes/No threshold questions (e.g.\ “Is it more or less than \(X\)?”) that places the value in an interval. The thresholds are typically randomized to mitigate anchoring and response-order effects. JusterSmith1997 and JusterSmithStafford1999 discuss the design and performance of this approach, which is now used in a range of surveys including HRS, PSID, ELSA and the Survey of Health, Ageing and Retirement in Europe.

This design naturally produces both point-and-interval data with respondent-specific intervals (Scenario 2). We use the 2018/19 wave of ELSA (see ELSA_IJE2013), focusing on households whose financial respondent is aged 50--74. This yields \(n=4{,}422\) observations. We study two measures of liquid savings: (i) a narrow single-question measure with relatively straightforward interval data structure and (ii) a broader composite measure that aggregates three components and therefore exhibits more interval complexity. For each measure we compare sharp bounds to commonly used imputation-based point estimates.

The most straightforward, and narrowest definition of liquid precautionary savings is balances in savings and checking accounts at banks or building societies (a mutual financial institution in the UK that serves a similar role as community banks or credit unions in the US). ELSA financial respondents are asked to give the value of their household's total current balances in bank and building society accounts in a single question response, returning a value of 0 if the household does not have such accounts. Of 4,422 observations in our sample, respondents were able and willing to give point values for this variable in 3,827 (86%) of the cases and sample statistics for this point value sample are given in panel A of Table (ref). Of the remaining 595 observations with individual-specific intervals, 233 had some kind of bounded interval data generated from the unfolding bracket procedure, and a further 362 had interval data generated from the unfolding brackets that was unbounded at the top. For simplicity here we have taken twice the value of the top band as the upper interval limit for all unbounded cases. Specifically, for this variable in ELSA, the final bracket point ends at \textsterling 150,000, so we have set the value for an open band to be \textsterling 300,000.\footnote{Alternative choices are straightforward to implement.} Overall, the 595 observations generate only 11 distinct intervals.\footnote{These unique intervals are $[0,999]$, $[0,4999]$, $[0,19999]$, $[0,300000]$, $[1001,4999]$, $[1001,300000]$, $[5001,19999]$, $[5001,300000]$, $[20001,149999]$, $[20001,300000]$, $[150001,300000]$. }

In reality, other types of financial assets are functionally equivalent to balances in bank and building society accounts, in terms of offering similar liquidity, lack of risk and comparable interest rates. In the context of the UK savings landscape, so-called National Savings products and tax-advantaged Individual Savings Accounts (that are held in the form of cash as opposed to stocks and shares) fall into this category. Thus a more comprehensive measure of liquid savings should include any such balances held by households. In cases such as these researchers need to aggregate over multiple (in this case three) interview questions, each with differing patterns of missingness, in order to arrive at a total. Analogously to the narrow definition above, a respondent may still have a point value for this variable but this now happens if and only if all three components have point values. If at least one of the variables in the set has interval values for a respondent, then this broader definition of liquid savings will take an interval value with its lower value being the sum of the lower bracket points across the three variables and the same being true for the upper value. If a component variable is reported with a point value then it is treated as both the lower and the upper bracket in this calculation.

The consequence of this aggregation of variables, each with their own pattern of missingness, leads to more instances of interval data and, importantly, considerably more complexity in the nature of the intervals. Out of the 4,422 observations in our sample, 3,708 have point values for the broader measure of liquid savings and 714 individuals have interval data. Panel B of Table (ref) gives summary statistics for the subset of individuals with point data. As for the observations with the interval data for this variable, there are 199 unique interval types; the five most frequent are $[0, 300000]$ (157 times), $[0, 340000]$ (109 times), $[0, 999]$ (39 times), $[1001, 4999]$ (34 times), $[5001, 19999]$ (34 times).

In order to put the bounds we will compute in a relevant context we begin by computing values for the Gini coefficient in this household savings data using commonly-used approaches to deal with missing data. The most simplistic of these is simply to drop the cases where point data are missing and calculate the Gini on the basis solely of the continuous part of the sample. Researchers concerned with whether the data are missing at random, however, typically want to include information from the whole sample and therefore use a variety of imputation based methods.

The most basic imputation approaches just assign either the midpoint of the relevant band, or else the mean value computed from the continuous sample lying within the band, to all observations where the continuous data is missing. These have the advantage of using the non-missing data but are somewhat inappropriate for the study of inequality since they, by definition, reduce the variance in the overall sample. Hot-deck imputation approaches, where a random donor observation from the pool of continuous observations lying within the relevant band is assigned as the imputed value for each missing observation, do not suffer from this problem. For the purposes of this exercise, we consider a single hot-deck imputation measure and a multiple hot-deck imputation measure based on ten imputations that reduces the dependency on any single draw (see Rubin1987 for an overview).

The results from these exercises for both the narrow and broad definitions of liquid savings are presented Table (ref). The level of inequality is high compared to that typically seen in income distributions but this is entirely to be expected, and particularly so for households at the older end of the cycle where differences in income and expenditure will have accumulated up over time. Looking across the various measures, mean and midpoint imputation methods yield a lower measure of inequality as expected, while the inequality in the continuous (non-missing) sample is higher than that measured by the hot-deck imputation methods that take into account the unfolding bracket information for the missing cases. The Gini coefficients (and their s.e.) for the broader savings definition presented in Panel B of the table are lower than for the narrow single-variable case, which is consistent with some substitutability between the three different types of assets in the class.

For each of the hot-deck imputation methods we ran the procedure 1000 times, reporting the mean values and bootstrap confidence intervals for the Gini coefficient in the table. For every hot-deck run of 1000 repetitions the distribution is more compressed in the multiple imputation case, as expected. However, the aggregated values across 1000 repetitions have similar bootstrap standard errors.

The sharp bounds on the Gini interval computed using our methodology for Scenario 2 are given in Table (ref).\footnote{As expected, the sharp upper bound on the Gini is sensitive to our treatment of the unbounded top interval. When we set the upper interval limit for entries with the open right bracket to 1,500,000 (instead of 300,000), the sharp upper bound on the Gini for the narrow definition of liquid savings changes to 0.8829.} The bounds are considerably wider than either the confidence intervals on Gini coefficients calculated on the imputed data or on the range of values for the Gini covered by the 1000 draws of the hot-deck imputation samples. The bootstrap standard error for the length of the sharp Gini interval in the narrow case is 0.0039 and the confidence interval for that length constructed on the basis of the bootstrap 2.5th and 97.5th percentiles is $[0.0699,0.0853]$. The bootstrap distributions of lower and upper sharp bounds for the narrow savings definition, as well as the sharp Gini interval length, are given in Figure (ref).

Note that while the imputation based methods in Table (ref) indicate a lower Gini and smaller standard error in the broad savings case than for narrow savings measure, comparing Panel A and Panel B of Table (ref) shows that the width of the bounds on the Gini for the broad definition are actually wider, reflecting the increasing missingness and more complex interval nature of the broader data --- the width of the bounds is 0.0895 in comparison to 0.0774 in the narrow case.

From the computational perspective we maximize the Gini index represented as $G_{2,n}(\widehat{\phi})$ in ((ref)) subject to the exhaustive list of inequality constraints in ((ref)) and one equality constraint capturing the total number of interval records. In the broad savings definition case $|\mathcal{U}|=961$ (in contrast, in the narrow definition $|\mathcal{U}|=521$), so our optimization is over 961 variables and, since $|\mathcal{B}|=240$ (in the narrow definition $|\mathcal{B}|=10$), we have 57,368 inequality constraints, in contrast to 88 constraints in the previous narrow definition case.

This shows how in our generic computational approach of Section (ref) the problem naturally gets more complicated as we incorporate more variables with potentially missing data. This happens through both an increase in constraints (driven by the number of unique intervals for individuals in set $Q$ for aggregated variable) and an increase in the number of unknowns (driven by $|\mathcal{U}|$) even though for the latter one could incorporate our results of Theorem (ref) to reduce that dimensionality. In practice, a measure of total net financial wealth might bring together ten to fifteen subcomponents, so the missing-data structure can become quite intricate. Even so, our algorithm works in exactly the same way. This mainly means it raises computation time but does not alter the core logic of the procedure.

Income inequality in the U.S. from historical data

In our second example we consider a different type of interval data corresponding to the more straightforward scenarios 1A and 1B in our previous theoretical analysis. This is the scenario where data are given in non-overlapping intervals, as is often the case when information on how a distribution breaks down into summary intervals is presented in distributional tables typically presenting the number of observations falling into different ranges, perhaps with additional subgroup means or medians. For our empirical example we consider the case of historical data on US income inequality in the early to mid twentieth century and use data from the distributional tables covering various years from 1929--1971 that were produced by government agencies prior to their public release of microdata for analysis by researchers.

Early descriptions of the income distribution in America were developed by Selma Goldsmith from the US Department of Commerce Office of Business Economics, who produced various summary distributional tables for select years between 1929 and 1950 (see Goldsmith1954 and Goldsmith1958). Subsequently, the Office of Business Economics (OBE) produced distributional tables as part of their regular outputs in many years from 1958 onwards (OBE1958USIncomeOutput). Inequality statistics were not produced, and in addition the distributional tables were in different formats for different years. The 1929 income distribution was presented in 8 income categories, with additional aggregate income ratios presented for the bottom 40%, 40-60%, 60-80%, 80-95% and top 5%. Data covering select years 1935 to 1950 (in Goldsmith1958) used the same 8 income categories but instead presented subgroup mean incomes for the bottom 40%, 40-60%, 60-80%, 80-95% and top 5% groups. In contrast to this, data from the OBE (OBE1958USIncomeOutput) presented the distribution in 13 categories along with subgroup means until 1954, then with the addition of the median, 20th, 40th, 60th, 80th and 95th quantiles for the years 1955--1962, before switching to 25 categories along with subgroup means within each category from 1964 onwards.

These data provided a wealth of information for scholars of historical inequality in the US (see, for example, Budd1970 or Lindert2004) but did not provide statistics, such as the Gini, which could be combined with calculations from modern microdata to produce long run series. While some individual studies have used assumptions on the underlying distribution to allow an estimate of the Gini from subsets of the available distributional statistics, our method allows for the computation of sharp bounds on Gini coefficients or quantile ratios for the US income distribution incorporating all the information from the published statistics, regardless of the fact that the nature of the information available changes over time.

Sharp Gini bounds computed from these historical US Income data are given in Figure (ref). The figure presents our calculations for two time-series of upper and lower bounds, depending on whether the source data come from Goldsmith1958 covering years from 1929 to 1950 or the various OBE sources covering for the period from 1944 onwards. In addition, and as a comparison, we present the point estimates for the Gini over this time period that have been estimated by Lindert2004.

The dramatic fall in the Gini coefficient from the late 1920s up to the beginning of the second world war is immediately apparent, despite our bounds being relatively wide in these early years due to the lack of granularity of the distributional tables that were produced. The lower bound on the Gini for 1929 is also higher than current estimates, despite other measures of inequality such as the top 1% share of income having already returned to, or surpassed, their 1920s Gilded Age levels (see, for example, Saez2018). The figure also shows how the width of our bounds reduces as additional information is added to the distributional tables. Particularly noticeable is that the addition of subgroup means as well as the frequencies in different income categories, from 1964 onwards, tightens the bounds significantly.

As discussed previously, our procedures can also deliver bounds for quantile ratios of the distribution. In the case when the frequencies in various intervals are all that is known then such bounds would be straightforward to compute but when additional information on, for example, subgroup means or income shares needs to be taken into account then it becomes a non-trivial operation. We apply our algorithm to produce sharp bounds for the 90/50 quantile ratio from these same historical US Income data, presented in Figure (ref) in order to examine the role of upper middle inequality in driving these changes in the Gini.\footnote{It is of course possible to compute bounds on any quantile ratios. We choose to present the 90/50 since the upper bound for the other commonly used 90/10 and 50/10 ratios is extremely large or even infinite in the earliest years of the data where published distributional information shows that it is possible that considerable numbers of families at the bottom of the distribution had extremely low income or even zero income.} Unlike the bounds for the Gini, the bounds on the 90-50 ratio do not narrow so much with time when the published distributional tables begin to be produced using more categories.

Conclusion

We develop a framework for sharp partial identification and inference on inequality indices when outcomes are observed only coarsely, through grouped data or individual-level intervals, possibly combined with auxiliary linear restrictions. For many commonly used indices, a linear-fractional representation after sorting reduces bound computation to tractable linear (or closely related) programs. Our solution-form characterizations also support the asymptotic analysis in Section (ref): under the conditions stated there, the bound map is directionally differentiable and the estimated endpoints admit $\sqrt{n}$ inference, with $m$-out-of-$n$ bootstrap validity in general, and Efron bootstrap validity if the optimizer is unique.

The applications illustrate two practical messages. First, conventional point estimators based on midpoints or imputation can understate the true identification uncertainty when interval information is substantial. Second, incorporating credible auxiliary information (e.g.\ subgroup means or Lorenz-curve restrictions) can materially tighten identification, sometimes turning wide bounds into informative intervals. The observational regimes studied here, such as grouped tables and individual-specific interval responses, are common in applied work, and additional linear restrictions can be incorporated within the same framework.

Extending these computational ideas to additional Schur-convex inequality indices is one promising direction for future work. More broadly, the linear constraint framework developed here is not specific to scalar inequality indices: the same feasibility characterization, finite-dimensional reduction, and inferential approach may be applicable to a wider class of distributional functionals (some of which could be linear combinations of inequality functionals) in settings where outcomes are only partially observed, and we view this as a fruitful avenue for future research.

\setcounter{page}{1}