The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
67,981 characters
Estimating Treatment Effects Under Bounded Heterogeneity
\author{Soonwoo Kwon\footnote{Department of Economics, Brown University, Email: soonwoo\[email removed]} \ and Liyang Sun\footnote{Department of Economics, University College London and CEMFI, Email: [email removed]}}
\date{March 2026}
\maketitle
\begin{abstract}
Specifications that impose constant treatment effects are common but biased,
while fully flexible alternatives can be imprecise or infeasible. Under a
bound on treatment effect heterogeneity, we propose a generalized ridge
estimator, \texttt{regulaTE}, that yields heterogeneity-aware confidence
intervals (CIs). The ridge penalty is chosen to optimally trade off worst-case
bias and variance in a Gaussian homoskedastic setting; the resulting CIs
remain tight more generally and are valid even
under lack of overlap. Varying the bound enables sensitivity analysis to departures from constant effects, which we illustrate in leading empirical applications of unconfoundedness and staggered adoption designs.
\end{abstract}
\smallskip
\noindent\texttt{Key words:} heterogeneous treatment effects, limited overlap, ridge regression
\smallskip
\noindent\texttt{JEL classification codes:} C12, C21, C51.
\newpage
\section{Introduction}
In many empirically relevant settings for causal inference, correctly specifying
the model requires allowing for fully flexible treatment effect heterogeneity as
a function of covariates. However, estimating such rich models often leads to
imprecise estimates, and can become infeasible under limited overlap in
covariate distributions. Consequently, estimating the fully flexible model may
not be desirable in empirical applications, particularly when researchers are
primarily interested in the average treatment effect (ATE) rather than the full
distribution of heterogeneous effects.
Consider the unconfoundedness design, where the treatment is assumed to be as
good as randomly assigned conditional on the covariates. Many applied studies
estimate a ``short'' regression of the form that assumes constant effects
\begin{equation*}
Y_i = D_i \beta_{\text{short}} + X_i'\gamma_{\text{short}} + u_i,
\end{equation*}
and interpret the estimated regression coefficient $\hat \beta_{\text{short}}$ as an estimate of the ATE \citep[e.g.,][]{martinez2014role,favara2015credit,michalopoulos2016long,bobonis2016monitoring,xu2018costs}.
There is growing recognition that treatment effect heterogeneity affects the interpretation of the ``short'' regression. In an influential study, \cite{angrist1998} shows that under unconfoundedness and discrete covariates, the ``short'' regression estimates a convex average of treatment effects associated with each covariate value, and therefore can differ from estimates based on matching that does not restrict treatment effect heterogeneity. \cite{sloczynski_interpreting_2022} further shows that this convex average corresponds to a weighted average of treatment effects on the treated and untreated, with counter-intuitive weights. Similar concerns have been raised in recent studies of the ``short'' regression for staggered adoption designs \citep{double_fe,goodman2021difference}.
However, revisiting \cite{angrist1998}, the textbook \cite[Chapter 3.3.1]{angrist2009mostly} asserts that if treatment effect heterogeneity is limited, the ``short'' regression estimate may still be preferable to unbiased but noisier semi- and nonparametric estimators of the average effect. First, when treatment effects are constant, the ``short'' regression is unbiased and most efficient under homoskedasticity. Second, \cite{angrist2009mostly} informally argue that the bias of the ``short'' regression remains small when treatment effects do not vary too much, and the resulting efficiency gain can still justify its use. An additional, and often overlooked, consideration is that in settings without overlap, semi- and nonparametric estimators may cease to be well defined, whereas the short regression remains feasible because it extrapolates treatment effect estimates to covariate values without overlap. Consequently, practitioners' intuition is that the confidence intervals (CIs) based on the ``short'' regression should have coverage close to the nominal level, making the ``short'' regression a pragmatic choice in applied work.
In this paper, we formalize the argument of \cite{angrist2009mostly}. Under a
restriction on treatment effect heterogeneity, specifically a bound on the variance of
conditional average treatment effect heterogeneity (VCATE), we first show how to
bias-correct the ``short'' regression CI such that the resulting CI has correct
coverage and is heterogeneity-aware. This approach is feasible even when
the average effect is not point identified due to a complete lack of
overlap, since the bound on treatment effect heterogeneity provides partial
identifying information on the treatment effect for the subsample without overlap.\footnote{When overlap fails completely, the ATE is not point identified. The VCATE bound provides partial identifying information that restricts the range of treatment effects in the non-overlapping sample, enabling valid inference even in this setting.} We denote this bound by $C^2$. As $C^2$ increases from zero, the associated bias-corrected
``short'' regression CI achieves correct coverage under each restriction,
allowing one to assess the sensitivity of the ``short'' regression estimates to
violations of constant effects.
However, as the ``short'' regression estimator does not adjust to the bound, bias correction can yield an excessively wide heterogeneity-aware CI. To provide a more informative sensitivity analysis, we propose a
generalized ridge regression, \texttt{regulaTE}, which coincides with the ``short'' regression when
heterogeneity is absent but otherwise can depart from the ``short'' regression if
its bias is large. We propose selecting the ridge penalty to minimize the heterogeneity-aware CI length under a normal, homoskedastic setting, which facilitates rapid computation. We show in simulations and empirical illustrations that under more general error distributions, the resulting CI remains tighter than the bias-corrected ``short'' regression CI.
We provide an implementation of sensitivity analysis based on \texttt{regulaTE} in our accompanying R
package.\footnote{An R package implementing the \texttt{regulaTE} estimator is
available online at \url{https://github.com/lsun20/regulaTEr}.}
An important input for our method is $C^{2}$, the bound on VCATE, which we primarily view as a sensitivity
parameter rather than a quantity to be estimated. In practice, the researcher
examines how the CIs change as $C^{2}$ varies over
a user-specified range, thereby assessing how sensitive the empirical conclusions
are to deviations from constant treatment effects. A useful summary is the
``breakdown value'' of $C^{2}$, the smallest value at which the ATE
becomes statistically insignificant, following \cite{kline2013sensitivity} and \cite{masten2020inference}. Because $C^{2}$ has a clear interpretation as a bound on VCATE, which is frequently used in variance decomposition and welfare analyses
\citep{kline2020leave, sanchezbecerra2023robust,
dechaisemartin2024estimatingtreatmenteffectheterogeneitysites} and about whose magnitude applied researchers often have intuition, the resulting sensitivity analysis is straightforward to
interpret. We conduct sensitivity analysis based on our proposed \texttt{regulaTE} in two
leading empirical applications under unconfoundedness and staggered
adoption. We show how subgroup analyses reported in the original empirical
studies can be used to inform a plausible range of bounds for the sensitivity
analysis. We emphasize, however, that choosing $C^{2}$ in a data-driven manner
generally leads either to invalid confidence intervals or to excessively wide
intervals, reflecting the impossibility results of \cite{armstrong_optimal_2018}.
\medskip
\noindent\textbf{Related literature. } Since the seminal work of \cite{angrist1998}, a growing body of literature has characterized the bias of the ``short'' regression estimator under treatment effect heterogeneity, including \cite*{gibbons_serrato_urbancic_jem2018}, \cite{sloczynski_interpreting_2022}, \cite*{goldsmithpinkham2024contamination}, \cite{double_fe}, \cite{goodman2021difference}. We show how to bias-correct the ``short'' regression under a restriction on treatment effect heterogeneity, and propose a sensitivity analysis to evaluate the impact of treatment effect heterogeneity on robust estimation of the average effect.
Several papers have proposed estimators for the average effect under restrictions on treatment effects. In the context of the unconfoundedness design, \cite{lechner2008note} assumes outcomes are bounded and derives the Manski bound for the ATE. \cite{lee2021bounding} improves the Manski bound by combining it with the inverse probability weighting (IPW) estimator using reference propensity scores. \cite*{athey_approximate_2018} assume sparsity of effects, \cite{armstrong2021finite} assume smoothness of effects with respect to covariates, and \cite{dechaisemartin2021tradingoff} impose bounds on effects. In the context of difference-in-differences, \cite{manski_how_2018} impose bounds on the variation in outcomes over time and across states. In contrast to these papers, our method is tied directly to a measure of treatment effect heterogeneity by restricting the variance of treatment effects. This does not entail a functional form assumption on the heterogeneity, mimicking the common belief that effects are broadly similar.
There are semi- and nonparametric estimators, such as the inverse propensity
score weighted estimator, that remain unbiased for the average treatment effect
even under unrestricted treatment effect heterogeneity. However, these
estimators rely on strong overlap for consistency and asymptotic normality. When
propensity scores approach 0 or 1, they can become highly variable, making the
conventional normal approximation unreliable, as demonstrated in
\cite{rothe2017robust} or \cite{heiler2021valid}. To address extreme values of propensity scores, \cite*{crump_dealing_2009} propose trimming the inverse propensity score weighted estimator by removing observations with estimated propensity scores outside the range $(0.1,0.9)$. Since this approach only estimates the average effect for a subpopulation with strong overlap, several papers have characterized the resulting bias and, under certain smoothness conditions, proposed bias correction methods; see, e.g., \cite{ma2020robust} and \cite{sasaki2022estimation}. Our method provides an alternative route for dealing with limited overlap, which includes settings with complete lack of overlap: rather than focusing on a subpopulation with strong overlap, we restrict treatment effect heterogeneity and use this restriction to infer the worst-case heterogeneity for covariate cells with limited overlap.
The remainder of the paper is organized as follows. Section~\ref{sec:setup} introduces the
setup and formalizes the restriction on treatment effect heterogeneity.
Section~\ref{sec:proposed-method} develops the proposed inference method and
illustrates its finite-sample performance in calibrated simulations. Section~\ref{sec: estimate VCATE}
applies the method to conduct sensitivity analysis for several empirical settings. Section~\ref{sec:conclusion} concludes. Proofs and comparison to the adaptive estimator of \cite{armstrong2023adapting} are found in the Appendix.
\section{Setup and treatment effect heterogeneity}\label{sec:setup}
Consider the unconfoundedness setting where we observe a random sample of $n$ units, and each unit is characterized by a pair of
potential outcomes $(Y_i(0), Y_i(1))$ under no treatment and treatment,
respectively; a set of (possibly including flexible transformations of) covariates $X_i \in \mathbb{R}^{k+1}$
that includes a constant; and a treatment indicator $D_i \in \{0,1\}$. We assume
unconfoundedness (also known as selection on observables):
${Y_i(0), Y_i(1)} \perp D_i \mid X_i$. As common in empirical practice \citep{wooldridge_econometric_2010}, we further assume that
$\ensuremath{\mathbb{E}}[Y_i(1) \mid X_i]$ and $\ensuremath{\mathbb{E}}[Y_i(0) \mid X_i]$ are both linear in $X_i$.\footnote{In principle, such linearity assumptions are not
required, and our approach to estimating average effects under bounded
treatment effect heterogeneity can be extended to more general settings. For
example, assuming linearity only of $\ensuremath{\mathbb{E}}[Y_i(0)\mid X_i]$ leads to a related
ridge-type estimation procedure. One could also adopt a fully nonparametric
approach by imposing a bound on treatment effect heterogeneity in addition to
the smoothness assumptions of \cite{armstrong2021finite}. Nevertheless, given the prevalence of linear models in empirical practice, and the intuitive interpretation and computational efficiency afforded by this assumption, we view this as a useful trade-off between generality and practicality.}
Let $\tau(x) = \ensuremath{\mathbb{E}}[Y_i(1)-Y_i(0)\mid X_i = x]$ denote the covariate-specific
average treatment effect (CATE). As discussed in Section~\ref{sec:challenge}, limited overlap remains a key empirical challenge even under the linearity assumption.
Unless stated otherwise, we condition on the realized values
$\{x_i, d_i\}_{i=1}^n$ of the covariates and treatment $\{X_i, D_i\}_{i=1}^n$,
where lowercase letters denote realizations. All probability statements are
therefore taken with respect to the conditional distribution of
$\{Y_i(0), Y_i(1)\}_{i=1}^n$. Under these assumptions, the DGP for the observed
outcome $Y_i$ can be written as
\begin{equation}
Y_i = d_i \tau(x_i) + x_i'\gamma + \varepsilon_i, \label{eq:dgp_selection}
\end{equation}
where $\{\varepsilon_i\}_{i=1}^n$ is a sequence
of mutually independent noise terms that are mean zero. This fixed-design framework follows, for example, \citet{rothe2017robust} and \citet{armstrong2021finite}, to investigate finite-sample properties. In Section~\ref{sec:application}, we extend to settings where the covariates that parameterize treatment effect heterogeneity differ from the confounders $x_i$, but focus first on the case where they coincide to describe the challenge cleanly.
Our goal is to estimate a weighted average of the CATEs with known weights,
where the weights reflect a predetermined target population. The leading example
is the ATE,
$\beta=\beta^{\text{ATE}}=\ensuremath{\mathbb{E}}_n[\tau(x_i)]=\frac{1}{n}\sum_{i=1}^n \tau(x_i)$.\footnote{This parameter is sometimes referred to as the conditional ATE (e.g., \citet{imbens2009recent}). Since our analysis is always conditional on the realized covariates, we simply refer to it as the ATE without confusion, and reserve the term CATE for the function
$\tau(x)$.} Here $\ensuremath{\mathbb{E}}_n$ denotes the ``empirical'' mean
(i.e., $\ensuremath{\mathbb{E}}_n a_i = n^{-1} \sum_{i=1}^n a_i$ for any fixed or random sequence
$\{a_i\}_{i=1}^n$). Under the linearity assumption, we can write
$\tau(x_i) = \alpha+x_{i,-1}'\delta$, where $x_{i,-1}$ denotes the vector of
covariates excluding the constant. Letting
$\widetilde x_{i,-1}=x_{i,-1} - \ensuremath{\mathbb{E}}_n[{x}_{i,-1}]$ denote the demeaned
covariates, deviations of the CATE from the ATE are given by
$\tau(x_i)-\beta^{\text{ATE}} = \widetilde x_{i,-1}'\delta$. Therefore, the DGP
\eqref{eq:dgp_selection} can be reparameterized as the \textit{long regression}
as in \citet{imbens2009recent}:
\begin{equation}
Y_i = d_i \beta^{\text{ATE}}_{\text{long}} + (d_i \tilde x_{i,-1} )'\delta + x_i'\gamma + \varepsilon_i,\ \text{where }\widetilde x_{i,-1}=x_{i,-1} - \ensuremath{\mathbb{E}}_n[{x}_{i,-1}]. \label{eq:long reg}
\end{equation}
Other leading target parameters include the average treatment effect on the
treated (ATT) and the average treatment effect on the untreated (ATU). Since the long regression~\eqref{eq:long reg} handles all three targets through appropriate demeaning of the covariates, we focus on the ATE for simplicity.
We define the short regression as the specification that omits the interaction terms between treatment and covariates:
\begin{equation}
Y_i = d_i \beta_{\text{short}} + x_i'\gamma_{\text{short}} + u_i, \label{eq:short reg}
\end{equation}
which is the most commonly used specification in empirical research.
\subsection{Challenges under treatment effect heterogeneity}\label{sec:challenge}
While the long regression~\eqref{eq:long reg} is the correct specification, estimating it can be practically challenging: the resulting estimates may be highly imprecise, and the model may not even be
estimable under limited overlap. First, when the covariates $x_i$ are mixed,
containing both continuous covariates and discrete covariates, the discrete
covariates typically enter the long regression as fixed effects capturing
confounding from, for example, location effects. If all units in a given
location are treated, contributing to lack of overlap, then the interaction
between treatment and the centered location fixed effect becomes multicollinear
with the treatment indicator and the location fixed effect itself, rendering the
long regression $\hat\beta_{\text{long}}$ undefined.
Second, if the covariates $\{x_i\}_{i=1}^n$ are generated by saturating discrete covariates, when overlap fails
entirely for a particular covariate value, the long regression $\hat\beta_{\text{long}}$ is again not well-defined. The long regression coefficient estimator is algebraically equivalent to the following IPW estimator:
\[
\hat\beta_{\text{long}}^{\text{ATE}}=\ensuremath{\mathbb{E}}_{n}\left[\frac{d_{i}-p(x_{i})}{p(x_{i})\left(1-p(x_{i})\right)}Y_{i}\right]
\]
where $p(x_i)$ denotes the empirical propensity score for covariate value $x_i$. Therefore, if $p(x_{i})\in\{0,1\}$ for some $x_{i}$ due to lack of overlap, then the long regression $\hat\beta_{\text{long}}$ is not well-defined. Furthermore, because $p(x_{i}) (1-p(x_{i}))$ appears in the denominator, the variability of
$\hat\beta_{\text{long}}^{\text{ATE}}$ is heavily influenced by observations for
which $p(x_{i})$ is close to zero or one.
In contrast, the short regression estimator is less affected by limited
overlap and remains well-defined even under complete lack of overlap at some covariate values. To see this, note that by the Frisch-Waugh-Lovell (FWL) theorem, the short regression
estimator in \eqref{eq:short reg} can be written as
\[
\hat\beta_{\text{short}}=\frac{\ensuremath{\mathbb{E}}_{n}\left[\left(d_{i}-p(x_{i})\right)Y_{i}\right]}{\ensuremath{\mathbb{E}}_{n}\left[\left(d_{i}-p(x_{i})\right)d_{i}\right]}.
\]
This representation highlights that the short regression places little weight on
observations with extreme propensity scores. Moreover, it remains well-defined
even when overlap fails for some covariate values, since such observations have
no influence on the estimator. However, by omitting interactions between $d_i$
and $\tilde{x}_{i,-1}$, the short regression is generally misspecified in the
presence of treatment effect heterogeneity. As a result, it may suffer from
omitted variable bias, with the bias measured relative to the ATE. From \cite{angrist1998}, the short regression estimand is:
\begin{equation}
\beta_{\text{short}}= \ensuremath{\mathbb{E}}[\hat{\beta}_{\text{short}}] = \ensuremath{\mathbb{E}}_{n}\left[\frac{p(x_{i})\left(1-p(x_{i})\right)}{\ensuremath{\mathbb{E}}_{n}\left[p(x_{i})\left(1-p(x_{i})\right)\right]}\tau(x_{i})\right]. \label{eq:short estimand}
\end{equation}
This estimand coincides with $\beta$ only if one of the following holds: (i)
$\tau(x_i)$ is constant, (ii) $p(x_i)$ is constant, or (iii) $p(x_i)(1-p(x_i))$
is uncorrelated with $\tau(x_i)$. In general, the short regression recovers a weighted average of heterogeneous treatment effects, with weights proportional to the conditional variance of treatment assignment. As discussed in \citet{poirier2024quantifyinginternalvalidityweighted}, the subpopulation for which the short-regression estimand corresponds to a true average treatment effect can be small, limiting the internal validity of the estimator.
Nevertheless, when treatment effect heterogeneity is limited, the short
regression remains nearly unbiased and typically yields more precise
estimates. This point is made informally by \citet[Chapter
3.3.1]{angrist2009mostly}: ``Of course, the difference in weighting schemes is of
little importance if $\tau(x)$ does not vary across cells (though weighting
still affects the statistical efficiency of estimators).'' This observation
helps explain why the short regression is substantially more common than the
long regression in empirical practice.
\subsection{Bounding the treatment effect heterogeneity via VCATE}\label{sec:vcate_bound}
To formalize the belief that treatment effects are broadly similar across
covariates, we impose an upper bound on the (sample) variance of the CATE
$\tau(x_i)$, the VCATE, which we define as
$\ensuremath{\mathbb{E}}_n [(\tau(x_i) - \ensuremath{\mathbb{E}}_n [\tau(x_i)])^2]$.
While any restriction that imposes that the vector of heterogeneous effects
$\{\tau(x_i)-\ensuremath{\mathbb{E}}_n[\tau(x_i)]\}_{i=1}^n$ lies in a bounded set would, in
principle, suffice to control the bias induced by misspecification of the short
regression or other linear-in-outcome estimators, we focus on bounds on VCATE
for several reasons. First, because variance is the most widely used and
canonical measure of dispersion, VCATE provides an intuitive measure of
treatment effect heterogeneity that captures how dispersed treatment effects are
across covariate values, rather than imposing potentially opaque geometric
restrictions on the entire vector of effects.\footnote{For example,
\cite{levy2021fundamental} and \cite{sanchezbecerra2023robust} introduce VCATE
as leading measures of treatment effect heterogeneity.} This makes the
restriction easy to interpret and to vary in sensitivity analyses. In contrast,
many alternative convex restrictions, such as $\ell_p$ balls, do not correspond to commonly interpreted notions of treatment effect
heterogeneity and are therefore harder to calibrate in empirical applications.
Second, VCATE is a well-studied object in the treatment effects literature, and
as a result researchers often have a clear sense of how to reason about its
magnitude.\footnote{See \citet{kline2020leave},
\citet{sanchezbecerra2023robust}, and
\citet{dechaisemartin2024estimatingtreatmenteffectheterogeneitysites} for
recent work on estimation and inference for VCATE.} VCATE arises naturally in
variance decompositions and distributional analyses, and recent work has
developed estimators and inference procedures for VCATE under a range of
identifying assumptions. This existing body of work makes VCATE a familiar and
interpretable sensitivity parameter: researchers can draw on prior results,
empirical benchmarks, and domain knowledge to assess whether a given bound on
VCATE is plausible in their setting. Moreover, Popoviciu's inequality gives a
simple (yet conservative) bound
\[
\ensuremath{\mathbb{E}}_n [(\tau(x_i) - \ensuremath{\mathbb{E}}_n [\tau(x_i)])^2] \le \frac{1}{4}(\max_i \tau(x_{i})-\min_i \tau(x_{i}))^2,
\]
providing a simple reference point for calibration.
Note that, under our maintained specification, VCATE can be expressed as a simple
quadratic form in the heterogeneity coefficient $\delta$:
\begin{equation}
\ensuremath{\mathbb{E}}_n\!\left[(\tau(x_i)-\ensuremath{\mathbb{E}}_n[\tau(x_i)])^2\right]
= \ensuremath{\mathbb{E}}_n\!\left[(\tilde{x}_{i,-1}'\delta)^2\right]
= \delta'\,\ensuremath{\mathbb{E}}_n[\tilde{x}_{i,-1}\tilde{x}_{i,-1}']\,\delta.\label{eq:wtd quadratic restriction}
\end{equation}
For notational convenience, define
\(
V_x := \ensuremath{\mathbb{E}}_n[\tilde{x}_{i,-1}\tilde{x}_{i,-1}'],
\)
which is the sample analogue of $\text{Var}(X_{i,-1})$. We assume that $V_x$ is
invertible, which is a minimal condition ensuring that the short regression is
well defined. Accordingly, imposing an upper bound on VCATE of the form
\(
\ensuremath{\mathbb{E}}_n\!\left[(\tau(x_i)-\ensuremath{\mathbb{E}}_n[\tau(x_i)])^2\right] \le C^2
\)
is equivalent to imposing the weighted quadratic constraint
\(
\delta'\,V_x\,\delta \le C^2
\)
on the parameter vector $\delta$. Throughout the paper, we report and plot results on the standard-deviation scale $C$ rather than the variance scale $C^2$, because $C$ shares the units of the outcome.
Although we do not pursue this direction here, the approach we introduce can be
extended straightforwardly to more general restrictions that require $\delta$ to
lie in a bounded, convex set. For example, one could replace the quadratic
constraint above with coordinatewise bounds of the form
$\max_{\ell}|\delta_\ell|\le C$ or other convex constraints on $\delta$ by following the general procedure described in
\citet{armstrong2023biasaware}.
\section{Heterogeneity-aware confidence intervals}
\label{sec:proposed-method}
In this section, we develop inference methods that explicitly account for
treatment effect heterogeneity. We begin by showing how to construct valid
confidence intervals based on estimators that are linear in outcomes, while
allowing for limited or even complete lack of overlap.
\subsection{Proposed method}\label{sec:bias-variance}
Building on~\eqref{eq:long reg}, we consider the regression model
\begin{equation}\label{eq: reg model fs}
Y_{i} = d_i\beta + x_{i}'\gamma + d_{i}\widetilde{x}_{i,-1}'\delta + \varepsilon_{i},
\end{equation}
where $Y_i$ is the outcome, $d_i$ is the (non-random) treatment indicator, and
$x_i=(1,x_{i,-1}')'$ is the vector of non-random covariates
including a constant. Define $D=(d_1,\ldots,d_n)'$ as the (realized) treatment vector and $X=(x_1,\ldots,x_n)'$
the matrix of (realized) covariates. We assume $X'X$ is invertible. To introduce our proposal
cleanly, we assume $\varepsilon_i \overset{i.i.d.}{\sim} N(0,\sigma^2)$ throughout this subsection, an
assumption we relax in Section~\ref{sec:implementation}. That is, treatment effect heterogeneity is fully captured by the covariates $x_i$
under this specification.
In a fixed-design setting, estimating the ATE under a VCATE bound is
equivalent to conducting inference on $\beta$ in~\eqref{eq: reg model fs}
subject to the constraint $\delta'V_x\delta \le C^2$ as explained in~\eqref{eq:wtd quadratic restriction}. Equivalently, the restriction on the parameter space is
$(\beta,\gamma',\delta')' \in \Theta_C$, where
\[
\Theta_C := \{(\beta,\gamma',\delta')' \in \mathbb{R}^{2+2k}:\ \delta'V_x\delta \le C^2\}.
\]
Let $Y=(Y_1,\ldots,Y_n)'$. We focus on linear estimators of the form
$\hat\beta_a = a'Y$, where $a\in\mathbb{R}^n$ is non-random. Since $(d_i,x_i')'$
are treated as fixed, the weights $a$ may depend on them; this class includes
typical estimators such as difference-in-means, the short regression, and the long
regression (when defined).
The variance of the estimator $\hat\beta_a=a'Y$ is
$\text{Var}(\hat\beta_a)=\sigma^2 a'a$ under homoskedasticity. From~\eqref{eq: reg model fs}, the worst-case
bias of $\hat\beta_a$ over $\Theta_C$ is
\begin{equation}\label{eq:worst-case bias}
\operatorname{\overline{bias}}_{a,C}
=
\sup_{(\beta,\gamma',\delta')\in\Theta_C}
\left\{a'(D\beta + X\gamma + (D\circ\tilde X_{-1})\delta) - \beta\right\},
\end{equation}
where $D\circ\tilde X_{-1}$ is the matrix with $i$th row equal to
$d_i\tilde x_{i,-1}'$. Hence, given the weights $a$, the worst-case bias of the
corresponding estimator can be calculated explicitly. Since $\Theta_C$ leaves $\beta$ and $\gamma$ unrestricted, it is necessary that the weights satisfy $a'D=1$ and $a'X=0$ to achieve finite worst-case bias. Under these conditions, the remaining bias $a'(D\circ\tilde X_{-1})\delta$ is linear in $\delta$, and since $\Theta_C$ is centrosymmetric in $\delta$ (i.e., $\delta\in\Theta_C$ implies $-\delta\in\Theta_C$), the one-sided supremum in~\eqref{eq:worst-case bias} equals $\sup_{\Theta_C}|\ensuremath{\mathbb{E}}[\hat\beta_a]-\beta|$.
A fixed-length confidence interval (FLCI) is constructed as
\begin{equation*}\label{eq:bias-corrected cv} \hat\beta_a \pm \chi_{a,C}, \qquad
\chi_{a,C} = \text{Var}(\hat\beta_a)^{1/2}\,
\text{cv}_\alpha\!\left(\operatorname{\overline{bias}}_{a,C}/\text{Var}(\hat\beta_a)^{1/2}\right),
\end{equation*} where $\text{cv}_\alpha(B)$ denotes the $(1-\alpha)$ quantile of the
folded normal distribution $|N(B,1)|$. The critical value therefore reflects the
potential worst-case bias of $\hat\beta_a$. Importantly, the validity of the
confidence interval does not rely on point identification of $\beta$. In cases
of complete lack of overlap, $\beta$ is set identified rather than point
identified under a VCATE bound; as shown in Lemma~\ref{thm:FLCI}, an immediate
consequence of \citet{armstrong_optimal_2018}, the FLCI attains correct coverage
uniformly over the identified set.
\begin{lemma}\label{thm:FLCI}
Under the assumption that VCATE is bounded by $C^2$ and the error terms in
\eqref{eq: reg model fs} are Gaussian,
\[
\inf_{(\beta,\gamma',\delta')\in\Theta_C}
\Pr\!\left(
\beta\in[\hat\beta_a-\chi_{a,C},\,\hat\beta_a+\chi_{a,C}]
\right)
\ge 1-\alpha.
\]
\end{lemma}
\begin{proof}
Under the stated assumptions,
$(\hat\beta_a-\beta)/\text{Var}(\hat\beta_a)^{1/2}\sim N(b,1)$, where the bias term
$b$ satisfies $|b|\le \operatorname{\overline{bias}}_{a,C}/\text{Var}(\hat\beta_a)^{1/2}$. The result follows by construction of
$\chi_{a,C}$.
\end{proof}
Write
$P_{X} = X(X'X)^{-1}X'$ and $H_{X} = I - P_{X}$. As shown in Lemma~\ref{claim:short bias} in the Appendix, the worst-case bias
of the short regression estimator is
\[
\operatorname{\overline{bias}}_{a_s,C}
=
C\sqrt{
a_s'(D\circ\tilde X_{-1})V_x^{-1}(D\circ\tilde X_{-1})'a_s},
\qquad
a_s=(D'H_XD)^{-1}H_XD.
\]
Because the short-regression weights $a_s$ do not adjust as $C$ increases,
bias-corrected short-regression confidence intervals can become excessively
long. Since the length of any FLCI of the form $\hat\beta_a\pm\chi_{a,C}$ equals
$2\chi_{a,C}$ and increases in both $\operatorname{\overline{bias}}_{a,C}$ and $\text{Var}(\hat\beta_a)$,
improving performance for $C>0$ requires a better bias-variance trade-off.
To achieve a better bias-variance trade-off, we propose solving the generalized
ridge least-squares problem
\begin{equation}
\label{eq:gen_ridge}
\min_{\beta, \gamma, \delta} \frac{1}{n} \norm{Y - D\beta - X\gamma - (D \circ
\tilde{X}_{-1})\delta}^{2} + \lambda\, \delta'V_{x}\delta
\end{equation}
where $\lambda>0$ is chosen as discussed later in Section~\ref{sec:lambda choice}. Let $\hat{\beta}_\lambda$ denote the resulting coefficient
estimator for $\beta$. Since $\delta'V_x\delta$ equals the sample VCATE as shown
in~\eqref{eq:wtd quadratic restriction}, the penalty directly shrinks the
heterogeneity coefficients $\delta$ toward zero, that is, toward the
constant-effects model, with the degree of shrinkage governed by $\lambda$. It is
not \textit{a priori} obvious that such a penalty on the outcome regression
should yield optimal inference for $\beta$, since penalized regression targets a
prediction error objective rather than the bias-variance trade-off specific to
$\beta$.\footnote{Under an $\ell_1$ bound on $\delta$, for instance, the
LASSO-penalized outcome regression does not produce optimal inference for
$\beta$. The coincidence between the ridge penalty on the outcome regression
and the optimal inference procedure for $\beta$ is specific to the quadratic
VCATE constraint, which has been pointed out in \cite{li1982minimaxity} and \cite{armstrong2023biasaware}.} However, the quadratic geometry of the VCATE restriction
implies that $\hat{\beta}_\lambda$ lies on the bias-variance frontier for every
$\lambda>0$: no linear estimator achieves simultaneously smaller variance and
smaller worst-case bias. Moreover, all linear estimators that lie on the frontier can be written in this form. Before we formalize this optimality result in
Theorem~\ref{thm:ridge frontier} in Section~\ref{sec:lambda choice}, we provide some intuition.
Lemma~\ref{lem:outcome ridge} shows that the generalized ridge regression coefficient estimator $\hat{\beta}_{\lambda}$ can be written as
\begin{equation}\label{eq:opt_sol}
\hat{\beta}_{\lambda} = a_{\lambda}'Y = \frac{\tilde{D}_{\lambda}'
Y}{\tilde{D}_{\lambda}'D},
\end{equation}
where the weights
$a_{\lambda} := \frac{\tilde{D}_{\lambda}}{\tilde{D}_{\lambda}'D}$ and the residuals
$\tilde{D}_{\lambda} := D - X\pi_{1,\lambda} - (D \circ
\tilde{X}_{-1})\pi_{2,\lambda}$ are from a penalized propensity score
regression:
\begin{equation}
\label{eq:prop_lag}
\min_{\pi_{1}, \pi_{2}} \frac{1}{n} \norm{D - X\pi_{1} - (D \circ
\tilde{X}_{-1})\pi_{2}}^{2} + \lambda \pi_{2}'V_{x}\pi_{2}.
\end{equation}
For $\lambda>0$ define the shrinkage matrix
\begin{equation*}
S_\lambda
:=
H_X (D\circ \tilde X_{-1})\,\bigl((D\circ \tilde X_{-1})'H_X (D\circ \tilde X_{-1})\\
+n\lambda V_x\bigr)^{-1} (D\circ \tilde X_{-1})'H_X,
\end{equation*}
which maps any vector $v \in\mathbb R^n$ to the fitted values from the
ridge regression of $H_X v$ on $H_X (D\circ \tilde X_{-1})$ with penalty matrix $n\lambda V_x$. With
this notation, the residualized treatment used by the generalized ridge
estimator can be written as
\begin{equation*}\label{eq:Dtilde_Slambda}
\tilde D_\lambda
=
H_XD - H_X(D\circ \tilde X_{-1})\pi_{2,\lambda}
=
(I-S_\lambda)H_XD,
\end{equation*}
and therefore
\[
\hat\beta_\lambda
=
\frac{\tilde D_\lambda'Y}{\tilde D_\lambda'D}
=
\frac{((I-S_\lambda)H_XD)'Y}{((I-S_\lambda)H_XD)'D}.
\]
Hence, the effect of $\lambda$ on $\hat\beta_\lambda$ operates entirely through
the shrinkage operator $S_\lambda$.
When $\lambda=0$ and $(D\circ \tilde X_{-1})'H_X(D\circ \tilde X_{-1})$ is invertible (i.e., the long regression is
well-defined),
\(
S_0
=
H_X(D\circ \tilde X_{-1})((D\circ \tilde X_{-1})'H_X (D\circ \tilde X_{-1}))^{-1}(D\circ \tilde X_{-1})'H_X
\)
is the orthogonal projection onto $\mathrm{span}(H_X (D\circ \tilde X_{-1}))$, and
$(I-S_0)H_XD$ equals the residual from regressing $D$ on $(X, D\circ \tilde X_{-1})$. Consequently,
$\hat\beta_0$ coincides with the long-regression coefficient on $D$.\footnote{When the long regression is not well-defined due to lack of overlap, in the case of discrete covariates, it is possible to analytically characterize the limit of $\hat\beta_\lambda$ as $\lambda\to0$, which is the long regression estimator $\hat{\beta}_{\text{long}}$ restricted to the trimmed subsample with overlap as shown in Lemma~\ref{lem:ridgeless}.} On the other hand, as $\lambda\to\infty$, it is clear that $\tilde D_\lambda\to H_XD$ and $\hat\beta_\lambda$ converges to the short
regression estimator.
To gain more intuition about our estimator, for fixed $\lambda>0$, suppose $X$ is generated from saturating discrete covariates. Lemma~\ref{thm:population ridge} shows the weights in~\eqref{eq:opt_sol} can
be written as $a_{\lambda,i}=\frac{1}{n}\frac{ (d_i-p(x_i))(1-\widetilde{x}_{i,-1}'\pi_{2, \lambda}) }{\ensuremath{\mathbb{E}}_n\left[ (d_i-p(x_i))(1-\widetilde{x}_{i,-1}'\pi_{2, \lambda}) d_{i}\right]}.$
The estimand of $\hat{\beta}_{\lambda}$ is therefore a weighted average of CATEs:
\begin{equation}
\beta_{\lambda}= \ensuremath{\mathbb{E}}[\hat{\beta}_{\lambda}] = \ensuremath{\mathbb{E}}_n\left[\frac{(p(x_i)(1-p(x_i)))(1-\widetilde{x}_{i,-1}'\pi_{2, \lambda})}{\ensuremath{\mathbb{E}}_n\left[(p(x_i)(1-p(x_i)))(1-\widetilde{x}_{i,-1}'\pi_{2, \lambda})\right]}\tau(x_{i})\right].\label{eq:ridge estimand}
\end{equation}
As shown in Lemma~\ref{lem:convex_weights},
expression~\eqref{eq:ridge estimand} further simplifies to
\begin{equation}\label{eq:explicit weights}
\beta_{\lambda} = \ensuremath{\mathbb{E}}_n\left[\frac{p(x_i)(1-p(x_i))/(p(x_i)(1-p(x_i))+\lambda)}
{\ensuremath{\mathbb{E}}_n\left[p(x_i)(1-p(x_i))/(p(x_i)(1-p(x_i))+\lambda)\right]}
\tau(x_{i})\right].
\end{equation}
The weight on each unit's CATE is proportional to
$p(x_i)(1-p(x_i))/(p(x_i)(1-p(x_i))+\lambda) \in [0,1]$, which equals
$1/(1+\lambda/[p(x_i)(1-p(x_i))])$. This factor is close to one for
well-overlapped cells where $p(x)(1-p(x))$ is large relative to $\lambda$,
and close to zero for poorly overlapped cells where $p(x)(1-p(x))$ is small
relative to $\lambda$. As
$\lambda\to 0$, the shrinkage factor approaches one for all cells and the
weights converge to cell shares, recovering the ATE. As $\lambda\to\infty$,
the weights converge to the short regression weights in~\eqref{eq:short estimand},
proportional to $p(x)(1-p(x))$. For intermediate $\lambda$, the ridge
estimand overweights well-overlapped cells relative to the ATE, but less so
than the short regression: cells with $p(x)(1-p(x)) \gg \lambda$ receive
weight close to their cell shares, while cells with
$p(x)(1-p(x)) \ll \lambda$ are heavily downweighted. Note that the weights
are always non-negative and sum to one. Later in Section \ref{sec:simulations}, we illustrate these weights in the context of
\cite{angrist1998} (see Figure \ref{fig:sim weights}).
\subsubsection{Choosing the penalty parameter}\label{sec:lambda choice}
It remains to choose the penalty parameter $\lambda$. We select it to minimize the half-length of the resulting FLCI in the homoskedastic Gaussian setting. By the same argument as Lemma~\ref{claim:short bias} in the Appendix, the worst-case bias of $\hat{\beta}_\lambda$ takes the form
\[
\operatorname{\overline{bias}}_{a_{\lambda},C}=C\cdot\sqrt{a_{\lambda}'(D\circ\tilde{X}_{-1})V_{x}^{-1}(D\circ\tilde{X}_{-1})'a_{\lambda}}.
\]
Let $\lambda^{\ast}_{C}$ denote the penalty parameter value that minimizes the half-length of the
corresponding fixed-length CI:
\begin{equation}\label{eq:opt_lam}
\lambda^{\ast}_{C} := \text{argmin}_{\lambda} \,\, \sigma\norm{a_{\lambda}} \cdot
\text{cv}_{\alpha}\left(\operatorname{\overline{bias}}_{a_{\lambda},C}/(\sigma \norm{a_{\lambda}})\right),
\end{equation}
and construct the corresponding CI $\hat{\beta}_{\lambda^*_{C}} \pm
\chi_{a_{\lambda^{\ast}_{C}},C}$. We refer to this inference procedure as \texttt{regulaTE}, for its ability to regularize treatment effect heterogeneity.
The \texttt{regulaTE} CI is in fact the shortest possible fixed-length CI based
on any linear estimator in the homoskedastic Gaussian setting. We formalize this
in the following theorem.
\newpage
\begin{theorem}\label{thm:ridge frontier}
The generalized ridge regression estimator $\hat{\beta}_{\lambda}$ solves the
bias-variance trade-off problem:
\begin{equation}
\label{eq:bias-var}
\min_{a \in \mathbb{R}^{n}} a'a \quad \text{s.t.} \quad \sup_{(\beta,
\gamma', \delta') \in \Theta_{C}} a'(D\beta + X\gamma + (D \circ
\tilde{X}_{-1})\delta) - \beta \leq B
\end{equation}
with $B = \operatorname{\overline{bias}}_{a_{\lambda}, C}$. Consequently, the
\textnormal{\texttt{regulaTE}} CI is the shortest fixed-length CI based on
linear estimators.
\end{theorem}
The first part of the theorem establishes that $\hat{\beta}_\lambda$ lies on the bias-variance frontier for every $\lambda > 0$: at bias level
$B = \operatorname{\overline{bias}}_{a_\lambda, C}$, no linear estimator achieves smaller variance.
This follows from the modulus of continuity characterization of optimal inference
\citep{donoho1994statistical,armstrong_optimal_2018}: the generalized ridge
regression solves the modulus problem associated with the VCATE constraint. The
second part follows because the family
$\{\hat{\beta}_\lambda : \lambda > 0\}$ spans the entire bias-variance frontier. Since
the FLCI half-length $\chi_{a,C}$ depends on the weights $a$ only through the
variance $a'a$ and worst-case bias $\operatorname{\overline{bias}}_{a,C}$, minimizing CI length over
all linear estimators reduces to minimizing over the ridge family, which is
what~\eqref{eq:opt_lam} achieves.
\subsection{Implementation and validity under more general error distributions}
\label{sec:implementation}
We begin with an initial estimator
$\hat{\theta}^{\mathrm{init}} = (\hat{\beta}^{\mathrm{init}},
\hat{\gamma}^{\mathrm{init}}{}', \hat{\delta}^{\mathrm{init}}{}')'$ for
$\theta = (\beta, \gamma', \delta')'$, obtained either from the long regression
(when defined) or from a cross-validated generalized ridge regression that
penalizes the weighted $\ell_{2}$ norm
$\ensuremath{\mathbb{E}}_n[(\widetilde{x}_{i,-1}'\delta)^2]$. Define the residuals from this initial
estimator
$\hat{\varepsilon}_{\mathrm{init},i} = Y_{i} -d_i\hat{\beta}^{\mathrm{init}} -
x_{i}'\hat{\gamma}^{\mathrm{init}}-
d_{i}\widetilde{x}_{i,-1}'\hat{\delta}^{\mathrm{init}}$, and the corresponding
initial variance estimator
$\hat{\sigma}^2 = \frac{1}{n} \sum_{i=1}^n \hat{\varepsilon}_{\mathrm{init},i}^2.$
For each $\lambda$, compute $\hat{\beta}_\lambda$ as before and obtain
$\lambda^*_{C}$ via \eqref{eq:opt_lam} by plugging in $\hat{\sigma}$. Form the
robust variance estimator
\[
\hat{V}_{\lambda^{\ast}_{C},\text{rob}} = \sum_{i=1}^n
a_{\lambda^{\ast}_{C},i}^2 \hat{\varepsilon}_{\mathrm{init},i}^2, \quad
\text{where } a_{\lambda^{\ast}_{C}} = \frac{\tilde{D}_{\lambda^{\ast}_{C}}}{\tilde{D}_{\lambda^{\ast}_{C}}'D}.
\]
The feasible CI is then $\hat{\beta}_{\lambda^*_{C}} \pm
\text{cv}_\alpha \left(
\operatorname{\overline{bias}}_{a_{\lambda^*_{C}},C}/\hat{V}_{\lambda^*_{C},\text{rob}}^{1/2} \right)
\cdot \hat{V}_{\lambda^*_{C},\text{rob}}^{1/2}$. Cluster-robust versions follow analogously.
\begin{remark}
Due to the heteroskedastic nature of the error terms, the exact optimality results
stated in Section~\ref{sec:proposed-method} no longer hold in general. Nonetheless,
the procedure mirrors the common practice of using weights that are optimal under
homoskedasticity (i.e., OLS weights) combined with robust standard errors.
\end{remark}
The asymptotic validity of the feasible CI is formally stated in
Appendix~\ref{sec:feasible CI primitive condition}, which follows from a
central limit theorem (CLT) applied to $\hat{\beta}_{\lambda^*_{C}}$. The key requirement is that
the maximal Lindeberg weight associated with the estimator,
\begin{equation*}
\operatorname{Lind}(a_{\lambda^*_{C}})\coloneqq \max_{1\le i\le n} \frac{a_{\lambda^*_{C},i}^2}{\sum_{j=1}^{n}a_{\lambda^*_{C},j}^2}
\end{equation*}
shrinks sufficiently quickly relative to the error of the initial estimator
used to form the residuals. While we provide formal conditions for asymptotic
validity in Appendix~\ref{sec:feasible CI primitive condition}, the Lindeberg
weights can be computed in any given application and serve as a diagnostic for
the reliability of the normal approximation; see \citet{noack2024bias} for an
analogous discussion in the context of fuzzy regression discontinuity designs.
The companion R package reports the maximal Lindeberg weight and issues a warning when it
is large. Note that, due to limited overlap, the convergence rate may be slower than the
parametric rate of $n^{-1/2}$.
\subsection{Necessity of bounding VCATE}
\label{sec:nec_bd_vcate}
The previous sections demonstrate how to construct \texttt{regulaTE} CIs given a
bound on VCATE. One might instead wish to construct a ``wider'' CI that remains valid under unrestricted treatment effect heterogeneity, or a CI that adapts to the true underlying VCATE, shrinking in length according to the true bound while still guaranteeing coverage over a broader class of heterogeneity. We show that neither approach is feasible, thereby establishing the necessity of imposing an \textit{a priori} bound on VCATE.
Intuitively, absent any restriction on the parameter space, the worst-case bias
of any linear estimator must be unbounded when overlap fails. The reason is that
the data contain no information about treatment effects for strata in which
only treated (or untreated) units are observed. Formally, recall
from~\eqref{eq:worst-case bias} that the bias is linear in the
parameters. Thus, only unbiased estimators have bounded bias (in fact, zero)
when the parameter space is unrestricted. As evident from the expression,
unbiasedness requires $a'D =1$, $a'X = 0$ and $a'(D \circ \tilde{X}_{-1}) =
0$. Suppose overlap fails for the binary $j$th covariate so that $X_{j+1} \leq D$,
where $\leq$ is interpreted elementwise. Writing
$\bar{x}_{j+1} := \ensuremath{\mathbb{E}}_{n}[x_{i, j+1}]\neq 0$, the conditions for unbiasedness
imply $a'D = 1$, $a'X_{j+1} = 0$ and $a'(X_{j+1} - D\bar{x}_{j+1}) = 0$, but it
is clear that no weight vector $a$ can satisfy these conditions.\footnote{This
is because due to lack of overlap $D \circ \tilde{X}_{-1} = X_{j+1} - D\bar{x}_{j+1}$.} Hence,
no unbiased linear estimator exists, and the worst-case
bias of \textit{any} linear estimator is unbounded. A formal statement and proof
of this result can be found in Appendix~\ref{appsec:imposs}.
Moreover, it turns out that it is impossible to construct a CI that adapts to the
true VCATE. By definition, an adaptive CI has length that automatically reflects the true
magnitude of VCATE while maintaining coverage under a conservative \textit{a priori}
bound on VCATE. However, Corollary 2.1 of
\cite{armstrong2023biasaware} implies that any valid CI must have expected length that
reflects the conservative \textit{a priori} bound $C^2$, even when VCATE is much smaller
than $C^2$. In other words, one cannot automate the choice of the VCATE
bound $C^2$ when constructing the CI.
Therefore, while $C$ is an important input to our method, it cannot be set to $C = \infty$, nor can it be selected in a data-driven way with the goal of adapting to the true VCATE.
In Section~\ref{sec: estimate VCATE}, we discuss how to conduct sensitivity analysis with respect to $C$.
\subsection{Calibrated simulations}\label{sec:simulations}
We illustrate the theoretical results so far in realistic
settings through simulations calibrated to the data generating process in \cite{angrist1998}.
\cite{angrist1998} links social security earnings records to
administrative data on a sample of U.S. military applicants from 1979 to 1982 to
estimate the effects of voluntary military service on veterans' earnings. The following discrete characteristics are controlled for confounding: year of application, year
of birth, education at the time of application, race, and Armed Forces
Qualification Test (AFQT) score group. The paper documents heterogeneity in treatment effects: military service modestly
increased the civilian earnings of non-white veterans while reducing those of
white veterans. Further heterogeneity is observed across background
characteristics such as education and AFQT scores, which prompted
\cite{angrist1998} to theoretically analyze the different estimands of the short
regression and the long regression (referred to therein as the controlled contrast
estimator). The estimates were found to be significantly different, both
statistically and economically.
For simplicity, our simulation exercises focus on inference for the average treatment effect on earnings in 1988 among white applicants. Within this population, there are approximately 400 covariate cells. The public replication data from \cite{angrist1998} report cell-level summary statistics, such as the mean and standard deviation of earnings, the share of veterans, and the cell frequency, constructed from administrative records covering roughly 100,000 individuals.
To calibrate the simulation to \cite{angrist1998} and to construct a micro-level dataset, we draw 2,000 individuals, which can be interpreted as a 2\% subsample of the original population. Treatment status is assigned using the cell-level share of veterans as the true propensity score, which ranges from 4.6\% to 81.2\%. Given a fixed realization of treatment assignments, the relatively small sample size leads to a lack of overlap at some covariate values, equivalent to 12.4\% of the sample size. As a result, the long regression is not well defined and the ATE is not point identified.
To preserve both heteroskedasticity and treatment effect heterogeneity from the original data of \cite{angrist1998} in the data-generating process, outcomes are generated by treating the cell-level summary statistics as the true means and standard deviations of earnings and assuming normally distributed earnings within each cell. We treat the cell-level standard deviations as known for the simulation. The true standard deviation of the CATEs is $C_0=1452.195$ dollars.
\begin{figure}[htbp!]
\begin{center}
\caption{\label{fig:sim weights}Weighted Average Interpretation of \texttt{regulaTE} Estimators under Lack of Overlap in DGP Calibrated to Angrist (1998)}
\includegraphics[width=0.7\linewidth]{fixed_design_heterosked_no_overlap_weights_notrim.pdf}
\end{center}
\end{figure}
Because all covariates are discrete in the data-generating process, the \texttt{regulaTE} estimand admits a transparent interpretation as a weighted average of the cell-level treatment effects $\tau(x)$, as characterized in~\eqref{eq:explicit weights}. Figure~\ref{fig:sim weights} plots these cell-level \texttt{regulaTE} weights, after excluding cell shares, under several values of the heterogeneity bound $C$ against the within-cell treatment variance $p(x)(1-p(x)).$ When $C=0$, to minimize variance \texttt{regulaTE} coincides with the short regression, placing weights proportional to $p(x)(1-p(x)).$ As $C$ increases, \texttt{regulaTE} moves away from the short regression, but still places relatively high weight on well-overlapped cells and aggressively shrinking contributions from poorly overlapped ones where $p(x)(1-p(x))$ is small. With $C=C_0$ where treatment effect heterogeneity can be large, the \texttt{regulaTE} weights become closer to cell shares due to the increasing importance of worst-case bias relative to variance in the bias-variance trade-off.
\begin{figure}[htbp!]
\begin{center}
\caption{\label{fig:sim no overlap}Sensitivity of Coverage and CI Length under Lack of Overlap in DGP Calibrated to Angrist (1998)}
\begin{subfigure}{0.45\linewidth}
\includegraphics[width=\linewidth]{fixed_design_heterosked_no_overlap_CI_ratio_notrim.pdf}
\end{subfigure}
\hfill
\begin{subfigure}{0.45\linewidth}
\includegraphics[width=\linewidth]{fixed_design_heterosked_no_overlap_coverage_notrim.pdf}
\end{subfigure}
\end{center}
\footnotesize{Notes: ``\texttt{regulaTE} CI'' refers to the bias-aware fixed-length confidence interval based on the \texttt{regulaTE} estimate. ``Short CI'' refers to the CI based on the short regression estimate. ``Short BC CI'' refers to the bias-corrected short regression CI. Both ``\texttt{regulaTE} CI'' and ``Short BC CI'' are heteroskedasticity-robust with 95\% coverage guarantees under each heterogeneity bound $C$ on the horizontal axis. The ratio of the CI lengths is relative to the length of the \texttt{regulaTE} CI under $C=C_0$.}
\end{figure}
On the horizontal axis of Figure~\ref{fig:sim no overlap}, we consider various heterogeneity bounds $C$ and use them both to bias-correct the short-regression CI and to construct the heteroskedasticity-robust \texttt{regulaTE} CI. When no correction is applied, the short-regression CI, while quite short, exhibits substantial undercoverage, reflecting bias induced by omitting treatment effect heterogeneity. As the bound $C$ increases, both the bias-corrected short-regression CI and the \texttt{regulaTE} CI achieve coverage close to the nominal level. But across the range of $C$, in this heteroskedastic setting, the bias-corrected short-regression CI is longer than the \texttt{regulaTE} CI, as shown in the left panel.
Correct coverage is attained for values of $C$ that are strictly smaller than the true heterogeneity level $C_0$ because the validity guarantee is derived under worst-case heterogeneity, whereas the data-generating process considered here is less adversarial.
In Appendix~\ref{sec:appendix simulations}, we illustrate the behavior of \texttt{regulaTE} in settings with overlap and compare it with the long regression, which is now well-defined. We also compare with the adaptive estimator of \cite{armstrong2023adapting}, which combines the short and long regression estimators for efficiency gain. \texttt{regulaTE} is constructed to adjust the CI based on the user-specified bound $C$ for sensitivity analysis. Therefore, it connects the long and short regression CIs as $C$ increases from zero. The adaptive estimator does not adjust to the user-specified bound $C$ and remains close to the long regression in this DGP, while the bias-corrected short-regression CI is again overly long.
\section{Sensitivity analysis and empirical illustrations} \label{sec: estimate VCATE}
Researchers often use the short regression to estimate the ATE under the
(frequently implicit) assumption that treatment effect heterogeneity is not too
large. Our method provides a formal way to assess this assumption. In this case,
rather than taking a definitive stand on $C$, one can begin with $C = 0$ (which
corresponds to using the short regression) and gradually increase $C$ until the
results become insignificant.\footnote{Note that we take $C$, rather than $C^2$, as our sensitivity parameter because it has the same units as the outcome.} We denote the smallest value of $C$ that renders the estimated treatment effect insignificant as $C^\ast$. One can then evaluate whether the corresponding
breakdown point $C^\ast$ is plausible, which is feasible since VCATE is a highly
interpretable quantity. This is analogous to the breakdown frontier analysis
considered in, e.g., \cite{kline2013sensitivity}, \cite{masten2020inference} and \cite{li2021linear}. Our R
package provides a plotting functionality that aids such sensitivity analysis.
\subsection{Unconfoundedness}
\cite{aizer2016} evaluate the long-run effects of early 20th-century Mothers' Pension (MP) cash transfers on children's lifetime outcomes, using administrative data. Their study compares children of approved MP applicants to those whose approvals were subsequently reversed, accounting for observed characteristics, to estimate causal effects. The original study
estimates the following short regression as in \citet[Equation (1)]{aizer2016}:
\[
\log(\textit{Age at Death})_{ifts} = \theta_0 + \theta_1 MP_f + \theta_2 \boldsymbol{X}_{if} + \theta_3\boldsymbol{Z}_{st} + \boldsymbol{\theta}_c + \boldsymbol{\theta}_t + \varepsilon_{if},
\]
where the outcome is the natural logarithm of the age at death for a given individual $i$ in family $f$ born in year $t$ living in a county $c$ (state $s$), and the treatment $MP_f$ indicates MP receipt. The authors control linearly for $\boldsymbol{X}$, a vector of relevant family characteristics (marital status, number of siblings, etc.) and child characteristics (year of birth and age at application), and $\boldsymbol{Z}_{st}$, a vector of county-level characteristics in 1910 and state characteristics in the year of application. The authors also control for county and cohort fixed effects ($\boldsymbol{\theta}_{c}$ and $\boldsymbol{\theta}_{t}$). To illustrate our method, we focus on child longevity, which is one of the authors' main outcome of interest.
We assess the sensitivity of the estimate for ATT, which is the average impact of MP receipt among MP recipients, as it is the most relevant target parameter for program evaluation.
Based on the short regression, the ATT is estimated to be $1.82$\% and is statistically significant at the $5$\% significance level. Note that the long regression is infeasible because in eleven counties, accounting for about 9\% of the sample, all applicants received the MP. This lack of overlap renders the long regression undefined due to collinearity among the MP receipt indicator, interaction terms and county fixed effects.
After reporting their ATT estimates based on the short regression~\cite[Table 4 Panel A]{aizer2016}, the authors examine heterogeneous effects of the MP program across subgroups defined by family income, the child's age, and urban residence. Some subgroup estimates are statistically insignificant, and the magnitudes are broadly similar, ranging between 1\% and 2\% \citep[Table 5]{aizer2016}.
\begin{figure}[htbp!]
\begin{center}
\caption{\label{fig:specification_3}Sensitivity Analysis for \cite{aizer2016}}
\includegraphics[width=0.8\linewidth]{log_age_at_death_vs_het_bound_spec_3_att_package_breakdown.png}
\label{fig:specification_3_att}
\end{center}
\footnotesize{Notes: ``\texttt{regulaTE} Point Estimate'' refers to the \texttt{regulaTE} estimate optimized under each heterogeneity bound $C$ on the horizontal axis. ``Short Point Estimate'' and ``Short BC Point Estimate'' both refer to the short regression estimate. ``\texttt{regulaTE} CI'' and ``Short BC CI'' refer to the bias-aware fixed length confidence interval based on the corresponding estimates, with 95\% coverage guarantees valid under each heterogeneity bound $C$.}
\end{figure}
To formalize the robustness checks based on the subgroup analysis originally conducted by the authors, Figure~\ref{fig:specification_3} conducts the sensitivity analysis based on the \texttt{regulaTE} CI. We also present the bias-corrected short regression CI for comparison. Standard errors are clustered at the
county-level, as in the original analysis. As $C$ increases, the \texttt{regulaTE} CI for ATT expands to account for the possibility of greater worst-case bias. Note that the bias-corrected short regression CI is substantially wider, underscoring the efficiency gains from \texttt{regulaTE}. The breakdown point is around 0.75\% ($C^\ast=0.0075$). To put this breakdown frontier value into perspective, note that given an average age at death of 72.44 years, even a 2\% increase represents a substantively meaningful effect. This consideration motivates a benchmark on the standard deviation of heterogeneous percentage effects of $C^{\text{ref}} = 1\%$, corresponding to bounding the effects to be between 0\% and 2\% and applying Popoviciu's inequality on variances. The breakdown point is only marginally smaller than $C^{\text{ref}}$, suggesting the statistically significantly positive impact of MP receipt on child longevity among receiving families is robust to economically meaningful treatment effect heterogeneity to some extent.
\subsection{Staggered adoption}\label{sec:application}
The recent literature on staggered adoption designs also emphasizes the
bias in common specifications from omitting treatment effect heterogeneity; see
\cite{roth2023s} for a review. Consider a staggered
adoption setting where we observe a sample of i.i.d. draws
$\left\{ \{Y_{it}, d_{it}\}_{t=0}^{T}\right\}_{i=1}^n$ over $\mathcal T =T+1$
total time periods. Define $e_i := \min\{ t : d_{it} = 1 \}$ as the period in
which unit $i$ first receives treatment. Let $\mathcal{E}$ denote the set of all
such first treatment periods, and define a cohort as the set of units treated
for the first time in the same period: $\{ i : e_i = e \}$.
Under standard assumptions of parallel trends and no anticipation, the DGP for the observed outcome can be written as
\begin{equation}
Y_{it} = \alpha_e + \lambda_t + \sum_{\ell \geq 0, \, e \in \mathcal{E}} \tau_{e,\ell} \cdot d_{it} \cdot \mathbf{1}\{ e_i = e , t - e_i = \ell \} + \varepsilon_{it},
\label{eq:saturated TWFE}
\end{equation} where $\alpha_e$ and $\lambda_t$ are cohort and time fixed
effects respectively. The framework developed earlier for unconfoundedness settings in~\eqref{eq:dgp_selection} assumes $\tau(\cdot)$ is a linear function of the confounders $x_i$. Here $\tau(\cdot)$ is a linear function of group and relative time indicators, rather than the confounders. Specifically,
$\tau_{e,\ell} = \ensuremath{\mathbb{E}}[ Y_{i,e+\ell}(1) - Y_{i,e+\ell}(0) \mid e_i = e ]$
is the conditional average treatment effect (CATE) for cohort $e$ at relative time $\ell$, which can vary across cohorts and over time in an unrestricted
fashion. Therefore, the original framework continues to hold for the following ``short'' and ``long'' regressions.
Rather than the unrestricted~\eqref{eq:saturated TWFE}, researchers often estimate a simpler specification also known as the ``static'' two-way fixed effects (TWFE) regression
\begin{equation}
Y_{it} = \alpha_e + \lambda_t + \beta_{\text{short}} \cdot d_{it} + \varepsilon_{it}
\label{eq:static TWFE}
\end{equation}
which omits the interactions between $d_{it}$ and cohort and relative-time
indicators. Analogous to the short regression in the unconfoundedness setting,
as argued in \cite{double_fe} and \cite{goodman2021difference}, the ``static''
specification~\eqref{eq:static TWFE} implicitly assumes constant effects across
cohorts and over time (i.e., $C = 0$). When the effects indeed vary, the ``static'' specification~\eqref{eq:static TWFE} can be severely biased for a class of reasonable estimands due to negative weighting of $\tau_{e,\ell}$ for large $\ell$.
One reasonable estimand in this setting is the average effect over treated
cohorts (ATT), defined as
\[
\text{ATT} = \frac{
\sum_{t \in \mathcal{T}} \sum_{e \in \mathcal{E}} \mathbb{P}_n\{d_{it}=1\}\cdot\mathbb{P}_n\{e_i = e \mid d_{it} = 1 \} \cdot \tau_{e,t-e}
}
{ \sum_{t \in \mathcal{T}} \mathbb{P}_n\{d_{it}=1\},
}
\]
where $\mathbb{P}_n$ denotes empirical frequencies. That is, the ATT is a weighted average of the cohort-by-event-time effects $\tau_{e,t-e}$, where the weights reflect the empirical distribution of treated observations across cohort-time cells. To recover this estimand, we can reparameterize~\eqref{eq:saturated TWFE} as the ``long'' regression
\[
Y_{it} = \alpha_e + \lambda_t + d_{it} \beta^{\text{ATT}}_{\text{long}} + \left( d_{it} \cdot \tilde{x}_{it} \right)'\delta + \varepsilon_{it}
\]
where $\tilde{x}_{it}$ collects all but one of the re-centered cohort and relative time indicators:
\[
\tilde{x}_{it} =\begin{pmatrix}
\vdots \\[4pt]
\mathbf{1}\{ e_i = e,\, t - e_i = \ell \}
-
\dfrac{
\sum_{t \in \mathcal{T}} \mathbb{P}_n\{d_{it}=1\}
\mathbb{P}_n\{e_i = e \mid d_{it} = 1\}
\mathbf{1}\{t-e=\ell\}
}{
\sum_{t \in \mathcal{T}} \mathbb{P}_n\{d_{it}=1\}
} \\[6pt]
\vdots
\end{pmatrix}.
\]
This long regression coincides with the extended TWFE specification from \cite{wooldridge2025twfe}. However, this long regression can be noisy or even infeasible in practice. Estimation of~\eqref{eq:saturated TWFE} corresponds to aggregating
cohort- and time-specific DID estimators $\hat{\tau}_{e,t-e}$ for each
$\tau_{e,\ell}$, which can be noisy when few units are treated at a given time. Moreover, if all units are treated after time $t$, the required counterfactuals are not identified, similar to a lack of overlap in cross-sectional settings, and the ATT is no longer identified without further assumptions on treatment effect heterogeneity. To address this, we impose a bound $C > 0$ on the variance of treatment effect heterogeneity, restricting the deviation of $\tau_{e,t-e}$ from the ATT. Under this assumption, we can construct valid confidence intervals using our \texttt{regulaTE} procedure, thereby extending the sensitivity analysis to the staggered adoption setting.
To illustrate, we revisit the empirical application in \cite{sun_estimating_2020}, which builds on \cite{dobkin_finkelstein_kluender_notowidigdo_aer2018} to estimate the average effect of an unexpected hospitalization on out-of-pocket medical spending among adults who experienced at least one such hospitalization between waves 8 and 11 of the Health and Retirement Study (HRS). Each wave corresponds to approximately two calendar years. Using a balanced panel from wave 7 through wave 11 ($T=5$), all individuals are treated in the final period ($t=5$), so the ATT from waves 8 to 11 is not point identified and the long regression~\eqref{eq:saturated TWFE} is not well defined. Nevertheless, the ``static'' specification delivers a statistically significant and positive estimate, as shown in Figure~\ref{fig:oop_spend}.
\begin{figure}[!htb]
\begin{center}
\caption{\label{fig:oop_spend}Sensitivity Analysis for \cite{dobkin_finkelstein_kluender_notowidigdo_aer2018}}
\includegraphics[width=0.8\linewidth]{oop_spend_vs_het_bound_spec_static_breakdown.png}
\end{center}
\footnotesize{Notes: ``\texttt{regulaTE} Point Estimate'' refers to the \texttt{regulaTE} estimate optimized under each heterogeneity bound $C$ on the horizontal axis. ``Short Point Estimate'' and ``Short BC Point Estimate'' both refer to the short/static regression~\eqref{eq:static TWFE} estimate. ``\texttt{regulaTE} CI'' and ``Short BC CI'' refer to the bias-aware fixed length confidence interval based on the corresponding estimates, with 95\% coverage guarantees valid under each heterogeneity bound $C$.}
\end{figure}
To evaluate the robustness of this estimate to treatment effect heterogeneity, Figure~\ref{fig:oop_spend} reports sensitivity analysis based on the \texttt{regulaTE} CIs. Again, we present the bias-corrected short regression CIs as well for comparison. Standard errors are clustered at the individual level, as in the original analysis of \cite{dobkin_finkelstein_kluender_notowidigdo_aer2018}. The CIs based on \texttt{regulaTE} remain significant up to a breakdown value of $C^\ast = \$1{,}583$. To interpret this breakdown value, consider the original analysis in \cite{dobkin_finkelstein_kluender_notowidigdo_aer2018}, which focuses on heterogeneity over time. They report an increase in out-of-pocket spending of roughly \$3{,}000 in the first year after hospitalization and \$1{,}000 by the third year. Under the assumption that heterogeneity operates primarily over time, treating these magnitudes as rough bounds and applying Popoviciu's inequality on variances, we get a reference treatment effect heterogeneity value of $C^{\text{ref}} = \$1{,}000$. Therefore, using sensitivity analysis based on the \texttt{regulaTE} CIs, we find the statistically significant average increase in out-of-pocket spending due to unexpected hospitalization remains robust to substantial treatment effect heterogeneity. In contrast, the bias-corrected short regression CIs widen substantially, and a sensitivity analysis based on them would suggest a lack of robustness to treatment effect heterogeneity, potentially overstating the fragility of the original results.
\section{Conclusion}\label{sec:conclusion}
Many specifications commonly used in empirical research implicitly restrict
treatment effects to be constant, often favoring precision at the expense of
robustness. While such a restriction yields narrower confidence intervals, the resulting
estimates may be substantially biased in the presence of heterogeneity. This paper develops a
sensitivity analysis based on the proposed \texttt{regulaTE} CIs by varying the
bound on treatment effect heterogeneity.
There are several directions for future work. One natural extension is to
consider alternative forms of restrictions on treatment effect
heterogeneity. While variance bounds naturally capture the dispersion in
heterogeneous effects, one could instead impose an absolute bound on the
deviation of individual treatment effects from the average effect if such prior
knowledge is available. Another promising direction is to generalize the
framework to accommodate multiple treatments or layered sources of
heterogeneity, as in \cite{goldsmithpinkham2024contamination} and
\cite{sun_estimating_2020}.
\bibliographystyle{chicago}
\bibliography{ref}
\newpage