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.
107,506 characters · 20 sections · 43 citation commands
IV regression with distribution-valued outcomes
\onehalfspacing
This paper introduces an instrumental variables (IV) framework for settings where the outcome of interest is a distribution. A leading example is when researchers are interested in the effect of a treatment that varies at the group level on the distribution of an outcome within groups. Analyzing effects on the entire distribution of an outcome within groups, rather than on simple group-level averages, provides a better understanding of the impact of a treatment and allows for studying important distributional issues, such as inequality and tail risk. For instance, building on autor2013china, chetverikov2016iv (CLP\xspace henceforth) estimated the effect of import competition on the distribution of local wages and found that low-wage earners were most affected by increased import competition. In this application, the treatment (import competition) varies at the group (commuting zone (CZ)) level, and the outcome is the distribution within groups (the wage distribution within a CZ). Other examples include policies implemented at the firm, hospital, school, county, state, or country level that affect the entire distribution of employee, patient, student, or population outcomes (see also van2025regression). In many applications, these group-level treatments are endogenous, which motivates the use of IVs for estimating causal effects.
We propose a flexible IV regression framework allowing researchers to estimate the effects of real-valued treatments on distribution-valued outcomes, using real-valued instruments. We refer to the proposed procedure as IV Fr\'echet regression (\ensuremath{IVFR}\xspace), because it can be viewed as an extension of the Fr\'echet framework for regression on metric spaces of petersen2019frechet to the IV setting. We will use the grouped data terminology throughout the paper, but emphasize \ensuremath{IVFR}\xspace is very general and can be applied to panel data as well as any other setting where the outcomes of interest are distribution-valued.
We consider a linear structural model for the (known or estimated) group-level quantile function corresponding to the distribution-valued outcome of interest. We show that this structural quantile function is identified as the solution to a Fr\'echet IV regression problem in 2-Wasserstein space. This identification result suggests a computationally straightforward plug-in estimation strategy: (i) construct plug-in estimates of the IV weights, (ii) compute IV-weighted average quantile curves at each covariate value, (iii) project each curve onto the space of valid quantile functions---solving the sample Fr\'echet IV problem---and (iv) recover coefficient functions by OLS. We show that the resulting estimator converges weakly to a zero-mean Gaussian process and propose uniform inference procedures based on the multiplier bootstrap.
A key feature of \ensuremath{IVFR}\xspace is the projection step (iii), which has several desirable properties. First, it ensures that each estimated conditional distribution is a valid probability measure and that the estimator has the interpretation of an IV-weighted Wasserstein barycenter. Second, it provably decreases the error in estimating the structural quantile function in finite samples, which in turn leads to a finite-sample improvement in coefficient estimation. Finally, under a weak monotonicity condition on the population IV-weighted quantile functions, the projection step does not affect the asymptotic distribution, so that bootstrap inference remains valid.
We evaluate the finite-sample performance of \ensuremath{IVFR}\xspace in Monte Carlo simulations. The projection reduces IMSE by up to 63% relative to existing methods under weak to moderate instruments, with the gain diminishing to below 1% under stronger instruments. Our pointwise and uniform inference confidence bands exhibit good finite sample coverage.
We demonstrate the usefulness of \ensuremath{IVFR}\xspace in two empirical applications. First, we revisit the distributional wage effects of Chinese import competition in CLP\xspace. We find that \ensuremath{IVFR}\xspace produces 9--10% narrower pointwise confidence bands than CLP. Using our new uniform confidence bands, however, we find no evidence that wages are reduced at the very bottom of the distribution, but only between the 10th and 35th quantile. Second, we revisit the estimation of the causal effect of food stamps on the county birth weight distribution in melly2025minimum and find no significant effects across the distribution.
\paragraph{Literature.} We contribute to several strands of the literature. Our first contribution is to the literature on distributional effects in panel and grouped data settings koenker2004quantile,canay2011simple,galvao2015efficient, galvao2016smoothed,chetverikov2016iv,galvao2020unbiased, chen2023group, gunsilius2023distributional, pons2024quantile, torous2024optimal, chen2025quantile,melly2025minimum. We show that when the object of interest is the effect of a group-level variable on a group-level distribution, the problem can be formulated naturally as a regression with distribution-valued outcomes. This perspective places existing grouped quantile IV estimators in a Fr\'echet regression framework: in the absence of individual-level covariates, the estimators of CLP\xspace and MP\xspace hausman1981panel coincide with the unprojected version of our estimator. The Fr\'echet formulation shows that the estimand is an IV-weighted Wasserstein barycenter, yields fitted distributions that are valid by construction, and leads to a natural projection step that improves the finite-sample performance while leaving the first-order asymptotic distribution unchanged. It also allows us to derive novel theoretical results: we establish the properties of \ensuremath{IVFR}\xspace under misspecification of the linear quantile model and propose a multiplier bootstrap for constructing uniform confidence bands. See Section (ref) for a formal discussion of the relationship of our model and CLP\xspace and MP\xspace and Section (ref) for a simulation comparison.
Our second contribution is to the literature on Fr\'echet regression petersen2019frechet and to the recent literature extending Fr\'echet regression and related methods to causal inference settings lin2023causal, katta2024interpretable, kurisu2024geodesic, hoshino2024functional, bhattacharjee2025doubly, van2025regression,zhou2025geodesic, kurisu2025regression, kurisu2026lee. We contribute to this literature by providing a linear IV framework for estimating treatment effects of endogenous group-level treatments; as well as by establishing asymptotic normality and deriving uniform confidence bands for this novel estimator, complementing related results for global and local Wasserstein-Fr\'echet regression under exogeneity petersen2021wasserstein,van2025regression, xu2025wasserstein, song2026inference.
Third, we contribute to the literature on quantile models with endogeneity more broadly abadie2002instrumental,chernozhukov2005iv,chernozhukov2006instrumental,chernozhukov2008instrumental,lee2007endogeneity,lee2007nonparametric,frolich2013unconditional,kaplan2017smoothed,vuong2017counterfactual,decastro2019smoothed,wuthrich2019closed,wuthrich2020comparison,kaido2021decentralization,beyhum2023instrumental, holovchak2025distributional. The goal of this literature is to estimate the effects on the quantiles of real-valued scalar outcomes, whereas we consider the functional IV regression case in which the outcomes themselves are distribution-valued.
A last paper that does not fit neatly in these three groups is qu2024distributionally, who use the Wasserstein space to study IV methods, but focus on distributional robustness for classical IV assumptions rather than distributional treatment effects.
\paragraph{Notation.} We use the following notation throughout. For a vector $v \in \mathbb{R}^p$, $\|v\| = (v^{\!\mathsf{T}} v)^{1/2}$ denotes the Euclidean norm. For a positive semi-definite matrix $A$, $\|v\|_A = (v^{\!\mathsf{T}} A v)^{1/2}$ is the $A$-weighted norm. For functions $f \in L^2([0,1])$, $\|f\|_{L^2} = (\int_0^1 f(u)^2\, du)^{1/2}$ is the standard $L^2$ norm. The 2-Wasserstein distance between two distributions with quantile functions $Q_1, Q_2$ is $W_2(\mu,\nu) = \|Q_1 - Q_2\|_{L^2}$. For a measurable function \(f\), let \(P f:=\mathbb{E}[f(W)]\) and let \(\mathbb P_n f:=n^{-1}\sum_{j=1}^n f(W_j)\) denote the empirical measure applied to \(f\). We also write $\leadsto$ for weak convergence in the sense of vaart1996weak, $\overset{p}{\to}$ for convergence in probability, $\hat{\mathbb{G}}_\beta^*\rightsquigarrow_{\mathbb{P}}\mathbb{G}_\beta$ for conditional weak convergence in probability vaart1996weak: $\sup_{h\in BL_1}|E^x[h(\hat{\mathbb{G}}_\beta^*)]-E[h(\mathbb{G}_\beta)]| \overset{p}{\to} 0$, where $BL_1$ is the set of functions $\ell^\infty([a,b])^{p+1}\to\mathbb{R}$ with Lipschitz constant and supremum norm both bounded by one. Finally, we write $\ell^\infty(T)$ for the space of bounded functions on $T$ equipped with the supremum norm. For a closed convex set $K$ in a Hilbert space, $\Pi_K$ denotes the metric projection onto $K$.
Consider a setting with $n$ groups indexed by $j=1, \dots, n$. We are interested in estimating the effect of group-level variables $X_j \in \mathbb{R}^p$ with support $\mathcal{X} \subseteq \mathbb{R}^p$ on a distribution-valued outcome $Y_j\in \mathcal{Y}$, where $\mathcal{Y}$ is the space of one-dimensional cumulative distribution functions (CDFs) with finite variance. Let $Q_{Y_j} \in Q(\mathcal{Y})$ denote the quantile function corresponding to the CDF $Y_j$, where $Q(\mathcal{Y})$ is the space of quantile functions corresponding to $\mathcal{Y}$. Note that for a fixed quantile level $u\in (0,1)$, $Q_{Y_j}(u)$ is a real-valued random variable. We also have access to a vector of group-level IVs, which includes an intercept, $Z_j \in \mathbb{R}^{l+1}$, with $l\ge p$. We assume throughout that $\{(X_j,Y_j,Z_j)\}_{j=1}^n$ are sampled i.i.d. We fix $0 < a < b < 1$ and state all asymptotic results for quantile levels $u \in [a,b]$, avoiding tail quantiles.
We consider a linear model for $Q_{Y_j}$, noting that our asymptotic results will allow for misspecification,
where $\mu_X:=E[X_j]$, $\beta_0(u)$ is a scalar intercept function, $\beta_1(u)$ is a $p$-dimensional vector of slope functions, and $\eta_j(u)$ is an unobserved error term. Given that we include an intercept, the assumption that $E[\eta_j(u)]=0$ is a normalization. In what follows, we will often refer to $q(x,u):=\beta_0(u) + \beta_1(u)^{\!\mathsf{T}}(x-\mu_X)$ as the structural quantile function.
Writing the model in terms of demeaned variables $\tilde{X}_j=X_j-\mu_X$ will be convenient for our theoretical analysis.\footnote{When we use the demeaned coefficient and instrument vectors, we drop the intercept from those vectors to guarantee non-singularity of the covariance matrices.} Note that due to the linearity, we can always reparametrize the model and write it in terms of $X_j$ as
where $\tilde\beta_0(u)=\beta_0(u)-\beta_1(u)^{\!\mathsf{T}}\mu_X$. Finally, define $\mathbf{X}_j \coloneqq \left( 1, (X_j - \mu_X)^T\right)^T$.
In this section, we discuss identification and estimation in the \ensuremath{IVFR}\xspace model.
We obtain identification of the structural quantile function $q(x,u)$ using the vector of IVs, $Z_j$. To do so, we impose the following standard assumptions. Let $\tilde{Z}_j:=Z_j-E[Z_j]$.
This is the standard IV orthogonality condition for the structural error in (ref). Under misspecification, we define the pseudo-true coefficient functions by the population 2SLS projection given below; in general, the corresponding pseudo-residual need not satisfy $E[\tilde Z_j \xi_j(u)]=0$ in the overidentified case.
Define the population matrices $\Sigma_{ZZ} := E[\tilde Z_j\tilde Z_j^{\mathsf{T}}]$, $\Sigma_{ZX} := E[\tilde Z_j \tilde{X}_j^{\mathsf{T}}]$, and $\Sigma_{XX} := E[\tilde{X}_j \tilde{X}_j^{\!\mathsf{T}}]$.
To motivate \ensuremath{IVFR}\xspace, consider first the case where $X_j$ is exogenous, $E[\tilde{X}_j\eta_j(u)]=0$. In this case, we can use conventional global Fr\'echet regression petersen2019frechet to obtain the structural quantile function $q(x,u)$. petersen2019frechet motivate global Fr\'echet regression as a generalization of standard linear regression. Specifically, suppose that $Y_j\in \mathbb{R}$ and the conditional expectation function is linear, $m(x):=E\left[Y_j\mid X_j=x\right]=\alpha_0+\alpha_1^{\!\mathsf{T}}(x-\mu_X)$. Then, the formal characterization of the conditional expectation is $m(x)=\operatorname*{arg\,min}_{m\in \mathbb{R}}E\left[s(X_j,x)d_E^2(Y_j,m)\right]$, where $d_E(\cdot)$ is the standard Euclidean norm and $s(z,x)=1+(z-\mu_X)^{\!\mathsf{T}} \Sigma_{XX}^{-1}(x-\mu_X)$ are the linear regression weights. The idea of global Fr\'echet regression is to replace the Euclidean norm $d_E$ with a more general norm $d$ suitable for outcomes taking values in general metric spaces.
Here, we extend Fr\'echet regression to settings where $X_j$ is endogenous and the outcomes are distribution-valued. To motivate our approach, suppose first that $Y_j\in \mathbb{R}$ satisfies the linear IV model $Y_j=\alpha_0 + \alpha_1^{\!\mathsf{T}}(X_j-\mu_X) + \eta_j$ with $E[\eta_j]=0$. A key observation underlying \ensuremath{IVFR}\xspace is that the linear model $m(X_j):=\alpha_0 + \alpha_1^{\!\mathsf{T}}(X_j-\mu_X)$ can be obtained by rewriting the canonical two-stage least-squares (2SLS) estimator as $$ m(x)=\operatorname*{arg\,min}_{m \in \mathbb{R}} E\left[s(Z_j,x)d_E^2(Y_j,m) \right], $$ where
See the proof of Lemma (ref) in Appendix (ref) for a derivation.
Now, replace the Euclidean norm with the 2-Wasserstein distance $W_2$ suitable for distribution-valued outcomes, which for two one-dimensional distributions $Y_1, Y_2 \in \mathcal{Y}$ is defined as, \[ W_2(Y_1, Y_2) \coloneqq \left( \int_{0}^1 \left( Q_{Y_1}(u) - Q_{Y_2}(u) \right)^2 \mathrm{d} u \right)^{\frac12}. \]
Then, the above characterization of IV-Fr\'echet regression leads to the following instrumental-variables version of Fr\'echet regression,
where, by construction, for each given $x$, $m^{\ensuremath{\textnormal{IVFR}}\xspace}(x)$ is a distribution function, i.e., it lies in $\mathcal{Y}$.
It follows from the following Proposition that the quantile function associated with $m^{\ensuremath{\textnormal{IVFR}}\xspace}(x)$, denoted $Q_{m^{\ensuremath{\textnormal{IVFR}}\xspace}(x)}$, is the $L^2$-projection of the IV-weighted quantile function, $\psi_x(u) \coloneqq E[s(Z_j,x)Q_{Y_j}(u)]$, onto the space of quantile functions $Q(\mathcal{Y})$:
By the linearity of $s(Z_j,x)$ in $x$, the function $\psi_x(u)$ is affine in $x$ and can be written as $\psi_x(u) = \beta_0^{\text{unc}}(u) + \beta_1^{\text{unc}}(u)^{\!\mathsf{T}}(x - \mu_X)$, where $\beta^{\text{unc}}(u)$ are the standard 2SLS coefficients applied quantile by quantile. Also, define the pseudo-true residual as, \[ \xi_j(u):=Q_{Y_j}(u)-\mathbf X_j^{\!\mathsf{T}}\beta^{\mathrm{unc}}(u). \] By construction, $\beta^{\mathrm{unc}}(u)$ satisfies the population 2SLS normal equations. In the overidentified case, however, this does not generally imply $E[\tilde Z_j\xi_j(u)]=0$. A different GMM weighting matrix would in general define a different pseudo-true coefficient function. Throughout the paper, we focus on the 2SLS choice, but other linear GMM estimators could be used as well.
In population and under correct specification, the projection $\Pi_{\mathcal{Q}}$ is inactive and hence $Q_{m^{\ensuremath{\textnormal{IVFR}}\xspace}(x)}(u)$ coincides with the solution to the linear 2SLS regression quantile by quantile, i.e., $E[s(Z_j,x)Q_{Y_j}(u)]$. Moreover, the functional object $m^{\ensuremath{\textnormal{IVFR}}\xspace}(x)$ has the interpretation of being the (signed) IV-weighted Wasserstein barycenter of the group-level distributions $Y_j$. In other words, for a given $x$, $m^{\ensuremath{\textnormal{IVFR}}\xspace}(x)$ is the IV-weighted average of each group's distribution in probability space---which implies it can rightfully be called the “average” instrumented distribution. For more discussion on (conditional) Wasserstein barycenters, we refer to agueh2011barycenters, fan2024conditional, and panaretos2020invitation.
The \ensuremath{IVFR}\xspace coefficients are then equal to,
the OLS coefficients corresponding to the projected quantile function $Q_{m^\ensuremath{\textnormal{IVFR}}\xspace}(x)(u)$.
The following lemma shows that, under correct specification, $Q_{m^\ensuremath{\textnormal{IVFR}}\xspace(x)}(u)$ and $(\beta^\ensuremath{\textnormal{IVFR}}\xspace_0(u),$ $ \beta^\ensuremath{\textnormal{IVFR}}\xspace_1(u))$ recover the structural quantile function $q(x,u)$ and coefficients $(\beta_0(u),\beta_1(u))$, respectively.
Under correct specification, the projection is inactive and $\psi_x(\cdot)=q(x,\cdot)$ is a valid quantile function for all $x\in\mathcal X$. Hence $\Pi_{\mathcal Q}(\psi_x)=\psi_x$ and the Fr\'echet IV problem (ref) recovers the structural quantile function exactly.
Under misspecification, $\psi_x(\cdot)$ may violate monotonicity for some $x \in \mathcal{X}$, so $\Pi_{\mathcal{Q}}(\psi_x)$ may differ from $\psi_x$ and consequently $\beta^{\ensuremath{\textnormal{IVFR}}\xspace}(u)$ may differ from $\beta^{\text{unc}}(u)$. The \ensuremath{IVFR}\xspace coefficients remain well-defined and interpretable: for each $x$, the Fr\'echet IV solution $\Pi_{\mathcal{Q}}(\psi_x)$ is the closest valid probability distribution to $\psi_x$ in $W_2$-distance, and $\beta^{\ensuremath{\textnormal{IVFR}}\xspace}(u)$ parametrizes the best linear approximation to these projected distributions.
Here, we discuss the interpretation of the coefficient $\beta_1(u)$. The formal results underlying this discussion are presented in Appendix (ref). The coefficient $\beta_1(u)$ in the group-level quantile regression identifies the causal effect of group-level treatments on the group outcome distribution. When treatment is assigned at the group level, this group-level effect is the natural causal parameter. It captures the total response of the group’s distribution to the treatment, incorporating all channels through which the treatment operates. Note that this causal interpretation does not require any rank invariance assumption, as it operates at the group and not at the individual level.
The group-level model (ref) is agnostic about the within-group structure. The group quantile function $Q_{Y_j}(u)$ is the primitive; no assumptions are made about what generates it within the group. In particular, the within-group distribution can arise from an arbitrary mixture of individual-level outcomes, and the composition of the group may itself respond to treatment. The parameter $\beta_1(u)$ captures all these channels.
To illustrate, suppose that the treatment of interest is binary. Using potential outcomes notation, denote the counterfactual group-level CDFs without and with treatment as $Y(0)$ and $ Y(1)$, respectively. Then, under exogeneity and assuming the linear model holds, the causal parameter \ensuremath{IVFR}\xspace targets is the total causal effect, \[ \beta_1(u) = E[Q_{Y(1)}(u) - Q_{Y(0)}(u)], \] the average of the group-specific quantile treatment effects $Q_{Y(1)}(u) - Q_{Y(0)}(u)$. This object has a direct counterfactual interpretation: it is the difference between the group quantile functions with and without the treatment, averaging over group-specific responses. Thus, the parameter $\beta_1(u)$ is a natural analog of the average treatment effect in the standard setting where the outcome is scalar-valued.
This is different from the direct effect denoted by $\delta(u)$, which we show in Appendix (ref) is the effect identified by approaches that control for individual-level covariates within groups (such as CLP\xspace and MP\xspace with individual controls). The direct effect $\delta(u)$ holds the group composition (e.g., the share of low-skilled and high-skilled workers) constant and measures only the within-type response (e.g., the wage response of low-skilled vs.\ high-skilled workers). By contrast, the total effect $\beta(u)$ captures both the within-type response and the compositional response, capturing changes in who is in the group. For the direct effect to identify a causal effect on individuals, a rank invariance assumption is needed.
The total causal effect is a policy-relevant parameter in many applications. If a policymaker evaluates the impact of import competition on the local wage distribution, they typically care about what actually happens to wages in the affected regions, not a hypothetical scenario where wages change but workers cannot move. Indeed, the latter is often not a feasible policy. The total effect captures the equilibrium displacement of the distribution, and correctly characterizes the counterfactual effect on the treated unit, i.e., the entire group.
That said, the direct effect $\delta(u)$ can be informative when the goal is to isolate the within-type response, for instance in decomposition exercises that aim to separate wage structure changes from compositional shifts chernozhukov2013inference. The distinction between $\beta_1(u)$ and $\delta(u)$ is related to the distinction between unconditional and conditional quantile effects in the distributional treatment effects literature firpo2009unconditional, frolich2013unconditional; see also Remark 1 in MP\xspace for a related discussion.
We now compare $\ensuremath{\textnormal{IVFR}}\xspace$ to the estimators of CLP\xspace and MP\xspace. There are two cases.
\paragraph{Without individual-level covariates.} When CLP\xspace and MP\xspace do not control for individual characteristics, their estimand coincides with ours: both identify the total effect $\beta(u)$. In this case, working directly with random distributions at the group level, as the proposed \ensuremath{IVFR}\xspace estimator does, has several advantages. First, it allows us to avoid imposing any restrictions at the individual level, requiring only that the group quantile functions are linear in group-level covariates, while the individual-level structure can be fully nonparametric. Moreover, our framework also explicitly allows for misspecification of this group-level linear structure. By contrast, the statistical results in CLP\xspace and MP\xspace depend on the linear structure at the individual as well as group level. Second, our Fr\'echet regression formulation clarifies that estimation can proceed in a single step: regressing the group-level quantile functions on the treatment. Third, by projecting the quantile functions and recovering coefficients by OLS, our estimator can exploit the functional nature of the data to improve precision, as formally shown in Theorem (ref) below. Fourth, our estimator for $\beta(u)$ is guaranteed to produce valid quantile functions at every observed covariate value, while CLP\xspace and MP\xspace are not.\footnote{One could, alternatively, achieve the last two points by using monotone rearrangement instead of projection chernozhukov2010quantile. However, unlike projection, rearrangement does not naturally arise in the Fr\'echet regression framework, and one would lose the (IV-weighted) barycenter interpretation.}
Indeed, the simulations of Section (ref) show projection gains of up to 63% in IMSE, driven by finite-sample non-monotonicity in the fitted quantile functions. In our replication of CLP\xspace's empirical application in Section (ref), the fitted conditional wage quantile functions implied by unprojected \ensuremath{IVFR}\xspace violate monotonicity in around 2% of CZ--decade cells, so the projected point estimates are close to the unprojected ones. However, the projected bootstrap still delivers meaningfully tighter confidence bands.
\paragraph{With individual-level covariates.} When CLP\xspace and MP\xspace control for individual characteristics $W_{ij}$, they identify the direct (within-type) effect $\delta(u)$. As discussed above, these are different estimands, and which estimand is of interest depends on the application and the question being asked.
Finally, we emphasize that, in the presence of individual-level covariates, the instrument exogeneity assumption in CLP\xspace and MP\xspace is not weaker than the analogous assumption in \ensuremath{IVFR}\xspace without individual-level covariates. When the instrument used to identify the group-level treatment varies at the group level, our condition $\mathbb{E}[\tilde Z_j\eta_j(u)] = 0$ is equivalent to the orthogonality condition between the instrument and the group-level unobservable imposed by CLP\xspace and MP\xspace. Individual-level covariates change the estimand from a total group-level effect to a conditional/direct effect but do not weaken the required exclusion restriction for the group-level IV variation. This differs from the standard logic that conditioning on covariates can make the instrument exogeneity assumption more credible. See Appendix (ref) for a formal derivation. Note that, due to its functional nature, \ensuremath{IVFR}\xspace does not support the inclusion of individual-level covariates. Formulating a version of distribution-valued regression with individual-level covariates is an interesting avenue for future research, see also the distribution-on-distribution regression approaches of oliva2013distribution and ghodrati2022distribution.
The \ensuremath{IVFR}\xspace estimator solves the sample analogue of the population Fr\'echet IV problem (ref) at each observed covariate value, and recovers the linear parametrization by OLS:
where $\hat{\mathbf{X}}$ is the $n \times (p+1)$ matrix with rows $(1, (X_j - \hat{\mu}_X)^{\!\mathsf{T}})$, and $\hat{\psi}_x(u) := \frac{1}{n}\sum_{i=1}^n \hat{s}_i(Z_i, x)\, Q_{Y_i}(u)$ is the IV-weighted average quantile curve at $x$, with plug-in IV weights $$ \hat{s}_{j}(Z_j,x) := 1 + (x-\hat{\mu}_X)^{\!\mathsf{T}} (\hat{\Sigma}_{ZX}^{\!\mathsf{T}}\hat{\Sigma}_{ZZ}^{-1}\hat{\Sigma}_{ZX})^{-1} \hat{\Sigma}_{ZX}^{\!\mathsf{T}}\hat{\Sigma}_{ZZ}^{-1} (Z_j-\hat{\mu}_Z). $$
Since $\hat{\psi}_x(u)$ is linear in $x$, it defines the unprojected coefficient functions \[\tilde{\beta}(u) = (\tilde{\beta}_0(u), \tilde{\beta}_1(u)^{\!\mathsf{T}})^{\!\mathsf{T}}\] with $\hat{\psi}_x(u) = \tilde{\beta}_0(u) + \tilde{\beta}_1(u)^{\!\mathsf{T}}(x - \hat{\mu}_X)$. As above, because the IV weights can be negative, $\hat{\psi}_{X_j}(\cdot)$ may not be monotone. The projection $\Pi_{\mathcal{Q}}$ maps each IV-weighted average onto the space of valid quantile functions, yielding $\hat{Q}(X_j, u) := \Pi_{\mathcal{Q}}(\hat{\psi}_{X_j})(u)$---the closest valid probability distribution (in $W_2$) to the IV-weighted average at covariate value $X_j$. This is the sample IV-Fr\'echet barycenter, the sample analogue of the population projection (ref). It can be computed as an isotonic regression by the pool-adjacent-violators algorithm (PAVA) ayer1955empirical, miles1959complete, kruskal1964nonmetric. If only samples from the distribution $Y_j$ are available, we replace $Q_{Y_j}(u)$ by the empirical quantile function $\widehat{Q}_{Y_j}(u)$. We show in Section (ref) below that under a weak condition on the relationship between within- and across-group sample sizes, this additional first-stage estimation does not affect our asymptotic results.
The projection guarantees that each $\hat{Q}(X_j, \cdot)$ is a valid quantile function and is closer to the true structural quantile function $q(X_j, \cdot)$ than the unprojected $\hat{\psi}_{X_j}$ in all $L^p$ norms (Lemma (ref) below). The OLS step in (ref) then recovers the linear parametrization from the projected quantile functions.
In this section, we study the finite sample properties of the \ensuremath{IVFR}\xspace projection step and show that it decreases the estimation error. The first lemma shows that the projection step decreases the estimation error in the quantile function. This is a classical property of isotonic projections Robertson1988 and is reproduced from van2025regression.
Applied to \ensuremath{IVFR}\xspace, this result implies that at each observed $X_j$, the projected $\hat{Q}(X_j, \cdot) = \Pi_{\mathcal{Q}}(\hat{\psi}_{X_j})$ is closer to any valid quantile function than the unprojected $\hat{\psi}_{X_j}$, in every $L^p$ norm.
Next, we show that this improvement translates into finite sample improvements for the proposed estimator of the coefficient functions, which are often the main objects of interest. Throughout, let $b=(b_0,b_1)$ with $b_0\colon[0,1]\to\mathbb{R}$ and $b_1\colon[0,1]\to\mathbb{R}^p$ denote reference coefficients such that $q_b(X_j,\cdot):=b_0(\cdot)+b_1(\cdot)^{\!\mathsf{T}}(X_j-\hat\mu_X)\in\mathcal Q$ for all $j=1,\ldots,n$. The results below hold for any such $b$, and no assumption on the linear model (ref) is needed. In practice, $b$ is set equal to the estimation target. Under correct specification of (ref), setting $b_0(u)=\beta_0(u)+(\hat{\mu}_X-\mu_X)^{\!\mathsf{T}}\beta_1(u)$ and $b_1=\beta_1$ gives $q_b(X_j,u)=q(X_j,u)\in\mathcal Q$, so the bounds measure the estimation error for the true coefficients $(\beta_0,\beta_1)$. Under misspecification, a natural choice is the pseudo-true parameter $b = \beta^{\mathrm{unc}}$: under Assumption (ref) below, it produces strictly increasing quantile functions for all $x\in\mathcal X$ with slope at least $\kappa>0$, so $q_b(X_j,\cdot)\in\mathcal Q$ for all $j$ with probability approaching one as $n\to\infty$. In either case, the projection brings the \ensuremath{IVFR}\xspace coefficients closer to the target than their unprojected counterparts.
We emphasize that Theorem (ref) does not assume the linear model (ref). The bound in Theorem (ref) is a weighted joint improvement for estimating $(b_0, b_1)$. It does not guarantee an improvement for estimating any $b_{1,k}$ separately. The following results provide conditions under which the slope coefficients also improve individually.
To isolate individual slope coefficients, we use the Frisch--Waugh--Lovell (FWL) decomposition and the variational characterization of the isotonic projection. Let $C$ be the $n\times p$ matrix with rows $(X_j-\hat\mu_X)^{\!\mathsf{T}}$, write $c_k$ for its $k$-th column and $C_{-k}$ for $C$ with column $k$ removed. Define the FWL residual $r_k:=M_{-k}c_k$ with $M_{-k}:=I_n-C_{-k}(C_{-k}^{\!\mathsf{T}} C_{-k})^{-1}C_{-k}^{\!\mathsf{T}}$, let $r_{jk}$ denote its $j$-th entry, and set $J_k:=\{j\le n:r_{jk}\neq 0\}$ and $\hat v_k:=\frac{1}{n}\sum_{j=1}^n r_{jk}^2>0$. Let $\pi_k$ denote the coefficient vector from the sample linear projection of $c_k$ on $C_{-k}$, so that $c_k=C_{-k}\pi_k+r_k$. For a given target $b$, define the nuisance estimation error
which collects all estimation errors except the error in the $k$-th slope along the residualized variation $r_{jk}$. This decomposition is justified by the algebraic identity,
which separates the unprojected fitted curve $\hat\psi_{X_j}$ into three components: the target quantile function $q_b(X_j,\cdot)$, the nuisance error $e_{jk}$, and the $k$-th slope error along the identifying variation $r_{jk}$.
The bound (ref) decomposes the change in $k$-th coefficient error into a projection gain $\frac{1}{n\hat v_k}\sum_{j\in J_k}\|D_{X_j}\|^2_{L^2}$, which is always non-negative, and a nuisance cross-term $\frac{2}{n\hat v_k}\sum_{j\in J_k}\langle D_{X_j}, e_{jk}\rangle_{L^2}$, which captures the interaction between projection corrections and nuisance estimation errors and can be negative or positive. An improvement obtains whenever the nuisance cross-term is positive, or when it is negative but does not overwhelm the projection gain. The following corollary gives a sufficient condition that eliminates the nuisance cross-term entirely, thus guaranteeing an improvement.
The decomposition (ref) clarifies what condition (ref) requires: the error in the $k$-th slope along the identifying variation $r_{jk}$ must be necessary for any monotonicity violations in $\hat\psi_{X_j}$. The nuisance errors $e_{jk}$ may erode the monotonicity margin of $q_b(X_j,\cdot)$, but as long as they do not destroy it entirely, every projection correction is fixing a violation attributable to the $\beta_{1,k}$ error. Put differently, (ref) requires $q_b(X_j,\cdot)$ to increase steeply enough to absorb whatever local decreases $e_{jk}$ introduces. Since every term in $e_{jk}$ is $O_p(n^{-1/2})$ uniformly in $u$ and $j$---the latter by Assumption (ref), which bounds $\|c_{j,-k}\|$---the condition holds asymptotically whenever the target quantile functions have slopes bounded away from zero. This is the same requirement needed for uniform convergence of the estimators, formalized in Assumption (ref) below. When $p=1$, there is no partialling out ($r_{j1}=X_j-\hat\mu_X$), the nuisance slope term in $e_{j1}$ vanishes, and (ref) reduces to $q_b(X_j,\cdot)\in\mathcal Q$: the assumed monotonicity of the target.
Theorem (ref) and Proposition (ref) sit at two ends of a spectrum. Theorem (ref) requires no additional assumptions beyond valid reference coefficients, but guarantees improvement only for all coefficients jointly. Proposition (ref) isolates improvements coefficient by coefficient, at the cost of requiring condition (ref) on which errors drive the monotonicity violations. A natural intermediate approach is to pick a subset $S$ of coefficients for which (ref) is expected to hold. The same argument then yields an improvement result for that subset, with the nuisance error collecting only estimation errors from the intercept and the slopes outside $S$.
In this section, we establish functional central limit theorems for the \ensuremath{IVFR}\xspace estimator under weak conditions. We do not require the linear model (ref) to hold in population. When the full group-level quantile functions are observed for each group, the unprojected results allow these quantile functions to be discrete. For the projected results, we impose below an average smoothness condition that rules out synchronized jumps at fixed quantile indices while still allowing individual group distributions to have atoms and flat parts. This contrasts with existing inferential results for global Wasserstein-Fr\'echet regression petersen2021wasserstein and grouped quantile methods (e.g., CLP\xspace; MP\xspace), which typically impose density or pathwise smoothness conditions. Throughout this section, we work on $[a,b] \subset (0,1)$, and the projection $\Pi_{\mathcal{Q}}$ should hence be understood with respect to the relevant cone of monotone functions in $L^2([a,b])$.
Assumption (ref) corresponds to the standard 4th moment condition in linear IV regression hansen2022econometrics.
To establish the weak convergence of the fully projected \ensuremath{IVFR}\xspace coefficients, we will additionally impose a boundedness restriction on the covariates.
Assumption (ref) requires i.i.d. sampling across groups, but importantly leaves the within-group sampling scheme fully unrestricted. Additionally, Assumptions (ref) and (ref) are implied by the assumptions of CLP\xspace and MP\xspace, who assume bounded support of both $X$ and $Z$.
The next assumption requires that $\psi_x(\cdot)$ is uniformly strictly increasing on $[a,b]$ with slope $\geq \kappa$, for all $x$ in the support.
Assumption (ref) is weaker than requiring every group to have a density bounded away from zero (as, e.g., in CLP\xspace). In particular, it allows for flat parts and atoms in individual group quantile functions. The assumption only requires that, on every quantile interval, a non-negligible IV-weighted share of groups contributes a positive increment. Technically, Assumption (ref) gives \(\psi_x\) a uniform monotonicity margin on \([a,b]\), so the projection is asymptotically inactive and the functional delta method applies. Note that Assumption (ref) does not require smoothness. It is implied by a fixed positive IV-weighted fraction of groups locally “moving” (heterogeneous supports/atoms). See Remark (ref) below for the case where Assumption (ref) fails.
For the projected asymptotic results, we also require a mild average smoothness condition on the cross-group average quantile process.
Assumption (ref) is an average continuity condition in the quantile index. It allows individual group distributions to have atoms, but rules out a positive mass of groups having a jump at the same quantile index. It is weaker than the pathwise Lipschitz conditions imposed in CLP\xspace on the underlying group-specific coefficient and error processes. Note that it is distinct from Assumption (ref). Assumption (ref) gives a lower bound on the IV-weighted population slope, while Assumption (ref) gives an upper bound.
With these assumptions in hand, we now first establish functional CLTs for the unprojected and projected estimators of the structural quantile function at a fixed $x$, respectively.
The next two theorems establish functional CLTs for the unconstrained and \ensuremath{IVFR}\xspace estimators of the coefficient vector, respectively.
Under these conditions, the \ensuremath{IVFR}\xspace estimator has the same asymptotic variance as its unprojected variant, while Theorem (ref) guarantees smaller finite-sample error.
Under the conditions of Theorem (ref), projected and unprojected \ensuremath{IVFR}\xspace share the limiting process $\mathbb G_\beta$. Therefore, one can use the unprojected influence functions as a starting point for developing inference procedures. We now describe pointwise and uniform procedures, distinguishing unprojected and projected variants. The uniform bootstrap confidence band results are new to the literature.
\paragraph{Pointwise inference.} Write $\bar{Q}_n(u):=n^{-1}\sum_{j=1}^n Q_{Y_j}(u)$ for the cross-group average quantile function. For unprojected estimator $\tilde\beta(u)$, the asymptotic variance at a fixed $u_0$ is $\Omega(u_0, u_0)/n$ from Theorem (ref). A consistent estimator is \[ \hat{\Omega}(u_0, u_0) = \hat{T} \left(\frac{1}{n}\sum_{j=1}^n \hat{\Phi}_j(u_0)\hat{\Phi}_j(u_0)^{\!\mathsf{T}}\right) \hat{T}^{\!\mathsf{T}}, \] where $\hat{T} = \mathrm{diag}(1, \hat{S}_{\mathrm{2SLS}})$, $\hat{S}_{\mathrm{2SLS}}$ is the sample analogue of $S_0$, and the score is \[ \hat{\Phi}_j(u) =
, \qquad \hat{\xi}_j(u) = Q_{Y_j}(u) - \hat{\mathbf{X}}_j^{\!\mathsf{T}}\tilde{\beta}(u). \] The sandwich pointwise $(1-\alpha)$ CI for $\beta_k(u_0)$ is $\tilde\beta_k(u_0)\pm z_{1-\alpha/2}\,\hat\sigma_k(u_0)/\sqrt{n}$, where $\hat\sigma_k^2(u_0):=\hat\Omega_{kk}(u_0,u_0)$.
For the \ensuremath{IVFR}\xspace estimator, the projected pointwise CI replaces $\tilde\beta$ with $\hat\beta^{\ensuremath{\textnormal{IVFR}}\xspace}$ and recomputes the residuals accordingly: define $\hat{\xi}_j^{\ensuremath{\textnormal{IVFR}}\xspace}(u):=Q_{Y_j}(u)-\hat{\mathbf{X}}_j^{\!\mathsf{T}}\hat\beta^{\ensuremath{\textnormal{IVFR}}\xspace}(u)$ and \[ \hat{\Phi}_j^{\ensuremath{IVFR}\xspace}(u) =
, \] with $\hat\Omega^{\ensuremath{\textnormal{IVFR}}\xspace}$ and $\hat\sigma_k^{\ensuremath{\textnormal{IVFR}}\xspace}$ defined analogously. Under the conditions of Theorem (ref), both variance estimators are consistent for $\Omega$. However, the projected residuals $\hat\xi^{\ensuremath{\textnormal{IVFR}}\xspace}_j$ are closer to the population residuals in finite samples (by Theorem (ref)), which can yield tighter pointwise CIs, as we demonstrate in the simulations and the empirical application below.
\paragraph{Uniform inference.} For uniform confidence bands over $u \in [a,b]$, we use a multiplier bootstrap. Let $\{\omega_j\}_{j=1}^n$ be i.i.d. multiplier weights independent of the data satisfying $E[\omega_j]=0$, $E[\omega_j^2]=1$, and $E|\omega_j|^{2+\delta}<\infty$ for some $\delta>0$. The unprojected bootstrap process is \[ \hat{\mathbb{G}}_\beta^*(u) = \hat{T}\,\frac{1}{\sqrt{n}}\sum_{j=1}^n \omega_j \hat{\Phi}_j(u). \] Throughout, $P^x$ and $E^x$ denote probability and expectation conditional on the data (i.e., with respect to the multiplier weights only). Recall that $\hat{\mathbb{G}}_\beta^*\rightsquigarrow_{\mathbb{P}}\mathbb{G}_\beta$ denotes conditional weak convergence in probability vaart1996weak. The following theorem establishes the conditional weak convergence of $\hat{\mathbb{G}}_\beta^*$.
The $(1-\alpha)$ unprojected uniform band for $\beta_k(u)$ is $\tilde{\beta}_k(u) \pm \hat{c}_{1-\alpha}\cdot\hat{\sigma}_k(u)/\sqrt{n}$, where $\hat{c}_{1-\alpha}$ is the $(1-\alpha)$-quantile of $\sup_{u\in[a,b]} |\hat{\mathbb{G}}_{\beta,k}^*(u)| / \hat{\sigma}_k(u)$ conditional on the data. In practice, we sample $\omega_j \sim N(0,1)$, which satisfies the general conditions on the weights and is the standard approach.
The projected bootstrap constructs the bootstrap process differently: for each draw of multiplier weights, it computes bootstrap unconstrained coefficients $\tilde\beta^*(u)$, evaluates $\hat\psi^*_{X_j}(u):=\tilde\beta^*_0(u)+\tilde\beta^*_1(u)^{\!\mathsf{T}}(X_j-\hat\mu_X)$ at each observed $X_j$, applies $\Pi_{\mathcal{Q}}$ to each curve, and recovers projected bootstrap coefficients $\hat\beta^{\ensuremath{\textnormal{IVFR}}\xspace,*}(u)$ according to Eq.\ (ref). The projected bootstrap process is $\hat{\mathbb{G}}_\beta^{\ensuremath{\textnormal{IVFR}}\xspace,*}(u):=\sqrt{n}(\hat\beta^{\ensuremath{\textnormal{IVFR}}\xspace,*}(u)-\hat\beta^{\ensuremath{\textnormal{IVFR}}\xspace}(u))$, and the projected uniform band for $\beta_k^{\ensuremath{\textnormal{IVFR}}\xspace}(u)$ is constructed as above but with $\hat{\mathbb{G}}_\beta^{\ensuremath{\textnormal{IVFR}}\xspace,*}$ replacing $\hat{\mathbb{G}}_\beta^*$, $\hat\beta^{\ensuremath{\textnormal{IVFR}}\xspace}$ replacing $\tilde\beta$, and $\hat\sigma_k^{\ensuremath{\textnormal{IVFR}}\xspace}$ replacing $\hat\sigma_k$. As in the pointwise case, the projected bootstrap generally yields tighter pointwise CIs, as confirmed by the simulations and applications below. The following corollary establishes its validity.
In the previous sections, we have treated the group-level quantile functions $Q_{Y_j}$ as observed objects. This is suitable for settings with population-level within-group data (e.g., when all workers in a commuting zone or firm are observed). Here, we consider the case where researchers only observe a sample of the within-group realizations.
Suppose that for each group $j$, we observe draws $V_{jk}$, $k=1,\ldots,m_j$, from the distribution $Y_j$, which importantly need not be sampled i.i.d. The empirical distribution function and the associated empirical quantile function are \[ \widehat Y_j(y):=\frac{1}{m_j}\sum_{k=1}^{m_j}\mathbf{1}\{V_{jk}\le y\}, \qquad \widehat Q_{Y_j}(u):=\inf\{y:\widehat Y_j(y)\ge u\}, \quad u\in[a,b]. \] The feasible IV-weighted quantile curve is then \[ \bar\psi_x(u):=\frac{1}{n}\sum_{j=1}^n \hat s_j(Z_j,x)\,\widehat Q_{Y_j}(u), \] and we let $\bar{\tilde\beta}(u)=(\bar{\tilde\beta}_0(u),\bar{\tilde\beta}_1(u)^{\!\mathsf{T}})^{\!\mathsf{T}}$ denote the corresponding unconstrained coefficient functions defined by \[ \bar\psi_x(u)=\bar{\tilde\beta}_0(u)+\bar{\tilde\beta}_1(u)^{\!\mathsf{T}}(x-\hat\mu_X). \] The projected feasible estimator is \[ \bar\beta^{\ensuremath{\textnormal{IVFR}}\xspace}(u) := (\hat{\mathbf X}^{\!\mathsf{T}}\hat{\mathbf X})^{-1}\hat{\mathbf X}^{\!\mathsf{T}}\Pi_{\mathcal Q}(\bar\psi_{X_j})(u), \] that is, the estimator obtained from (ref) after replacing each $Q_{Y_j}$ by $\widehat Q_{Y_j}$. All other objects in Sections (ref)--(ref) are modified analogously.
Under the following high-level condition, this additional within-group sampling step is asymptotically negligible.
Assumption (ref) is stated in two parts. The first part requires that the stochastic fluctuations of the empirical quantile estimator are controlled uniformly over groups by a deterministic rate $\max_j r_{m_j}$, which follows from standard results under appropriate uniformity conditions, see Remark (ref) below. Note that this allows for the within-group sampling to exhibit dependence by using appropriate quantile function estimators rio2017asymptotic. Note that $Q_{Y_j}$ is itself a random target, but convergence here is conditional on a given group $j$. The randomness of $Q_{Y_j}$ across groups is accounted for in the across-group asymptotics above, and Assumption (ref) only controls the additional within-group estimation error $\widehat Q_{Y_j}-Q_{Y_j}$. The second part is a purely arithmetic growth condition ensuring that rate is fast enough relative to $n$. Together they imply $\sqrt{n}R_n = o_p(1)$, i.e., the first-stage quantile estimation error is uniformly smaller than the $\sqrt{n}$ rate governing the across-group IV problem.
The maximum over $j$ appears because Proposition (ref) below requires control uniformly across all groups. For instance, if $m_{\min}:=\min_{1\le j\le n} m_j$ and the empirical quantiles satisfy \[ R_n = O_p(m_{\min}^{-1/2}), \] then Assumption (ref) reduces to the sufficient condition $n/m_{\min}\to 0$. More generally, any quantile estimator whose uniform error over groups is $o_p(n^{-1/2})$ can be used.
The next proposition shows that, under Assumption (ref), replacing the true group quantile functions by empirical ones has no first-order effect.
As an immediate consequence, all of the limit theories in Section (ref) continue to hold verbatim for the feasible estimators.
To empirically validate the theoretical results above, we now compare unprojected and projected \ensuremath{IVFR}\xspace across several DGP configurations. Unprojected \ensuremath{IVFR}\xspace coincides with the CLP\xspace and MP\xspace estimators without individual covariates, though algorithmically our estimator, by working directly with the quantile functions as outcomes, does not require first-stage quantile regression. For each, we report coefficient integrated mean squared error (IMSE), the average squared Wasserstein distance $W_2^2$ between estimated and true conditional quantile functions, (uniform) confidence band coverage, and computational benchmarks. Across the reported configurations, the IMSE improves by up to 63%. Our uniform and projected bootstraps also achieve correct nominal coverage with up to 1.4% smaller width in the simulations (10% in the empirical application). Finally, computational benchmarks show our one-step estimator is faster than the two-step estimators of CLP\xspace and MP\xspace without individual covariates (see Appendix (ref)).
The simulations use the grouped-data IV design in CLP\xspace as the benchmark. Given our focus on unconditional group-level treatment effect, we drop CLP\xspace's individual-level covariate term (see the discussion in Section (ref)). The data-generating process (DGP) is
where $U_{ij}\sim U(0,1)$, $\zeta_j\sim U(0,1)$ and \[ X_j^0=(\pi_Z+\delta(\zeta_j-\tfrac{1}{2}))Z_j+\zeta_j+\nu_j, \qquad Z_j,\nu_j\sim \exp(0.25\,\mathcal{N}(0,1)). \] The terms $U_{ij}/2$ and $(\zeta_j-1/2)U_{ij}$ are written separately to match CLP's location term and group error. Together they equal $\zeta_j U_{ij}$.
Panel D in Table (ref) is the exact no-individual-covariate CLP benchmark. If $b_0\equiv0$, $\gamma(u)=\sqrt u$, $\pi_Z=1$, $\delta=0$, $p=1$, and $X_j=X_j^0$, then \[ y_{ij}=X_j\sqrt{u_{ij}}+\zeta_j u_{ij},\qquad X_j=Z_j+\zeta_j+\nu_j, \] which gives the DGP in CLP Appendix A with the individual covariate term omitted. The derivative of this population quantile function is $X_j/(2\sqrt u)+\zeta_j$, so the benchmark is strongly monotone because $X_j>0$ and $\zeta_j\ge0$.
Panels A--C keep the same first-stage endogeneity but add features that are common in applications: moderate instruments, controls, heterogeneous first stages, and treatment effects that need not be monotone in $u$. In these panels we standardize the endogenous first stage, $\bar X_j^0=(X_j^0-\bar X^0)/s_{X^0}$, and set $X_j=B_X\tanh(\bar X_j^0/B_X)$ with $B_X=1.5$ to make sure $X$ has bounded support without artificially clipping the range. Controls, when included, are bounded in the same way: $W_{jk}=B_W\tanh(\tilde W_{jk}/B_W)$ with $B_W=1.5$. Panel B sets their true effects to zero, so it isolates the finite-sample cost of estimating extra nuisance coefficient functions. Panel C uses mixed control effects: odd controls have $h_k(u)=u$ and even controls have $h_k(u)=0.1\sin(k\pi u)$. The treatment coefficient is \[ \gamma(u)=\sqrt{u}+\beta_{\text{slope}}\sin(2\pi u), \] so $\beta_{\text{slope}}=0$ gives the CLP coefficient and $\beta_{\text{slope}}>0$ adds curvature.
The extra baseline term $b_0(u)=\sigma q_0(u)$ lets us vary the shape of the within-group outcome distribution while preserving valid population quantile functions. For Panels A--C, we choose $\sigma$ analytically. Let $m_0=\inf_{u\in[a,b]}q_0'(u)$, $L_\gamma=\sup_{u\in[a,b]}|\gamma'(u)|$, and $L_k=\sup_{u\in[a,b]}|h_k'(u)|$. We set \[ \sigma m_0>B_XL_\gamma+B_W\sum_k L_k. \] Then every population quantile function generated by the DGP is strictly increasing on $[a,b]$. Any monotonicity violations reported below are therefore finite-sample violations in the fitted 2SLS curves, not failures of the population DGP.
Table (ref) reports coefficient IMSE, $W_2^2$, and the fraction of groups with non-monotone fitted curves, all averaged over 500 replications with 19 quantile grid points from $u = 0.05$ to $u = 0.95$. Unless otherwise noted, the designs use $n = 50$ groups and $N = 50$ observations per group.
\paragraph{Instrument strength (Panel A).} With a lognormal base distribution, a single endogenous regressor, and no controls, the largest gains occur when the first stage is weak or moderate. At $F \approx 5$, \ensuremath{IVFR}\xspace reduces IMSE by 63% and $W_2^2$ by 47%, with 15% of groups having non-monotone fitted distributions. At $F \approx 10$ the corresponding gains fall to 8% and 5%. With a stronger first stage ($F \approx 21$), violations are rare and gains are close to zero.
\paragraph{Controls and heterogeneous first stage (Panel B).} Panel B adds controls to the estimating equation while setting their true direct effects to zero. This isolates the finite-sample cost of estimating more nuisance coefficient functions. With $\delta=0$, increasing $p$ from 1 to 5 raises the share of non-monotone fitted curves from 1.4% to 2.6%, while projection gains remain small. With first-stage heterogeneity ($\delta=0.5$), violations rise from 4.7% to 6.5%. The projection then matters for estimation: with $p=2$, IMSE falls by 29% and $W_2^2$ by 20%; with $p=5$, the reductions are 25% and 15%. The gains need not be monotone in $p$, but the invalid-rate column shows the dimensionality channel directly. Each additional regressor adds a dimension to the fitted curve $\hat{\psi}_{X_j}(u)=\tilde{\beta}_0(u)+\sum_k\tilde{\beta}_{1,k}(u)(X_{jk}-\hat\mu_k)$, making it harder for all dimensions to satisfy monotonicity at once.
\paragraph{Realistic combination (Panel C).} Panel C combines small cells ($N=25$), $p=5$ controls, a heterogeneous first stage ($\delta=1$), lognormal outcomes, treatment-effect curvature ($\beta_{\text{slope}}=0.2$), and median first-stage $F\approx 11$. This conservative stress test gives a 20% IMSE reduction and a 12% $W_2^2$ reduction. About 13% of fitted quantile functions are non-monotone under unprojected \ensuremath{IVFR}\xspace.
\paragraph{CLP benchmark (Panel D).} Panel D is the exact no-individual-covariate CLP parameter point described above. In this all-positive-treatment design, population monotonicity is very strong: $Q_{Y_j}'(u)=X_j/(2\sqrt{u})+\zeta_j>0$. Projection gains are therefore concentrated at the smallest group counts, where sampling noise and weak first stages still create occasional fitted-curve violations. At $(n,N)=(25,25)$, IMSE falls by 17%; by $(n,N)=(50,50)$ the gain is essentially zero.
Theorem (ref) establishes that projected and unprojected \ensuremath{IVFR}\xspace share the same asymptotic Gaussian process under the stated assumptions, so standard inference (sandwich SEs, multiplier bootstrap) applies without modification.
Table (ref) reports pointwise and uniform coverage alongside band widths from 500 replications with $B = 500$ bootstrap draws. For unprojected \ensuremath{IVFR}\xspace we use sandwich SEs and the unprojected bootstrap; for projected \ensuremath{IVFR}\xspace we use the projected bootstrap, which runs PAVA inside each bootstrap draw. The first three rows use the simple $p=1$ design and vary $n$ and $N$ while keeping the median first-stage $F$ at about 10 or 20. The last three use the realistic design from Panel C and vary $n$ and $N$, with a median $F$-statistic of 11 and 23.
Coverage is near-nominal for both methods across all configurations, confirming the asymptotic equivalence in Theorem (ref). The projection does not distort inference. Because \ensuremath{IVFR}\xspace has lower IMSE at comparable standard errors, it has weakly higher power for testing $\beta_1(u) = 0$.
The projected bootstrap produces similar or narrower uniform bands for \ensuremath{IVFR}\xspace. In these coverage simulations, the reduction is below 1.5%; in the empirical CLP application it is about 10% (Section (ref)). The projected bootstrap captures the finite-sample regularization from the projection, yielding a weakly tighter critical value for the sup-statistic.
We present two applications. Section (ref) revisits the distributional effects of Chinese import competition on wages studied by autor2013china and CLP\xspace, illustrating the finite-sample gains from the \ensuremath{IVFR}\xspace projection and inference. Section (ref) studies the Food Stamp Program's effect on birth weights analyzed by almond2011inside and MP\xspace, using the decomposition from Appendix (ref) to compare the \ensuremath{IVFR}\xspace estimates to the conditional quantile treatment effect estimates.
We revisit the distributional effects of Chinese import competition on U.S.\ wages, following CLP\xspace and autor2013china (ADH henceforth). Their setting is well-suited for illustrating \ensuremath{IVFR}\xspace: there is one endogenous group-level treatment (import exposure), a standard instrument (other-country imports), and a distribution-valued outcome (the wage distribution within each commuting zone (CZ)).
\paragraph{Setting and data.} CLP\xspace estimate the effect of rising Chinese imports on the distribution of log weekly wages across U.S.\ commuting zones. The unit of observation is a CZ--decade pair ($j = 1, \ldots, 722$ CZs, two periods: 1990--2000 and 2000--2007, giving $n = 1{,}444$ observations). For each CZ--decade, the outcome is the vector of 19 within-CZ empirical quantiles $\hat Q_{Y_j}(u_q)$ for $u_q \in \{0.05, 0.10, \ldots, 0.95\}$. The endogenous variable $X_j$ is the change in Chinese import exposure per worker in commuting zone $j$. Following autor2013china, $Z_j$ is a shift-share instrument that assigns Chinese exports to other high-income countries to each CZ using lagged local industry employment shares, isolating China-side supply variation. The IV regression of $\hat Q_{Y_j}(u)$ on $X_j$ at the CZ level includes six continuous controls (manufacturing employment share, college share, foreign-born share, female employment share, routine task intensity, outsourcing exposure), eight census region dummies, and a period dummy, with the same controls included in the instrument vector $Z_j$. Regressions are weighted by CZ start-of-period population and standard errors are clustered by state (using a cluster multiplier bootstrap for the uniform bands). The first-stage $F$-statistic is 533, so the instrument is strongly relevant.
We apply \ensuremath{IVFR}\xspace to this data and compare the results with the original CLP\xspace quantile-by-quantile 2SLS estimates, which do not include individual-level covariates. In their replication package, CLP estimate a 2SLS regression with the decade-equivalent change in CZ wage quantiles as outcome. To match this approach, we project the level quantile functions in each period and then recover the slope coefficients by OLS on the long difference between these projected quantile functions. Since both estimators share the same asymptotic distribution under the conditions of Theorem (ref), any differences between them reflect the finite-sample effect of the \ensuremath{IVFR}\xspace projection step.
\paragraph{Results.} Figure (ref) overlays the CLP\xspace estimates (black circles) with the \ensuremath{IVFR}\xspace estimates (red triangles) for the full sample. Both series tell the same qualitative story: import competition depresses wages across the distribution, with the largest effects at lower quantiles (around $-1.4$ log points at the 10th quantile) and smaller effects in the upper tail ($-0.4$ to $-0.5$). The two coefficient functions are visually indistinguishable. The reason, as mentioned, is that CLP\xspace's specification is in long differences so the natural object whose validity the \ensuremath{IVFR}\xspace projection enforces is the implied level CZ wage quantile function for each (CZ, decade) cell. Anchoring those level quantile functions at the IPUMS-derived 1990 baseline quantile, fewer than $2\%$ violate monotonicity (0.0% in the full sample, 0.8% for females, 1.4% for males), and the projection alters the slope by at most $0.0012$ log points at any quantile. This does not mean the projection is irrelevant: applying it inside the bootstrap still smooths the finite-sample distribution of the coefficient process and tightens the confidence bands reported below.
Figure (ref) displays four layers of confidence bands around the \ensuremath{IVFR}\xspace estimates. The two red bands show pointwise 95% CIs: the outer (lighter) band uses cluster-robust sandwich SEs (identical to those of CLP\xspace); the inner (darker) band uses projected-bootstrap SEs that apply PAVA inside each multiplier-bootstrap draw. The two blue bands show cluster-robust uniform 95% bands over $u \in [0.05, 0.95]$: the outer (lighter) band uses the unprojected multiplier-bootstrap critical value; the inner (darker) band uses the projected-bootstrap critical value. The uniform bands are new---CLP\xspace did not provide a procedure to construct uniform confidence bands. The projected pointwise CIs are on average 9.4% narrower than the sandwich CIs (with reductions of up to 22% at quantiles where the unconstrained coefficient function is most non-monotone, e.g.\ $u = 0.85$); the projected uniform bands are 10.1% narrower at every $u$. Thus, the projection does meaningful empirical work for precision even though it does not change the CLP point estimate. Finally, we find that the unprojected multiplier-bootstrap SEs are within 1% of the sandwich SEs at every quantile, confirming that both estimate the same asymptotic variance.
The uniform bands refine CLP\xspace's conclusions about where effects are significant. CLP\xspace reported pointwise CIs, which reject the null of no effect at 12 of 19 quantiles. The projected pointwise CIs reject at 15 quantiles, illustrating the power gain coming from the narrower confidence bands. Pointwise inference, however, does not account for simultaneous testing across the quantile grid. The projected uniform band, which does, rejects the null at quantiles concentrated in the 10th--35th percentile range, plus $u = 0.75$. The very bottom of the distribution ($u = 0.05$), where the CLP\xspace pointwise CI barely excludes zero, does not survive the uniform correction. Above the median, effects are negative but imprecisely estimated. The projection adds modest power here as well: the unprojected uniform bands reject at 6 quantiles, and the projection flips $u = 0.75$ to significant. Import competition thus has its strongest and most robust effects on lower-middle wages, rather than at the very bottom of the distribution.
\paragraph{Subsampling exercise.} The full CLP\xspace sample has 722 CZs, causing unprojected and projected \ensuremath{IVFR}\xspace to produce near-identical point estimates, in line with the simulation evidence. As shown, the projection nonetheless still provides benefits in the form of tighter confidence bands. To further assess the practical gains from the projection step at smaller sample sizes, we subsample the data at various CZ counts and compare IMSEs.
Concretely, for each subsample size $M \in \{75, 100, 150, 200, 350, 500\}$, we draw 500 random subsets of $M$ commuting zones (keeping both decades for each CZ), re-estimate both unprojected and level-projected \ensuremath{IVFR}\xspace on each subsample, and compute the $\text{IMSE} = \int (\hat{\beta}_1(u) - \beta_1^*(u))^2\, du$, where $\beta_1^*(u)$ is the full-sample \ensuremath{IVFR}\xspace estimate. We decompose the IMSE into integrated squared bias and integrated variance.
Figure (ref) reports the results. The left panel shows the bias-variance decomposition; both methods are variance-dominated at every sample size, with squared bias accounting for less than 10% of total IMSE. The projection's effect is concentrated at the smallest sample sizes, where weaker first stages create larger fitted-CDF excursions for it to fix. At $M = 75$ (150 observations for 17 regressors), \ensuremath{IVFR}\xspace reduces IMSE by 14%, almost entirely through variance reduction. At $M = 100$ the gain shrinks to 1%, and at $M \ge 150$ the projection has essentially no IMSE effect. The simulation results in Table (ref) show the same pattern across DGPs: projection improves the IMSE materially when instruments are weak or fitted curves are oscillatory.
We now turn to the empirical application in MP\xspace, who study the effect of the Food Stamp Program (FSP) on birth weights. This setting provides a natural laboratory for the decomposition developed in Appendix (ref): the treatment (food stamps) is individually targeted but varies at the group level, creating a gap between the conditional quantile treatment effect $\delta(u)$ estimated by MP\xspace and the group-level distributional effect $\beta(u)$ estimated by \ensuremath{IVFR}\xspace.
\paragraph{Setting and data.} Following almond2011inside and MP\xspace, we study the staggered county-level introduction of the FSP between 1964 and 1975. Groups are county--trimester cells; within-group units are births in a county--trimester. The outcome is birth weight in grams. We use natality microdata from the NCHS (1968--1977) merged with county-level food stamp adoption dates from almond2011inside, restricting to births by Black mothers ($n \approx 2.8$ million individual births in approximately 17{,}000 county--trimester groups with at least 25 observations). Controls include per capita income, government transfers, and 1960 county characteristics interacted with a linear time trend. All regressions include county, state$\times$year, and trimester fixed effects, with standard errors clustered by county.
The treatment indicator $\mathit{fsp}_{ct}$ equals one if a food stamp program was in place at least three months before birth in county $c$ in trimester $t$. Identification in MP\xspace treats variation in FSP rollout as quasi-random, conditional on controls, so both estimators use $Z = X$.
\paragraph{Estimands.} MP\xspace estimate a conditional quantile model with individual-level covariates (child sex, mother's age and its square, legitimacy status),
where $\mathit{fsp}_{ct}$ was introduced above, $ x_{1ict}$ are individual-level controls, $x_{2ct}$ are the county-level controls mentioned above, and $\alpha(u, v_{ct})$ is the county-level unobservable. The coefficient $\delta(u)$ is the direct within-type effect: the shift in the $u$-th conditional quantile for a given type of mother. The \ensuremath{IVFR}\xspace estimand $\beta_1(u)$ targets the effect on the realized group quantile---the actual $u$-th percentile of birth weights in a county--trimester cell.
\paragraph{Results.} We replicate the MP\xspace estimates using their mdqr package mdqr and estimate $\beta_1(u)$ via \ensuremath{IVFR}\xspace on the same data. Figure (ref) shows the decomposition from equation (ref): \[ \beta_1(u) = \delta(u) + \underbrace{\mathbb{E}[\bar W_j(1) - \bar W_j(0)]'\gamma(u)}_{\text{composition}} + \underbrace{\mathbb{E}[\Delta_j(u;1) - \Delta_j(u;0)]}_{\text{re-ranking}}. \] The composition-fixed quantile $Q_{Y_j}^{\oplus}(u)$ is computed as the within-group mean of the first-stage fitted values from mdqr, and the re-ranking gap $\Delta_j(u) = Q_{Y_j}(u) - Q_{Y_j}^{\oplus}(u)$ is the difference between the realized group quantile and this average. Each component is then regressed on $\mathit{fsp}$ with the same fixed effects.
At the 5th percentile, $\delta(0.05) \approx 26$ grams (s.e.\ $9.5$): holding mother type fixed, FSP raises the conditional 5th percentile of birth weight by about 26 grams, a statistically significant effect.\footnote{Our point estimates of $\delta(u)$ are slightly smaller than those reported in MP\xspace---e.g., 26 vs.\ nearly 30 grams at the 5th percentile---most likely reflecting minor differences in the NCHS natality vintage and the county crosswalk used to merge births with the almond2011inside FSP rollout data, which leave us with 18{,}865 black county--trimester cells versus 19{,}482 in MP\xspace.} Yet $\beta_1(0.05) \approx -7$ grams (s.e.\ $10.8$)---the actual 5th percentile of the county birth weight distribution does not significantly move. The 33-gram gap is almost entirely accounted for by the re-ranking component, while the composition component is negligible.
The \ensuremath{IVFR}\xspace coefficient function $\beta_1(u)$ slopes upward: it is essentially zero at most quantiles and rises to roughly $+13$ grams at the 95th percentile. The pointwise 95% CI excludes zero at $u=0.95$ and is borderline at the next three grid points, consistent with FSP shifting the right tail of the birth weight distribution modestly upward. The uniform 95% confidence band, however, contains zero at every quantile. The pointwise significance at the top should therefore be read as suggestive evidence rather than a statistically robust finding at the conventional uniform level.
The large re-ranking term can be explained as follows. An individual at the group's 5th percentile is generally not at her own conditional 5th percentile. A young unmarried mother (whose child has lower baseline birth weight) sitting at the group's 5th percentile may be at, say, her conditional 15th percentile, where the treatment effect $\delta(0.15) \approx 6$ grams is far smaller than $\delta(0.05) \approx 26$ grams. The re-ranking term aggregates these within-type percentile shifts across all types at the group quantile cutoff: because $\delta(u)$ is steeply decreasing at the left tail, the effective treatment effect at the group's 5th percentile is a density-weighted average of $\delta$ evaluated at higher within-type ranks, where the effect is much smaller.
The composition channel---whether FSP changes who gives birth---is negligible throughout the distribution. This is consistent with food stamps affecting nutrition rather than fertility decisions.
While the model in (ref) imposes a common $\delta(u)$ across types, we can run the MP\xspace estimator separately within each of four demographic cells (mother age $<24$/$\geq 24$ $\times$ legitimate/illegitimate) to examine heterogeneity. Figure (ref) shows the type-specific conditional QTEs $\delta_k(u)$. At the 5th percentile, the effects range from $60$ grams for older married mothers to $-35$ grams for younger married mothers, with younger unmarried mothers---the dominant type at the group's left tail---showing an intermediate effect of approximately $40$ grams.
The overall $\delta(0.05) \approx 26$ grams is thus an average across heterogeneous type-specific effects, implicitly weighted by each type's conditional density at the group quantile cutoff. Types that are more concentrated at the left tail of the birth weight distribution---young unmarried mothers, who account for 46% of the density weight at the 5th percentile but only 36% of births---receive disproportionate weight. This density-weighting mechanism is inherent to any estimand based on conditional quantiles evaluated at a common quantile index.
\paragraph{Discussion.} The conditional quantile treatment effect is large and positive at the left tail, while the aggregate distributional effect is close to zero and statistically insignificant, because the types who drive the conditional estimate are not the types who are marginal at the group quantile cutoff. Young unmarried mothers---who comprise 36% of births but 46% of the probability density weight across types at the group's 5th percentile---drive the conditional estimate with a type-specific effect of approximately 40 grams, while the aggregate effect on the county's actual 5th percentile is indistinguishable from zero.
The two estimands answer different questions. The conditional effect $\delta(0.05) \approx 26$ grams documents that FSP delivers large nutritional benefits to the most vulnerable births, holding observed maternal type fixed. The \ensuremath{IVFR}\xspace estimate $\beta_1(0.05) \approx -7$ grams (insignificant) shows that these within-type gains do not translate to detectable shifts in the left tail of the county's birth-weight distribution. A county health department tracking the actual low-birth-weight rate---rather than conditional quantiles within demographic cells---would not detect a significant effect of FSP on the left tail of the birth weight distribution. This is because the babies at the county's 5th percentile are predominantly of young unmarried mothers sitting at much higher within-type ranks (around their 12th--15th conditional percentile), where the treatment effect $\delta$ is only 6--10 grams.
More broadly, this application illustrates when the two approaches diverge. For treatments that operate at the individual level but are identified through group-level variation (food stamps, school vouchers, Medicaid expansions), the conditional effect $\delta(u)$ captures the individual-level mechanism, while $\beta_1(u)$ captures the aggregate distributional impact. The re-ranking channel---which is first-order whenever $\delta(u)$ varies with $u$ and types are heterogeneously distributed across the group's outcome distribution---can dramatically attenuate the aggregate effect even when the conditional effect is large. For treatments that operate at the group level, like the import competition shock in Section (ref), the group quantile $\beta_1(u)$ is the natural estimand, and the decomposition is not needed.
This paper develops IV Fr\'echet regression (\ensuremath{IVFR}\xspace), a framework for estimating the effect of endogenous group-level treatments on distribution-valued outcomes. The approach recasts grouped quantile IV regression as an instrumental-variables problem in Wasserstein space. This perspective yields a simple estimator: construct IV-weighted average quantile curves, project them onto the space of valid quantile functions, and recover coefficient functions by OLS.
The paper makes three main contributions. First, it provides an identification result showing that, under standard quantile IV conditions, the structural distributional effect is the solution to an IV-weighted Fr\'echet problem. This gives the fitted object a clear interpretation as an instrumented average distribution. Second, it introduces a monotone projection step that guarantees valid fitted distributions and weakly improves finite-sample estimation error, while leaving the first-order asymptotic distribution unchanged under mild conditions. Third, it establishes functional asymptotic normality and multiplier-bootstrap procedures for pointwise and novel uniform inference over quantile indices.
Simulations and two empirical applications illustrate the practical value of the method. In finite samples, the projection can substantially reduce the integrated mean squared error (IMSE) relative to existing grouped quantile IV estimators. In an application to Chinese import competition, the method delivers tighter confidence bands and lower IMSE. Additionally, using our novel uniform confidence bands, we show that the evidence for wage losses is concentrated away from the very bottom of the distribution. In a second application to the effect of county-level food stamp programs on the birth weight distribution, we find no evidence for distributional effects using our uniform bands. More broadly, our results suggest that directly modeling outcomes as random distributions can sharpen both estimation and inference in settings where policy effects are inherently distributional.