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.
99,463 characters · 17 sections · 36 citation commands
A Machine-Learning-Compatible Omnibus Test for Treatment Effect Heterogeneity
\def\spacingset#1{ {#1}} \spacingset{1}
\if11 \fi
\if01 {
} \fi
\thispagestyle{empty}
{\it Keywords:} Average Treatment Effects, Heterogeneity, Machine Learning
{\it JEL classification:} C12, C14, C21, C23, C26
\spacingset{1.1}
Understanding whether treatment effects vary is central to policy evaluation because average effects may conceal differences in how policies, programs, or interventions affect different groups. Heterogeneous effects are relevant in many fields, including economics, political science, and health sciences, where researchers routinely investigate whether an intervention benefits some groups more than others. For example, one may want to assess whether a job-training program improves employment differently across skill levels, whether the effectiveness of a medical treatment varies by the age of patients, or whether an educational intervention influences learning gains differently across socioeconomic groups. In all these cases, a valid hypothesis test for treatment effect heterogeneity is essential.
In recent years, machine learning (ML) methods have become increasingly popular for uncovering such heterogeneity \citetext{see, e.g., athey2019machine and chern2024}. Compared with traditional econometric approaches such as interacted models or subgroup analysis, ML-based techniques typically require fewer parametric assumptions, offer more flexible modeling with safeguards against overfitting, and can accommodate large sets of control variables, often high-dimensional ones. Many of these methods were originally developed for prediction problems. In many empirical applications, however, researchers are interested in causal parameters that are directly relevant for policy evaluation, such as average treatment effects (ATEs) and their variation across groups. Traditional econometric methods therefore place particular emphasis on identification and valid statistical inference for such quantities. Recent advances in Causal-ML combine the flexibility of ML with the inferential objectives of econometrics by using ML methods to estimate high-dimensional nuisance functions while preserving valid inference for causal effects. Despite this progress, there remains a need for formal procedures that can assess treatment-effect heterogeneity while retaining the flexibility of modern ML methods and that can be readily implemented in practice.
This study contributes to the small but growing literature on detecting systematic treatment effect heterogeneity with ML (see literature review below). It develops a nonparametric omnibus test for treatment-effect heterogeneity that can be combined in a modular way with a wide range of estimators, including modern ML methods. The test provides a single global assessment of whether treatment effects exhibit systematic variation with respect to a chosen low-dimensional set of covariates, while allowing identification to rely on high-dimensional controls. We construct a locally robust test statistic that remains valid under flexible nuisance estimation and pair it with a computationally efficient bootstrap. The bootstrap requires nuisance functions to be estimated only once, avoiding repeated first-stage estimation in each iteration. Together, these features yield a practical tool for large empirical applications and for settings where heterogeneity is of interest but conventional low-dimensional methods are inadequate.
The broad applicability of the procedure is illustrated across several common empirical designs. We first consider heterogeneity in ATEs under conditional unconfoundedness, the standard identifying assumption in observational studies. The same framework naturally covers randomized experiments as a special case in which treatment is independent of covariates and nuisance functions are either known or estimated only to improve precision (for example, list2025, shows that ML-based covariate adjustment in experiments can enhance efficiency without affecting identification). We then extend the approach to difference-in-differences (DiD) designs to test heterogeneity in average treatment effects on the treated (ATT) under conditional parallel trends, accommodating both panel data and repeated cross-sections. Finally, we adapt the method to instrumental-variables settings to test for heterogeneity in local average treatment effects (LATE) when instruments are valid conditional on covariates.
For each design, we derive the corresponding test statistic and establish its asymptotic null distribution. The construction follows Integrated Conditional Moment principles in the spirit of bierens2016econometric. We reformulate the conditional null of no treatment-effect heterogeneity as a continuum of unconditional moment restrictions and integrate these moments to obtain a feasible test. Because the limiting distribution is not directly useful for implementation, we introduce a bootstrap approximation that is asymptotically valid under the null. The procedure requires estimating nuisance components only once, which we do using ML methods with cross-fitting. Neyman orthogonality of the statistic ensures local robustness to regularization bias, so small first-stage errors do not affect the asymptotic distribution. Estimating nuisances a single time and relying on simple matrix operations within the bootstrap makes the method computationally efficient.
Simulation results confirm the practical performance of the test. Across a range of sample sizes, covariate dimensions, and ML estimators for the nuisance components, the procedure attains close-to-correct size under the null and exhibits high power against relevant heterogeneous alternatives.
We illustrate the empirical relevance of the procedure using two applications. The first revisits the effect of eligibility for a 401(k) retirement savings plan on household financial assets under conditional unconfoundedness, following poterba_401k_1995 and chernozhukov2018double. The second considers treatment effect heterogeneity in a difference-in-differences setting studying the impact of tariff reductions on corruption along the South Africa--Mozambique trade corridor, building on sequeira_corruption_2016 and chang2020double. Across both applications, the proposed test detects economically meaningful heterogeneity while remaining computationally feasible in settings with flexible ML-based nuisance estimation.
Related Literature. This study contributes to a growing body of work on testing and estimating heterogeneous treatment effects using semi- and nonparametric tools, often combined with modern ML. We situate our approach within three related strands of the literature: a smaller literature on global inference for systematic heterogeneity, a smaller literature on tests of identical unit-level effects, and a large literature that uses ML under unconfoundedness to estimate heterogeneous treatment effects through conditional average treatment effects (CATEs) in high-dimensional settings.
A first strand develops global procedures for testing whether treatment effects vary systematically with observed covariates. Early work such as crump2008nonparametric proposes nonparametric homogeneity tests across covariate-defined subpopulations, but these methods are mainly designed for low-dimensional covariates and classical estimators. More recent contributions combine global heterogeneity assessment with ML methods. chernozhukov2023fisherschultzlecturegenericmachine develop a general sample-splitting framework for inference on low-dimensional summaries of treatment-effect variation, including best-linear predictors and sorted effects. In this framework, heterogeneity is assessed through an ML-generated score, either by testing whether the score is linearly related to treatment-effect variation or by comparing average effects across groups ranked by the score. Building on related ideas, Imai2025a propose a nonparametric heterogeneity test for randomized experiments based on sorted effects. Our approach differs from these contributions in three main respects. First, our objective is to test directly whether treatment effects vary systematically with a chosen set of covariates, without requiring heterogeneity to be summarized by a linear projection or a small number of ranked groups. Second, while these approaches focus primarily on randomized experiments or settings under unconfoundedness, our framework extends to several empirical designs, including difference-in-differences and instrumental-variables settings. Third, our procedure remains computationally efficient because the bootstrap avoids repeated re-estimation of nuisance components.
A second strand of research studies randomization-based tests built for the sharp hypothesis of identical unit-level effects. Key contributions include heck97 and ding2016randomization,ding2019decomposing. These tests are designed to detect any deviation from constant individual effects, including purely idiosyncratic variation. In contrast, our target is a weaker but empirically relevant null, namely the absence of explainable heterogeneity with respect to observed covariates. This distinction matters in applications where the goal is to determine whether effects differ across observable groups in a learnable and actionable way, rather than whether unit-level heterogeneity is driven by potentially idiosyncratic and unpredictable variation.
A third strand examines effect heterogeneity under unconfoundedness using flexible ML methods, focusing primarily on the estimation and inference of CATEs and related functionals in high-dimensional settings. This literature includes tree-based methods such as causal trees and causal forests ath2016,wager2018estimation,ath2019,Fried21,lech2022, debiased, orthogonalized, and partialling-out approaches for high-dimensional CATE inference semenova2021debiased,chernozhukov2018,chernozhukov2021,fan2022,lech2019,semenova2017,nekipelov2022,Foster2019,kennedy2023towards,Nie2020,scheidegger2026inference, approaches imposing additional structure on effect modifiers, such as sparsity in targeted undersmoothing Hansen2017, and Bayesian models for CATE heterogeneity based on hierarchical or prior-driven structure green2012,tad2016,hahn2020. Taken together, this strand provides powerful tools for estimating and interpreting heterogeneous effects, primarily under unconfoundedness, but it does not provide a general omnibus test of whether systematic heterogeneity is present with respect to observed covariates.
Finally, our paper is related to escanciano2026kernel, who develops a specification test for moment restrictions using Reproducing Kernel Hilbert Spaces, in settings where the null hypothesis is characterized in terms of high-dimensional nuisance parameters.\footnote{This paper was presented under the title "Kernel-Based Specification Testing with High-Dimensional Nuisance Parameters" at the ISNPS Conference in Thessaloniki, June 2026. A written version of the paper was not yet available at the time of writing.} Our contribution differs from his along three dimensions. First, we focus on a specific class of moment restrictions arising from treatment-effect heterogeneity, rather than on conditional moment restrictions more broadly. Second, we establish an expansion of the empirical process underlying our test statistic that is uniform with respect to the indexing parameter $t$, which underlies our asymptotic results. Third, whereas escanciano2026kernel obtains critical values through a multiplier bootstrap based on the influence-function representation of the test statistic, we propose a weighted bootstrap that departs from the multiplier scheme; in our empirical applications, this construction avoids the computation of certain terms that would otherwise complicate implementation.
Outline. Section (ref) provides an intuitive overview of the proposed testing framework and its motivation. Section (ref) develops the test for heterogeneity in ATEs under unconfoundedness, including the locally robust statistic and bootstrap procedure. Section (ref) extends the framework beyond unconfoundedness to difference-in-differences designs with panel and repeated cross-section data, as well as to LATEs identified through instrumental variables. Section (ref) reports simulation evidence on size and power. Section (ref) presents two empirical applications illustrating the proposed procedure in observational and difference-in-differences settings. The final section concludes.
Researchers are often interested not only in whether an intervention affects outcomes on average, but also in whether its effects differ across groups. The potential-outcomes framework provides a natural way to formalize this question. Let $D_i\in\{0,1\}$ denote a binary treatment indicator, where $D_i=1$ if unit $i$ receives the treatment and $D_i=0$ otherwise. Depending on the application, the treatment may represent participation in a job-training program, eligibility for a tax-advantaged savings plan, access to a microcredit program, or exposure to a policy reform. Let $Y_i(d)$ denote the potential outcome under treatment status $d\in\{0,1\}$, such as employment, household income, savings, or firm performance. Thus, $Y_i(1)$ is the outcome that would be observed if unit $i$ received the treatment, whereas $Y_i(0)$ is the outcome that would be observed if the same unit did not receive the treatment. The difference $Y_i(1)-Y_i(0)$ is the individual treatment effect. Let $X_i$ denote a vector of observed pre-treatment covariates, such as age, education, earnings and employment histories, household characteristics, or pre-treatment measures of business performance. For notational convenience, let $(Y(1),Y(0),D,X)$ denote generic population random variables corresponding to the unit-level quantities $(Y_i(1),Y_i(0),D_i,X_i)$.
A fundamental challenge of causal inference is that only one of the two potential outcomes is observed for any given unit. Individual treatment effects are therefore not directly observed. Nevertheless, averages of these individual effects can be identified from observed data under appropriate assumptions. A central example is the ATE, $\delta_0=\operatorname{E}\{Y(1)-Y(0)\}$, which measures the average impact of the intervention in the population. For example, in an evaluation of a job-training program, $\delta_0$ captures the average change in employment that would occur if everyone in the population participated in the program rather than if no one did.
In many empirical applications, researchers face two related but distinct tasks. The first is identification: estimating causal effects may require adjusting for a rich set of covariates. The second is interpretation of treatment-effect heterogeneity: substantive interest is often focused on whether effects differ across a smaller set of policy-relevant characteristics. We therefore distinguish between the full covariate vector $X$, which is used for identification, and the lower-dimensional subvector $X_c\subseteq X$, which contains the covariates with respect to which heterogeneity is assessed.
This distinction arises naturally in practice. In evaluations of job-training programs, identification may require controlling for detailed labor-market histories, demographic characteristics, regional indicators, and prior earnings, while the main question may be whether effects differ by age or educational attainment. In studies of medical treatments, adjustment may require many clinical and behavioral variables, while interest may center on whether treatment effects vary with baseline health status. In educational applications, researchers may control for a rich set of student, family, and school characteristics, while examining whether an intervention has different effects by gender, migration background, or socioeconomic status. In each case, the full vector $X$ may be high-dimensional, whereas the heterogeneity question concerns a lower-dimensional set of covariates $X_c$.
Treatment effects may also be summarized conditionally on observed characteristics. We write \[ \delta(X):=\operatorname{E}\{Y(1)-Y(0)\mid X\} \] for the CATE conditional on the full covariate vector $X$. Since our test focuses on heterogeneity with respect to the lower-dimensional covariate vector $X_c$, we also define \[ \delta_c(X_c):=\operatorname{E}\{Y(1)-Y(0)\mid X_c\}. \] Thus, $\delta(X)$ describes treatment-effect variation with respect to the full set of covariates, whereas $\delta_c(X_c)$ describes treatment-effect variation with respect to the covariates used in the heterogeneity test.
Our objective is to test whether treatment effects vary systematically with $X_c$. Formally, for the ATE, we consider the null hypothesis
where $\mathcal{H}_0^c$ denotes the logical complement of $\mathcal{H}_0$. Under the null hypothesis, the CATE indexed by $X_c$ is constant and equals the ATE, $\delta_0$. The null does not require $\delta_0$ to be zero, nor does it require all individual treatment effects to be identical. Rather, it states that the conditional average effect with respect to $X_c$ is constant. Under the alternative, the conditional average treatment effect varies systematically with at least one component of $X_c$.
This null is the analogue, for the chosen heterogeneity covariates $X_c$, of the treatment-effect homogeneity null considered by crump2008nonparametric, under which CATEs are identical across covariate-defined subpopulations. In contrast to their setting, our procedure allows identification to rely on a potentially high-dimensional vector of controls $X$, with the required nuisance functions estimated flexibly using ML methods.
A common alternative way to study treatment-effect heterogeneity is to estimate the CATEs, because they describe how treatment effects differ by observed characteristics and can therefore help identify which groups benefit more or less from an intervention. Modern examples include tree- and forest-based methods, as well as debiased and orthogonalized approaches for high-dimensional settings ath2016,wager2018estimation,ath2019,semenova2021debiased,fan2022. Our goal is different. Estimating the full conditional effect function is more ambitious than testing whether systematic heterogeneity is present. In addition, inference on CATEs is often used to test whether the treatment effect is zero at particular covariate values, for example \[ \mathcal{H}_0^{\mathrm{zero}}(x): \operatorname{E}\{Y(1)-Y(0)\mid X=x\}=0. \] This is not a general test of heterogeneity. A treatment effect may be constant and nonzero for all groups. In that case, a zero-effect null can be rejected even though there is no treatment-effect heterogeneity. Our test instead asks whether the conditional average effect differs across values of $X_c$, while allowing the average effect itself to be nonzero.
Recent ML-based approaches also provide formal tests for treatment-effect heterogeneity through summaries of treatment-effect variation chernozhukov2023fisherschultzlecturegenericmachine,Imai2025b. Best-linear-predictor (BLP) methods first use ML to construct a score, denoted by $S(X)$, that is intended to predict the CATE $\delta(X)$. This score is not the true CATE, which is unobserved, but a first-stage proxy for heterogeneity. The BLP approach then asks whether actual treatment effects are systematically related to this score by considering the population projection \[ \operatorname{BLP}\{\delta(X)\mid S(X)\} = \beta_1+\beta_2\{S(X)-\operatorname{E}[S(X)]\}. \] The corresponding null hypothesis is \[ \mathcal{H}_0^{\mathrm{BLP}}:\beta_2=0. \] This null asks whether the ML score is linearly related to treatment-effect variation. Thus, a rejection provides evidence of heterogeneity captured by the score $S(X)$.
Sorted-effect methods use the same general logic but test for heterogeneity across groups rather than through a linear slope. Let $G_1,\ldots,G_K$ denote groups formed by sorting observations according to the ML score $S(X)$, and let \[ \delta_k := \operatorname{E}\{\delta(X)\mid X\in G_k\} \] denote the average treatment effect in group $k$. The corresponding homogeneity null is \[ \mathcal{H}_0^{\mathrm{sorted}}:\delta_1=\cdots=\delta_K. \] This null asks whether average treatment effects are the same across the groups induced by the ML score. Both approaches, BLP and sorted-effects methods, therefore test whether heterogeneity is detectable through a particular ML-generated summary, a linear projection in the BLP case or ranked groups in the sorted-effect case. The objective of our test is more basic. We ask whether the data provide evidence that systematic heterogeneity with respect to $X_c$ is present at all. The test therefore does not require heterogeneity to be captured by a linear projection or by a small number of ranked groups.
The null hypothesis considered here also differs from sharp hypotheses of identical unit-level effects studied in randomization-based tests, such as \[ \mathcal{H}_0^{\mathrm{sharp}}: \operatorname{Var}\{Y(1)-Y(0)\}=0 \] heck97,ding2016randomization,ding2019decomposing. Those hypotheses ask whether every individual treatment effect is the same. They can therefore be violated even when treatment-effect variation is idiosyncratic and unrelated to observable characteristics. By contrast, our null is weaker and more targeted. It asks whether the average treatment effect varies in a way that is explainable by the observed covariates $X_c$. This distinction is important in empirical work, where researchers are often interested in whether effects differ across observable and policy-relevant groups rather than whether all individuals have exactly the same effect.
Thus, the proposed test addresses a simple question: do treatment effects vary systematically with $X_c$? The test is omnibus in the sense that it is designed to detect many forms of heterogeneity, including nonlinearities, threshold effects, and interactions among the components of $X_c$. At the same time, it does not require researchers to estimate or interpret the full heterogeneity function. The test can therefore serve as a first step before more detailed analyses of which covariates drive heterogeneity and which groups are most affected.
Consider the potential-outcomes framework introduced in Section (ref). We now formalize the proposed test for treatment-effect heterogeneity under unconfoundedness. Our objective is to construct a test that remains valid when treatment effects are identified using a potentially high-dimensional set of controls and auxiliary quantities are estimated using flexible ML methods.
The development of the test proceeds in four steps. First, we rewrite the null hypothesis of homogeneous treatment effects as a conditional moment restriction involving observable quantities. Second, we transform this conditional restriction into an equivalent collection of unconditional moment conditions that can be tested in practice. Third, we construct locally robust moments that remain valid when the required auxiliary quantities are estimated using ML methods. Fourth, we combine these moments into a test statistic and develop a computationally efficient bootstrap procedure for inference. The resulting procedure provides an omnibus test for treatment-effect heterogeneity that is flexible, computationally tractable, and compatible with a broad class of modern ML estimators.
Throughout this section, we use the notation introduced in Section (ref). The observed outcome for unit $i$ is $Y_i=D_iY_i(1)+(1-D_i)Y_i(0)$, and we observe a random sample $\{(Y_i,D_i,X_i)\}_{i=1}^n$, where $n$ denotes the sample size. We first state the identifying assumptions and then rewrite the null hypothesis in a form that can be tested using observed data.
Assumption (ref) requires that, after conditioning on the observed covariates $X$, treatment assignment is independent of the potential outcomes. For this to be the case, $X$ must contain all variables that jointly influence treatment assignment and potential outcomes.
Let
denote the propensity score, that is, the conditional probability of treatment given the observed covariates.
Assumption (ref) is a common-support condition. It requires that, for every value of $X$, both treated and untreated units occur with positive probability. This ensures that meaningful comparisons between treated and untreated units can be made.
The null hypothesis in (ref) is formulated in terms of potential outcomes, which are only partially observed. We therefore derive an equivalent representation based on quantities that can be estimated from the data. Define
Under Assumption (ref), these conditional outcome regressions identify the corresponding conditional potential outcomes:
Hence, the CATE is identified as
and the ATE is identified as
Let $X_c$ be a subvector of $X$, with $X_c\in\mathbb{R}^c$. The full covariate vector $X$ is used for identification, whereas $X_c$ contains the covariates with respect to which treatment-effect heterogeneity is evaluated. Applying the Law of Iterated Expectations yields
Consequently, the null hypothesis in (ref) is equivalent to
Equation (ref) states that, after removing the average treatment effect $\delta_0$, the remaining treatment-effect variation is conditionally mean independent of $X_c$. This conditional moment restriction is the starting point for the test.
Directly testing (ref) is difficult because it involves a conditional expectation given $X_c$. We therefore use a strategy from the conditional moment literature and replace the conditional restriction by a collection of unconditional moment restrictions. Intuitively, if treatment effects do not vary systematically with $X_c$, then $\{\mu_1(X)-\mu_0(X)-\delta_0\}$ should be uncorrelated with any function of $X_c$.
Assume, without loss of generality, that $X_c$ has bounded support. Then, by bierens2016econometric, the null hypothesis (ref) is equivalent to
where $\varphi$ is an analytic non-polynomial function satisfying $\partial^l\varphi(0)\neq0$ for all $l\in\mathbb{N}$, and $\mathcal{T}$ is a subset of $\mathbb{R}^{\dim(X_c)}$ containing a neighborhood of the origin. See Assumption (ref) for the precise conditions.
To understand this equivalence, note that (ref) is equivalent to
for every square-integrable function $f$. Verifying this condition directly is infeasible because it requires considering an infinite collection of arbitrary functions. The result of bierens2016econometric shows that it is sufficient to focus on the structured family, $\{\varphi(t^T\cdot):t\in\mathcal{T}\}$. This reduction is useful because it leads to a test statistic with a simple closed-form representation.
The moment conditions in (ref) provide a testable characterization of the null hypothesis. A direct implementation, however, would require estimating the unknown functions $\mu_0$ and $\mu_1$. If estimation errors from these first-stage quantities enter the test statistic too strongly, standard inference can fail. To address this issue, we replace the treatment-effect residual in (ref) by a locally robust score. This construction also uses the propensity score $m^D$.
Consider the doubly robust score
The quantity $V$ combines outcome-regression adjustment through $\mu_0$ and $\mu_1$ with propensity-score adjustment through $m^D$. Such doubly robust scores are widely used in modern causal inference because they are stable under flexible first-stage estimation \citetext{see, for example, chernozhukov2018double}.
The score $U$ is constructed so that its conditional expectation recovers the treatment-effect residual appearing in (ref). Specifically,
This follows from
Thus, $U$ captures deviations of the conditional treatment effect from the average treatment effect. Under the null hypothesis of homogeneous treatment effects, these deviations should not vary systematically with $X_c$.
Combining this result with (ref) and the Law of Iterated Expectations yields
Equation (ref) provides the locally robust moment representation used below. Unlike the direct moments in (ref), the moments based on $U$ are locally robust with respect to $\mu_0$, $\mu_1$, and $m^D$. Intuitively, small errors in estimating these functions have only a negligible effect on the resulting moments, which makes the construction well suited for flexible ML-based estimation.
Let us now construct a feasible test statistic based on the robust moments in (ref). If the score $U$ were observed, testing $\mathcal{H}_0$ would be straightforward. In particular, we could measure the extent to which the sample moments \[ \operatorname{E}_n\{U_i\varphi(t^TX_{c,i})\} \] deviate from zero across values of $t$, where $\operatorname{E}_n$ denotes the empirical mean operator, \[ \operatorname{E}_ng(\zeta_i) := \frac{1}{n} \sum_{i=1}^n g(\zeta_i), \] $\zeta_i:=(Y_i,D_i,X_i)$. Following the integrated conditional moment literature, a natural test statistic is
where $\rho$ is a probability measure supported on $\mathcal T$.
The challenge is that $U$ depends on the unknown functions $\mu_0$, $\mu_1$, and $m^D$, which must be estimated from the data. To estimate these quantities while avoiding overfitting, we employ cross-fitting. This is particularly important when the nuisance functions are estimated using flexible ML methods. The basic idea is that each observation is evaluated using models that were estimated without using that observation. This reduces overfitting and helps ensure valid inference.
Formally, let $\{\zeta_i\}_{i=1}^n$ denote the observed sample and partition it into $N$ mutually exclusive folds $\{K_1,\ldots,K_N\}$.\footnote{Each observation $\zeta_i$ will belong to a unique fold and the union of such folds gives the entire sample. Formally, $\cup_{s=1}^N K_s=\{\zeta_i\}_{s=1}^n=:K$ and $K_s\cap K_j=\emptyset$ for $s\neq j$.} Let $K(i)$ denote the fold containing observation $i$, and let $K(-i):=K - K(i)$ denote the remaining observations. We estimate the nuisance functions using only the observations in $K(-i)$. Specifically, let $\widehat m^D_{K(-i)}$, $\widehat\mu_{0,K(-i)}$, and $\widehat\mu_{1,K(-i)}$ denote the corresponding estimators of $m^D$, $\mu_0$, and $\mu_1$.
Using these estimators, we construct
where
Note that the estimator $\widehat\delta$ is not simply the sample analogue of \[ \delta_0=\operatorname{E}\{\mu_1(X)-\mu_0(X)\}. \] A natural plug-in estimator would replace $\mu_0$ and $\mu_1$ with estimated versions and then average the resulting treatment effects. However, when these quantities are estimated using flexible ML methods, estimation errors can propagate into later stages of the analysis and invalidate standard inference procedures. We therefore use the estimator $\widehat\delta$ because it is locally robust to small errors in the estimation of $\mu_0$, $\mu_1$, and $m^D$; see chernozhukov2018double. This robustness property plays a key role in ensuring valid inference while retaining the flexibility of ML estimation.
Our feasible test statistic is
where, for notational convenience, we define $\phi_t(\cdot):=\varphi(t^T\cdot)$.
Intuitively, the statistic measures the extent to which the estimated treatment-effect residuals $\widehat U_i$ are systematically related to functions of the heterogeneity covariates $X_c$. Under the null hypothesis, these empirical moments should be close to zero for all values of $t$. Large values of $S_n$ therefore provide evidence against treatment-effect homogeneity.
Under suitable choices of $\rho$ and $\phi_t$, the statistic admits a simple closed-form representation. In particular, suppose that $\rho$ has a symmetric Fourier transform, denoted by $\mathcal{F}_\rho$, and let \[ \phi_t(\cdot)=\exp(\bm{i}t^T\cdot), \] where $\bm{i}$ denotes the imaginary unit. Then
where $\bm{\widehat U} = (\widehat U_1,\ldots,\widehat U_n)^T$ and $\bm{F_\rho}$ denotes the $n\times n$ matrix with generic element $\mathcal{F}_\rho(X_{c,i}-X_{c,j})$. Equation (ref) is particularly attractive from a computational perspective. The statistic depends only on the estimated quantities $\widehat\mu_0$, $\widehat\mu_1$, and $\widehat m^D$, and does not require any nonparametric smoothing with respect to $X_c$.
We now introduce the regularity conditions required for the asymptotic analysis of the test statistic.
Parts (i) and (ii) of Assumption (ref) are standard regularity conditions. Parts (iii) and (iv) ensure that the continuum of unconditional moments in (ref) provides a complete characterization of the original conditional moment restriction. In particular, they guarantee the equivalence between the null hypothesis and the integrated moment conditions used by the test; see bierens2016econometric. The bounded-support assumption on $X_c$ is without loss of generality, since any unbounded covariate can be transformed through a bounded one-to-one mapping without changing the underlying conditional moment restriction.
For any function $f(\zeta)$, with $\zeta:=(Y,D,X)$, define the $L^2(P)$ norm as \[ \|f\|_{L^2(P)}^2 = \int f(\zeta)^2 dP(\zeta), \] where $P$ denotes the probability law of $\zeta$. Let $\widehat \mu_{1,K-K_j}$ denote the estimator of $\mu_1$ computed on the subsample $K-K_j$, that is, on the sample excluding the observations contained in fold $K_j$. The estimators $\widehat \mu_{0,K-K_j}$ and $\widehat m^D_{K-K_j}$ are defined analogously. Finally, let $\mathcal X$ denote the support of $X$.
Assumption (ref) concerns the quality of the estimators used for the nuisance functions $\mu_0$, $\mu_1$, and $m^D$. Part (i) requires these estimators to be sufficiently accurate in large samples. Formally, the estimation error must converge to zero at a rate faster than $n^{-1/4}$ in the $L^2(P)$ norm. Similar conditions are standard in the double ML literature. For example, the required convergence rates have been established for neural networks by farrell2021deep and for random forests by scornet2015consistency. Part (ii) requires the estimated outcome regressions $\widehat \mu_{d,K-K_j}$ to remain uniformly bounded over the support $\mathcal X$. Part (iii) requires the estimated propensity score $\widehat m^D_{K-K_j}$ to remain bounded away from zero and one. This condition ensures that the inverse-probability weights appearing in (ref) remain well behaved and do not become excessively large.
The following proposition characterizes the large-sample behavior of the test statistic $S_n$. To state the result, let
and recall that $\zeta_i=(Y_i,D_i,X_i)$.
Part (a) of Proposition (ref) establishes the asymptotic distribution of the test statistic under the null hypothesis of homogeneous treatment effects. Part (b) shows that the test statistic diverges under the alternative hypothesis, implying that the proposed test is consistent against fixed alternatives. The proof is provided in the supplementary material.
Although the limiting distribution in Part (a) is well defined, it depends on unknown features of the data-generating process and is therefore difficult to use directly for inference. In the next subsection, we develop a bootstrap procedure that provides a practical way to approximate the distribution of the test statistic and obtain critical values.
To obtain the critical values necessary for testing, we bootstrap our test statistic. Let $\{\xi_i\}_{i=1}^n$ be a sequence of iid bootstrap weights independent from the sample data, with $\operatorname{E}\{\xi\}=1$ and $\operatorname{Var}\{\xi\}=1$. We define the bootstrap counterpart of our statistic as
where
When $\rho$ has a symmetric Fourier transform, $\mathcal{F}_\rho$, and $\phi_t(\cdot)=\exp(\bm{i}t^T\cdot)$, the bootstrapped statistic $S_n^*$ has a simple expression given by
where $\bm{\widehat U^*}=(\widehat U_1^*,\ldots,\widehat U_n^*)^T$ and $\bm{\widetilde U^*}=(\widetilde U_1^*,\ldots,\widetilde U_n^*)^T$.
In practice, to obtain the critical value necessary for testing, we can apply the following steps:
The above procedure is based on vector and quadratic forms computations, so it is straightforward to implement. Note that it does not require to compute the ML estimators at each bootstrap iteration. This is computationally convenient. \\ We next establish the validity of this bootstrap procedure. Let us denote with $\Pr^*$ the probability that considers as random only the bootstrap weights $\{\xi_i\}_{i=1}^n$, that is, the probability conditional to the sample data. We define formally the $1-\alpha$ quantile of the bootstrap distribution of $S_n^*$, $\widehat q_{1-\alpha}$, as
The above result is shown in the supplementary material. Part (a) of the above proposition shows that our bootstrap test controls the size asymptotically. Part (b) shows that it can detect violations of the null hypothesis.
Section (ref) develops our testing procedure under identification by unconfoundedness. The same ideas, however, extend naturally to other widely used causal inference frameworks. In this section, we consider Difference-in-Differences (DiD) models with panel and repeated cross section data, as well as Local Average Treatment Effects (LATE) identified through instrumental variables.
The construction closely parallels that developed in Section (ref). In each case, we reformulate the null hypothesis of homogeneous treatment effects in terms of observable quantities, construct locally robust moment conditions, and build a test statistic based on cross-fitted estimators of the required nuisance functions. The resulting procedures retain the same ability to accommodate flexible ML estimators while maintaining valid inference.
To avoid repetition, we focus on the framework-specific assumptions and formal results. Detailed derivations of the moment conditions, test statistics, asymptotic distributions, and bootstrap procedures are provided in the Appendix.
We first consider a Difference-in-Differences (DiD) design with two periods, $t=0,1$, and panel data. Let $G_i\in\{0,1\}$ indicate whether unit $i$ belongs to the treated group. Treated units receive treatment only in period $1$, while control units are never treated. Let $Y_{i,t}(d)$ denote the potential outcome in period $t$ under treatment status $d\in\{0,1\}$, and let $Y_{i,t}$ denote the observed outcome. We observe
where $X_i$ is a vector of pre-treatment covariates. Define
Our objective is to test whether the average treatment effect on the treated varies systematically with the heterogeneity covariates $X_c\subseteq X$. Formally, we consider
where $\mathcal{H}_0^{DiD,c}$ denotes the logical complement of $\mathcal{H}_0^{DiD}$.
Identification relies on the standard conditional parallel-trends assumption.
Assumption (ref) requires that, after conditioning on $X$, treated and control units would have experienced the same average outcome trend in the absence of treatment. Assumption (ref) ensures that treated units exist and that suitable control units are available across the relevant support of $X$. Assumption (ref) ensures that the conditional treatment effect for the treated is well defined for the values of $X_c$ considered in the test.
Define
Following the same logic as in Section (ref), the null hypothesis can be written as a collection of locally robust moment restrictions. Let
Then
A detailed derivation of Equation (ref) is provided in Appendix (ref). As in Section (ref), we estimate the nuisance functions using cross-fitting. Let $\widehat \mu^{Y_1-Y_0}_{0,K(-i)}$ and $\widehat m^G_{K(-i)}$ denote estimators obtained without using the fold containing observation $i$. Define
The resulting test statistic is
where $\phi_t(X_c)=\varphi(t^TX_c)$. As in Section (ref), choosing
yields the computationally convenient representation
where
and $\bm{F_\rho}$ is defined in Section (ref). We impose the following regularity conditions. Let $\zeta:=(Y_0,Y_1,G,X)$, let $P$ denote the distribution of $\zeta$, and let $\mathcal X$ denote the support of $X$. The norm $\|\cdot\|_{L^2(P)}$ is defined as in Section (ref).
Assumptions (ref) and (ref) are the DiD analogues of the regularity conditions imposed in the unconfoundedness setting. They ensure that the integrated moment representation characterizes the null hypothesis and that the nuisance estimators are sufficiently accurate for locally robust inference. As before, the required convergence rates can be satisfied by a broad class of modern ML estimators.
Part (a) establishes the asymptotic distribution of the test statistic under the null hypothesis of homogeneous average treatment effects on the treated. Part (b) shows that the statistic diverges under the alternative, implying that the test is consistent against fixed alternatives. The proof is provided in the supplementary material.
We obtain critical values using a bootstrap procedure analogous to that in Section (ref). The nuisance functions are estimated only once and are kept fixed across bootstrap iterations. The bootstrap recomputes the low-dimensional quantities that enter the DiD score, namely the treated share and the ATT estimate. Let $\{\xi_i\}_{i=1}^n$ be an iid sequence of bootstrap weights independent from the sample data, with $\operatorname{E}\{\xi\}=\operatorname{Var}\{\xi\}=1$. We define
Then, the bootstrap counterpart of the test statistic $S_n^{DiD}$ is
By setting $\phi_t(\cdot)=\exp(\textbf{i} t^T\cdot)$, the bootstrapped statistic $S_n^{DID *}$ can be computed similarly to $S_n^{DiD}$, as follows
where $\bm{\widehat U^{DiD*}}=(\widehat U^{DiD*}_1,\ldots,\widehat U^{DiD*}_n)^T$ and $\bm{\widehat U^{DiD}}=(\widehat U^{DiD}_1,\ldots,\widehat U^{DiD}_n)^T$. The bootstrap test can be implemented by the Monte Carlo procedure described in Section (ref).
Let $\widehat q_{1-\alpha}^{DiD}$ denote the empirical $(1-\alpha)$ quantile of the bootstrap statistics generated by the procedure above. The following proposition establishes the validity of this design-specific bootstrap critical value.
Proposition (ref) establishes that the bootstrap test has asymptotically correct size under the null hypothesis and is consistent under fixed alternatives. The proof is provided in the supplementary material.
We next consider a Difference-in-Differences design with repeated cross sections. In contrast to the panel-data setting, the same unit is not observed before and after treatment. Instead, the data contain different units observed either in the pre-treatment period or in the post-treatment period. Let $T_i\in\{0,1\}$ indicate whether unit $i$ is observed in the post-treatment period, and let $G_i\in\{0,1\}$ indicate whether unit $i$ belongs to the treated group. We observe
where $Y_i$ is the observed outcome and $X_i$ is a vector of pre-treatment covariates. We maintain the conditional parallel-trends and support conditions from Assumptions (ref), (ref), and (ref). As before, let
The repeated-cross-section sampling scheme is formalized as follows.
Assumption (ref) states that the observed data combine two cross sections: one from the pre-treatment period and one from the post-treatment period. In particular, $\lambda=\Pr(T=1)$ denotes the population share of observations from the post-treatment period.
Our objective is again to test whether the treatment effect on the treated varies systematically with $X_c$. Formally, we consider
where $\mathcal{H}_0^{DiDrcs,c}$ denotes the logical complement of $\mathcal{H}_0^{DiDrcs}$.
Following the same logic as in Section (ref), the null hypothesis can be represented through locally robust moments. Define
Let
Then
A detailed derivation of Equation (ref) is provided in Appendix (ref). The score $U^{DiDrcs}$ plays the same role as $U^{DiD}$ in the panel-data case. It converts the null hypothesis into moments that can be estimated using observed data and remains locally robust with respect to the first-stage nuisance estimates. The additional feature in the repeated-cross-section setting is that the score depends on the sampling share $\lambda$, which also has to be estimated.
We estimate $\mu_0^Y$, $\mu_0^{YT}$, and $m^G$ by cross-fitting. Let $\widehat\mu^Y_{0,K(-i)}$, $\widehat\mu^{YT}_{0,K(-i)}$, and $\widehat m^G_{K(-i)}$ denote estimators obtained without using the fold containing observation $i$. Define
The resulting test statistic is
where $\phi_t(X_c)=\varphi(t^TX_c)$. As before, choosing
gives the closed-form expression
where
We impose the following regularity conditions. Let $\zeta:=(Y,G,T,X)$, let $P$ denote the distribution of $\zeta$, and let $\mathcal X$ denote the support of $X$. The norm $\|\cdot\|_{L^2(P)}$ is defined as in Section (ref).
Assumptions (ref) and (ref) are the repeated-cross-section analogues of the regularity conditions used above. They ensure that the integrated moment representation remains valid, that the nuisance estimators are sufficiently accurate for locally robust inference, and that the inverse-probability weights remain stable.
Part (a) establishes the asymptotic distribution of the test statistic under the null hypothesis of homogeneous treatment effects on the treated, accounting for estimation of $\lambda$ and $\pi_0$. Part (b) shows that the statistic diverges under the alternative, implying consistency against fixed alternatives. The proof is provided in the supplementary material.
We obtain critical values using a bootstrap procedure analogous to that in Section (ref). The nuisance functions are estimated only once and are kept fixed across bootstrap iterations. Compared with the panel-data case, the bootstrap must also account for the estimation of the post-treatment sampling share $\lambda$, which enters the score $W(\lambda)$. Let $\{\xi_i\}_{i=1}^n$ be a sequence of iid bootstrap weights independent from the sample data, with $\operatorname{E}\{\xi\}=1$ and $\text{Var}\{\xi\}=1$. We define
and
Then, the bootstrap counterpart of our statistic is
By setting $\phi_t(\cdot)=\exp(\textbf{i} t^T\cdot)$, we obtain the following simple expression for the bootstrapped statistic
where $[\bm{\widehat U ^{DiDrcs*}}=(\widehat U ^{DiDrcs*}_1,\ldots,\widehat U ^{DiDrcs*}_n)^T$ and $\bm{\widetilde U ^{DiDrcs*}}=(\widetilde U ^{DiDrcs*}_1,\ldots,\widetilde U ^{DiDrcs*}_n)^T$. The bootstrap scheme can then be implemented by the Monte Carlo procedure described in Section (ref). We remark that the above bootstrap scheme is not based on a multiplier bootstrap version of the influence function representation of the empirical process at the basis of the test statistic. Such a version would have a rather complicated expression. The bootstrap scheme in the present section can be seen as a hybrid between a weighed bootstrap and a multiplier bootstrap.
Let $\widehat q_{1-\alpha}^{DiDrcs}$ denote the empirical $(1-\alpha)$ quantile of the bootstrap statistics generated by the procedure above. The following proposition establishes the validity of this repeated-cross-section bootstrap critical value.
Proposition (ref) establishes that the bootstrap test has asymptotically correct size under the null hypothesis and is consistent under fixed alternatives. The proof is provided in the supplementary material.
We finally consider settings with a binary instrumental variable. In contrast to the unconfoundedness case, treatment assignment may be endogenous even after conditioning on $X$. Identification instead relies on an instrument $Z\in\{0,1\}$ that affects treatment take-up while satisfying the conditions stated below. Let $Y(d)$ denote the potential outcome under treatment status $d\in\{0,1\}$, and let $D(z)\in\{0,1\}$ denote the potential treatment status under instrument value $z\in\{0,1\}$. We observe an iid sample
where $D_i=D_i(Z_i)$ and
Define the instrument propensity score
The target parameter is the Local Average Treatment Effect (LATE), that is, the average treatment effect for compliers, defined as units whose treatment status is shifted by the instrument. Our objective is to test whether this local treatment effect varies systematically with the heterogeneity covariates $X_c\subseteq X$.
Assumption (ref) requires the instrument to be conditionally exogenous given $X$. Assumption (ref) rules out defiers. Assumption (ref) ensures that compliers exist for each value of $X_c$, so that the conditional LATE is well defined.
Our null hypothesis is
where $\mathcal{H}_0^{LATE,c}$ denotes the logical complement of $\mathcal{H}_0^{LATE}$. Under the null hypothesis, the conditional LATE is constant in $X_c$ and coincides with the unconditional LATE, $\delta_0$.
We also impose a common-support condition for the instrument.
Define, for $z\in\{0,1\}$,
Following the same logic as in Section (ref), the null hypothesis can be represented through locally robust moments. Let
Here, $V^Y$ is the locally robust analogue of the reduced-form effect of the instrument on the outcome, while $V^D$ is the corresponding first-stage effect on treatment take-up. The null hypothesis is equivalently characterized by
We estimate $m^Z$, $\mu^Y_{Z=z}$, and $\mu^D_{Z=z}$ using the same cross-fitting procedure as in Section (ref). Given the fitted values \[ \left\{ \widehat \mu^Y_{Z=1,K(-i)}(X_i), \widehat \mu^Y_{Z=0,K(-i)}(X_i), \widehat \mu^D_{Z=1,K(-i)}(X_i), \widehat \mu^D_{Z=0,K(-i)}(X_i), \widehat m^Z_{K(-i)}(X_i) \right\}_{i=1}^n, \] define
The estimator $\widehat\delta^{LATE}$ is the locally robust analogue of the unconditional LATE estimator. It equals the ratio of the reduced-form and first-stage components, while allowing both components to be estimated flexibly.
The resulting test statistic is
where $\phi_t(X_c):=\varphi(t^TX_c)$. As before, choosing
yields the closed-form expression
where
We impose the following regularity conditions. Let $\zeta:=(Y,D,Z,X)$, let $P$ denote the distribution of $\zeta$, and let $\mathcal X$ denote the support of $X$. The norm $\|\cdot\|_{L^2(P)}$ is defined as in Section (ref).
Assumptions (ref) and (ref) are the LATE analogues of the regularity conditions used above. They ensure that the reduced-form and first-stage nuisance functions are estimated accurately enough and that the estimated instrument propensity score remains bounded away from zero and one.
Part (a) establishes the asymptotic distribution of the test statistic under the null hypothesis of homogeneous local treatment effects. Part (b) shows that the statistic diverges under the alternative, implying consistency against fixed alternatives. The proof is provided in the supplementary material.
We obtain critical values using a bootstrap procedure analogous to that in Section (ref). The nuisance functions are estimated only once and are kept fixed across bootstrap iterations. Since the LATE parameter is a ratio, the bootstrap recomputes the reduced-form and first-stage components in each iteration. Let $\{\xi_i\}_{i=1}^n$ be an iid sequence of bootstrap weights independent from the sample data, with $\operatorname{E}\{\xi\}=\operatorname{Var}\{\xi\}=1$. We define
The bootstrap counterpart of the test statistic is
By setting $\phi_t(\cdot)=\exp(\textbf{i}\cdot)$, the above statistic has the following simple expression
where $\bm{\widehat U^{LATE*}}:=(\widehat U_1^{LATE*},\ldots,\widehat U_n^{LATE*})^T$ and $\bm{\widehat U^{LATE}}:=(\widehat U_1^{LATE},\ldots,\widehat U_n^{LATE})^T$. In practice, to implement the bootstrap test, we can use the Monte Carlo procedure described in Section (ref).
Let $\widehat q_{1-\alpha}^{LATE}$ denote the empirical $(1-\alpha)$ quantile of the bootstrap statistics generated by the procedure above. The following proposition establishes the validity of this LATE bootstrap critical value.
Proposition (ref) establishes that the bootstrap test has asymptotically correct size under the null hypothesis and is consistent under fixed alternatives. The proof is provided in the supplementary material.
This section evaluates the finite-sample performance of our bootstrap test for DiD with repeated cross-sections. Appendix (ref) provides additional simulations for the ATE under conditional unconfoundedness.
We use a data-generating process (DGP) similar to that in chang2020double, adapted to feature conditional, but not unconditional, parallel trends, different levels of treatment-effect heterogeneity, and common support.
We first simulate $X\sim \mathcal{N}(0,I_{10})$ and set $\theta_1=\left(1,\frac{1}{2},\frac{1}{3},\frac{1}{4},\frac{1}{5},0,0,0,0,0\right)'$. We then define
Thus, assignment to the treated group depends on the first five covariates, with the first covariate receiving the largest weight. By contrast, whether a unit is observed in the post-treatment period is independent of its covariates.
Next, we generate $\varepsilon_{i,k} \sim \mathcal{N}(0,0.01)$ for $k \in \{1,2,3\}$ and set $\theta_2=\left(\frac{3}{2},1,\frac{5}{6},\frac{3}{4},\frac{7}{10},0,0,0,0,0\right)'$. Potential and observed outcomes are generated as follows:
When $\xi\neq 0$, the parallel-trends assumption holds only conditionally. We therefore fix $\xi=1/2$. The parameter $\gamma$ controls the degree of heterogeneity in the ATT. In particular, for any subset $X_c$, if $\gamma=0$, then \[ \operatorname{E}\left[Y_{i,1}(1)- Y_{i,1}(0)\mid G_i=1, X_{c,i}\right]=3, \] so the null hypothesis of no ATT heterogeneity with respect to $X_c$ holds. If $\gamma\neq 0$ and $\theta_{2,c}$ denotes the subvector of $\theta_2$ corresponding to the covariates in $X_c$, then \[ \operatorname{E}\left[Y_{i,1}(1)- Y_{i,1}(0)\mid G_i=1, X_{c,i}\right] = 3 + \gamma\left(1+1.5X_{c,i}'\theta_{2,c}\right). \] Thus, whenever $X_c$ contains at least one of the first five covariates in $X$, the null hypothesis is violated.
We evaluate rejection rates under the null and alternative hypotheses across 1,000 Monte Carlo simulations. We consider sample sizes $n\in\{500,1000\}$ and use $B=499$ bootstrap iterations. We test for heterogeneity with respect to subsets $X_c$ containing the first $d\in\{2,5\}$ components of $X$, that is, $X_c=(X_1,X_2)$ or $X_c=(X_1,\ldots,X_5)$. For nuisance estimation, we use lasso and random forests. Cross-fitting is implemented with five equally sized random folds. We set $\phi_t(\cdot)=\exp(\bm{i}t^\top\cdot)$ and let $\rho(\cdot)$ be the standard Gaussian density, so that $\mathcal{F}_\rho$ corresponds to a Gaussian kernel with bandwidth $h=1$.
We first evaluate the validity of the bootstrap test under the null hypothesis, corresponding to $\gamma=0$. Table (ref) reports rejection rates under the null for significance levels of 5% and 10%. For $d=2$, the test is close to correctly sized. For $d=5$, it becomes slightly conservative. This pattern is stable across sample sizes and nuisance-estimation methods.
We next examine power under the alternative by evaluating rejection rates for different values of $\gamma$. Figure (ref) displays the corresponding power curves at the 5% significance level. As expected, rejection rates increase with $\gamma$ and approach one. Power is lower for $n=500$, but including more heterogeneity covariates, $d=5$, does not substantially affect rejection rates. Finally, lasso and random forests yield nearly identical results.
We illustrate the proposed test using two established empirical applications. The first studies heterogeneity in the effect of 401(k) eligibility on household financial assets under conditional unconfoundedness. The second studies heterogeneity in the effect of tariff reductions on corruption in a difference-in-differences setting with repeated cross-sections. In both applications, our goal is not to revisit the original identification arguments in detail, but to demonstrate how the proposed test can be implemented with flexible ML-based nuisance estimation.
The first application considers the effect of eligibility for a 401(k) retirement savings plan on household financial asset accumulation. We follow the empirical setting of poterba_401k_1995, which was later revisited in the double ML framework of chernozhukov2018double. Eligibility for a 401(k) plan is not randomly assigned, since workers employed at firms offering such plans may differ systematically from other workers. Following this literature, we treat eligibility as conditionally exogenous after controlling for income and demographic characteristics related to job choice and saving behavior.
The data come from the 1991 Survey of Income and Program Participation (SIPP), using the sample analyzed in chernozhukov2018double. The outcome is net financial assets, defined as household financial assets net of non-mortgage debt, and the treatment is an indicator for 401(k) eligibility. The covariates include age, income, education, family size, marital status, two-earner status, defined benefit pension status, IRA participation, and home ownership. We estimate the nuisance functions using lasso and random forests.
Figure (ref) reports p-values as the set of heterogeneity variables is expanded sequentially, beginning with age and income and then adding further demographic and financial characteristics. The results provide evidence of treatment-effect heterogeneity across observable household characteristics for both nuisance-estimation methods. The p-values are generally smaller in lower-dimensional specifications and increase as more heterogeneity variables are included, especially with lasso. This pattern is consistent with the simulation evidence, where increasing the dimension of $X_c$ reduces power. The random forest results are more stable across specifications.
The findings are economically plausible. The effect of access to a tax-advantaged retirement account may differ by income, age, and pre-existing savings opportunities. The application therefore illustrates how the proposed procedure can detect systematic heterogeneity in an observational setting while allowing for flexible adjustment for confounding.
The second application studies treatment-effect heterogeneity in the difference-in-differences setting analyzed by sequeira_corruption_2016. The application examines how reductions in import tariffs affected corruption along the South Africa--Mozambique trade corridor. The empirical design exploits the Southern African Development Community (SADC) Trade Protocol, which gradually reduced tariffs on South African imports into Mozambique and generated a large tariff reduction in 2008.
We use the sample from chang2020double, which combines disaggregated customs data with audit data on shipments passing through the port of Maputo and the Mozambique--South Africa border crossing. The outcome is the log bribe amount associated with a shipment. Treatment is defined by whether a product category experienced a tariff reduction in 2008, and the post-treatment indicator captures the period after the reform. The covariates include shipment characteristics, product categories, inspection indicators, shipment values, tariff levels, and other controls related to the import process. We estimate the nuisance functions using lasso and random forests in the repeated-cross-section DiD framework.
Figure (ref) reports p-values as the set of heterogeneity variables is expanded sequentially. The random forest specification detects significant heterogeneity across several specifications, whereas the lasso specification generally does not. This difference suggests that flexible nuisance estimation may be important in this application, where the relevant adjustment functions may involve nonlinearities or interactions. As in the 401(k) application, p-values tend to increase as the dimension of $X_c$ grows, in line with the finite-sample patterns documented in the simulations.
The evidence is consistent with economically meaningful heterogeneity. The effect of tariff reductions on bribe payments may vary with shipment characteristics, product types, tariff exposure, and monitoring intensity. The application shows that the proposed test can be implemented in repeated-cross-section DiD designs and can reveal systematic heterogeneity while preserving the use of flexible ML-based nuisance estimation.
This paper develops a nonparametric omnibus test for assessing whether treatment effects vary systematically with observed characteristics. The procedure is designed for settings in which identification may rely on a rich, potentially high-dimensional set of controls, while heterogeneity is evaluated with respect to a lower-dimensional subset of covariates. The test is compatible with flexible nuisance estimation, including modern ML methods, and is based on locally robust moments that reduce sensitivity to first-stage estimation error.
The framework applies to several common empirical designs. We first develop the test for average treatment effects under conditional unconfoundedness and then extend it to treatment effects on the treated in difference-in-differences designs and to local average treatment effects in instrumental-variables settings. For each design, inference is implemented using a computationally efficient bootstrap that keeps the nuisance estimates fixed across bootstrap iterations.
Monte Carlo simulations show that the proposed test attains near-nominal size under the null and has good power against heterogeneous alternatives. The empirical applications illustrate how the procedure can be used to detect systematic treatment-effect heterogeneity in observational and difference-in-differences settings with flexible ML-based nuisance estimation.
Several extensions are left for future work. In particular, it would be useful to study settings in which both identification and heterogeneity are governed by high-dimensional covariates, as well as data structures beyond the designs considered here. Overall, the proposed method provides a practical tool for researchers who want to formally assess whether treatment effects differ across observable and policy-relevant characteristics.