Dimension Reduction for Conditional Density Estimation with Applications to High-Dimensional Causal Inference

Dimension Reduction for Conditional Density Estimation with Applications to High-Dimensional Causal Inference

\onehalfspacing

\onehalfspacing

Introduction

Understanding the relationship between a response variable and a set of covariates is of central interest in statistics and econometrics. The conditional probability density function fully characterizes this relationship. Beyond being a fundamental object of interest, the conditional density plays a pivotal role in a wide range of estimation and inference methods. For instance, in causal inference, the propensity score--the conditional density of a binary treatment on covariates--is central to inverse propensity weighting (IPW; e.g., Hirano2003), doubly robust, and propensity score matching (e.g., abadie2016matching) estimators that are designed to adjust for confounding and mitigate selection bias in observational studies. The literature on this topic is vast, and we refer readers to recent comprehensive treatments by ImbensRubin, ding2024first, and chernozhukov2024applied. The importance of conditional densities is also pronounced in modern econometrics, with applications including the estimation of structural models (e.g., GuerreEtal, Matzkin2013, PerrigneVuong), nonparametric estimation of nonseparable models (e.g., AltonjiMatzkin, BlundellEtal2017), and semiparametric estimation of discrete choice models with endogeneity (e.g., Lewbel2000), among others.

We investigate the estimation of conditional densities in high-dimensional settings, a task often challenged by the “curse of dimensionality”. To address this issue, existing approaches commonly impose strong distributional or functional form assumptions (e.g., MaZhu2013). However, in real-world data, these assumptions can be rather restrictive. In this paper, we propose an alternative approach that does not depend on these constraints. Instead, we introduce a “sparsity” condition, where only a small subset of covariates meaningfully affects the distribution of the response variable. While this sparsity condition limits the applicability of our method in cases with many weak signals, it is inspired by the observation in the well-known work of HallRacineLi that irrelevant covariates are surprisingly common in practice. Nonetheless, the primary advantage of our method is that it does not rely on distributional or functional form assumptions.

This paper is further motivated by a practical challenge that arises in the doubly robust estimation of average treatment effects (ATE) in high-dimensional settings. While practitioners may tend to include many covariates to strengthen the ignorability assumption, doing so undermines the overlap condition (because the treatment status becomes more predictable, and the propensity score is therefore more likely to approach 0 or 1) and exacerbates the curse of dimensionality for nonparametric estimators. This, in turn, can hinder the application of modern double/debiased machine learning (DML) methods for robust causal inference. Our identification results, summarized in Proposition (ref), offer a direct solution to this challenge. Under the above sparsity condition, these results establish a “double dimension reduction”: recovering the sparsity structure of the treatment variable not only simplifies the propensity score model but also reduces the dimensionality of the outcome regression models through a refined ignorability condition. This insight, combined with our proposed dimension-reduced conditional density estimator, enables the identification and use of a lower-dimensional set of truly relevant covariates. As a result, it enhances the plausibility of the overlap condition and facilitates flexible nonparametric estimation of both the propensity score and outcome regressions in high-dimensional contexts.

To set the stage, let $Y$ denote a response variable, $\boldsymbol{W}$ be a set of pre-selected covariates known to be relevant to $Y$, and $\boldsymbol{X}$ be a high-dimensional vector of candidate covariates. The covariates in $\boldsymbol{W}$ are typically selected based on domain knowledge, theoretical models, or prior selection in earlier analysis stages; $\boldsymbol{W}$ may also be empty.\footnote{If data sparsity or model interpretability is not a key concern, $\boldsymbol{W}$ may also include practitioner-constructed sufficient statistics, such as leading principal components.} Beyond its intrinsic interest, including well-chosen $\boldsymbol{W}$ can simplify the dimension reduction task when many covariates in $\boldsymbol{X}$ influence $Y$ only through $\boldsymbol{W}$, as illustrated in our empirical application. Our framework allows the dimension of $\boldsymbol{X}$ to grow at a polynomial rate. The goal is to identify the covariates among $\boldsymbol{X}$ that are relevant to $Y$ given $\boldsymbol{W}$. To this end, we propose a new, model-free measure of conditional dependence between $Y$ and each covariate in $\boldsymbol{X}$, conditional on $\boldsymbol{W}$. This measure is non-negative and equals zero if and only if $Y$ is independent of the covariate given $\boldsymbol{W}$. This conditional dependence measure underpins our dimension reduction approach and constitutes one of the core contributions of the paper.

Building on this measure, we develop a novel, two-stage procedure for dimension-reduced conditional density estimation. In the first stage, we adapt the “sure independence screening” (SIS) framework of FanLv2008 to design a screening procedure that efficiently identifies relevant covariates in a high-dimensional setting. Specifically, we compute the conditional dependence measure for each candidate covariate individually, rank them, and retain those that exhibit the strongest conditional dependence with $Y$ given $\boldsymbol{W}$. This screening step is particularly well-suited for nonparametric density estimations, as directly measuring dependence or regressing $Y$ on high-dimensional covariates without parametric or functional form assumptions is generally infeasible.

In the second stage, we refine the variable selection using a modified cross-validation (CV) procedure, building on HallRacineLi, to determine the final set of covariates in an objective and data-driven manner. The key insight of the CV-based refinement is that, in kernel conditional density estimation, the optimal bandwidths for irrelevant covariates, chosen to minimize the integrated squared error (ISE), tend to diverge to their upper extremes, effectively smoothing them out of the final density estimate. Our modification significantly reduces the computational burden of the original algorithm by HallRacineLi via restricting the search space, thereby improving numerical stability and computational efficiency, as demonstrated in our simulation studies. The final conditional density estimate is obtained using the selected covariates and the optimal bandwidths from this refinement step.

We demonstrate the practical utility of our dimension reduction method by applying it to the doubly robust estimation of ATE. Our primary contribution in this section is the identification results summarized in Proposition (ref). Simulation results show that the doubly robust ATE estimators that incorporate the insights of Proposition (ref) and leverage our dimension-reduced propensity score estimator exhibit satisfactory finite-sample performance, in sharp contrast to estimators that use the full high-dimensional covariate set. We further illustrate this approach through revisiting the empirical analysis of the effect of 401(k) eligibility on asset accumulation, using rich longitudinal data from the 1996 Survey of Income and Program Participation (SIPP).

Our measure of conditional dependence is a novel generalization of the Kolmogorov-Smirnov (KS) distance-based measure proposed by MaiZou2015, and it differs from theirs in several important aspects. First, while their measure is designed exclusively for assessing unconditional (in)dependence (i.e., when $\boldsymbol{W}$ is empty), our measure naturally extends to conditional settings. Second, by introducing fixed reference values when computing the KS distance, our measure substantially reduces computational burden; see Remark (ref) for further discussion.

Two other widely recognized dependence measures in the literature are the (conditional) distance correlation introduced by Szekely_et_al and WangEtal2015, and the (conditional) rank-based correlation proposed by Chatterjee2021 and azadkia2021simple. In terms of computation, our method is significantly faster than distance correlation but slower than the rank-based approach, as detailed in Remark (ref). As expected, our simulations indicate that the statistical power of our method in small samples lies between the two: lower than that of distance correlation, but higher than that of the rank-based approach. We therefore view our method as particularly useful for moderate sample sizes (e.g., in the thousands to tens of thousands), striking a balance between computational efficiency and small sample statistical power.

This paper also contributes to the recent advances in estimating conditional density functions using a variety of machine learning (ML) methods; see, for example, izbicki_converting_2017, rothfuss2019conditional, gao2022lincde, and the references therein. However, all these studies essentially assume that the researcher has prior knowledge of the finite-dimensional set of relevant covariates. In contrast, our approach provides a data-driven solution to identifying relevant covariates in high-dimensional settings. Additionally, our method complements the rapidly growing literature on causal inference using DML methods, where the propensity score or conditional density is often a key input for constructing robust estimators. Examples include farrell2015robust, ChernozhukovEtal2018, chang2020double, FarrellEtal2021, zhang2024continuous, and haddad2024difference, among many others reviewed in chernozhukov2024applied.

The remainder of this paper is organized as follows. Section (ref) introduces our new conditional dependence measure, the variable screening and CV refinement procedures, and the dimension-reduced conditional density estimator. Section (ref) demonstrates how the proposed methods can be applied to estimate ATE using doubly robust estimators, with a focus on dimension-reduced propensity score estimation. Section (ref) evaluates the finite-sample performance of our approach through comprehensive simulation studies, where we also provide practical implementation guidance. Section (ref) demonstrates the application of our approach to the analysis of the effect of 401(k) eligibility on savings. Section (ref) concludes the paper.

Supporting materials are provided in the appendices. Appendix (ref) contains additional results and further discussion of related methods in the literature. Appendix (ref) presents a computationally efficient algorithm to calculate our proposed dependence measure. Proofs of the main theoretical results are provided in Appendix (ref), while technical lemmas and their proofs are deferred to Appendix (ref). Finally, Appendix (ref) collects tables summarizing the results of our Monte Carlo simulations and empirical illustration.

\noindentNotation. All vectors are column vectors. Capital letters denote random elements, and the corresponding lowercase letters denote their realizations. Boldface letters represent vectors, while regular (non-bold) letters represent scalars. The notation $\textrm{dim}(\boldsymbol{z})$ denotes the dimension of a vector $\boldsymbol{z}$. Superscripts “$\mathsf{c}$” and “$\mathsf{d}$” are used to indicate continuous and discrete random variables, respectively. We use $\Pr(\cdot)$ and $\mathbb{E}[\cdot]$ to denote probability and expectation, respectively. The function $\mathbf{1}(\cdot)$ is the indicator function, which equals one if the event in the parentheses occurs, and zero otherwise. For two random vectors $\boldsymbol{U}$ and $\boldsymbol{V}$, the notation $\boldsymbol{U}|\cdot \overset{d}{\sim} \boldsymbol{V}|\cdot$ indicates that $\boldsymbol{U}$ and $\boldsymbol{V}$ have identical distributions conditional on $\cdot$. The notation $\boldsymbol{U} \perp \boldsymbol{V}$ ($\boldsymbol{U} \perp \boldsymbol{V}|\cdot$) denotes stochastic independence (conditional on $\cdot$). We use $F(\boldsymbol{u}|\boldsymbol{v})$ and $f(\boldsymbol{u}|\boldsymbol{v})$ to denote the conditional cumulative distribution function (CDF) and probability density function (PDF) of $\boldsymbol{U}$ given $\boldsymbol{V}$, respectively. As $n \to \infty$, the notations $\overset{P}{\rightarrow}$ and $\overset{d}{\rightarrow}$ denote convergence in probability and convergence, respectively. The symbol $C$ represents a generic positive constant whose value may vary from line to line. For deterministic sequences $\left\{ a_{n}\right\} _{n=1}^{\infty}$ and $\left\{ b_{n}\right\} _{n=1}^{\infty}$, we write $a_{n}\propto b_{n}$ (or $a_n=O(b_n)$) to mean that $0<C_{1}\leq\lim\inf_{n\rightarrow\infty}\left\vert a_{n}/b_{n}\right\vert \leq\lim\sup_{n\rightarrow\infty}\left\vert a_{n}/b_{n}\right\vert \leq C_{2}<\infty$ for some $C_1, C_2 > 0$. We write $a_n \lesssim b_n$ if $\limsup_{n \to \infty} |a_n / b_n| \leq C$ for some $C>0$, and $a_n \gtrsim b_n$ if $b_n \lesssim a_n$. We write $a_n \ll b_n$ (or $a_n = o(b_n)$) if $a_n / b_n \to 0$, and $a_n \gg b_n$ if $b_n \ll a_n$. The symbol $\setminus$ denotes set difference.

Dimension Reduction for Conditional Densities

This section develops our dimension reduction methodology for conditional density estimation. We consider a setting with a response variable $Y$, a set of pre-selected covariates $\boldsymbol{W}$ known to be relevant, and a high-dimensional vector of candidate covariates $\boldsymbol{X} = (X_1, \dots, X_p)'$. The goal is to estimate the conditional density $f(Y|\boldsymbol{X}, \boldsymbol{W})$ under a sparsity assumption: only a small, unknown subset of $\boldsymbol{X}$ is truly relevant for this density.

To achieve this, we first develop a new measure of conditional dependence (Section (ref)) and apply it within an SIS procedure to screen for relevant variables (Section (ref)). We then propose a CV refinement step and construct the post-selection conditional density estimator (Section (ref)). In Section (ref), we discuss how this two-stage procedure can be specifically adapted for propensity score estimation. In what follows, we decompose $\boldsymbol{W}$ into its continuous and discrete components as $\boldsymbol{W} = \left( \boldsymbol{W}^{\mathsf{c}\prime}, \boldsymbol{W}^{\mathsf{d}\prime} \right)'$, where $\boldsymbol{W}^{\mathsf{c}}$ and $\boldsymbol{W}^{\mathsf{d}}$ are $q^{\mathsf{c}} \times 1$ and $q^{\mathsf{d}} \times 1$ subvectors, respectively.

New Measure and Statistic

This section introduces our measure of conditional dependence. We focus on the dependence between a scalar response variable $Y$ and a single candidate covariate $X$, conditional on $\boldsymbol{W}$.

Note that $X\perp Y|\boldsymbol{W}=\boldsymbol{w}$ holds if and only if the conditional distribution of $X$ given $(Y, \boldsymbol{W})$ is invariant with respect to the value of $Y$. That is,

for any values $y$ and $y'$. To quantify departures from this condition, we use the KS distance. Let $\Lambda_X(Y=y,\boldsymbol{W}=\boldsymbol{w})$ be the KS distance between the conditional CDF of $X$ and the same CDF evaluated at a user-specified reference value $y^*$: \[ \Lambda_{X}\left(Y=y,\boldsymbol{W}=\boldsymbol{w}\right)=\sup_{x}\left|F\left(x|y,\boldsymbol{w}\right)-F\left(x|y^{*},\boldsymbol{w}\right)\right|. \] Here, $y^*$ can be set to a fixed value, such as the sample median of $Y$. Clearly, $X\perp Y\mid \boldsymbol{W}=\boldsymbol{w}$ if and only if $\sup_{y}\Lambda_{X}(Y=y,\boldsymbol{W}=\boldsymbol{w})=0$.

It is worth noting that we do not need to use a statistic of the following form:

In fact, ((ref)) is strictly positive if and only if $\sup_{y}\sup_{x}\left|F\left(x|y,\boldsymbol{w}\right)-F\left(x|y^{*},\boldsymbol{w}\right)\right|>0$. To see this, suppose there exists a pair $(y,y')$ such that \[ \sup_{x}\left|F(x\mid y,\boldsymbol{w})-F(x\mid y',\boldsymbol{w})\right|\geq\delta>0. \] Then, by the triangle inequality, we must have \[ \max\left\{ \sup_{x}\left|F(x\mid y,\boldsymbol{w})-F(x\mid y^{*},\boldsymbol{w})\right|,\sup_{x}\left|F(x\mid y',\boldsymbol{w})-F(x\mid y^{*},\boldsymbol{w})\right|\right\} \geq\frac{\delta}{2}>0. \] The reverse direction is obvious. Using $\Lambda_{X}(Y=y,\boldsymbol{W}=\boldsymbol{w})=0$ (i.e., fixing $y'$ to $y^{*}$ in ((ref))) can substantially reduce the computational burden, as further discussed in Remark (ref).

Since $X\perp Y|\boldsymbol{W}$ requires the distributional equivalence in ((ref)) to hold for all values of $\left(y,\boldsymbol{w}\right)$, we define our overall (single summary) measure of dependence by taking the expectation of $\Lambda_X$ over the joint distribution of $(Y, \boldsymbol{W})$:

By excluding events with zero probability under the distribution of \((Y, \boldsymbol{W})\), we have that \(X \perp Y \mid \boldsymbol{W}\) if and only if \(\rho = 0\). To see this, note that

Thus, \(\rho > 0\) if and only if there exists a set of \(\boldsymbol{w}\) with nonzero measure such that, conditional on these \(\boldsymbol{w}\), there exists a set of \(y\) (also with nonzero measure) for which \(\Lambda_X(y, \boldsymbol{w}) > 0\).

The sample counterpart of $\Lambda_{X}(Y=y,\boldsymbol{W}=\boldsymbol{w})$ can be estimated as

where the estimated CDFs are given by

Here, $w_{l}^{\mathsf{c}}$ and $w_{l}^{\mathsf{d}}$ denote the $l$-th elements of the continuous and discrete parts of $\boldsymbol{w}$, respectively. $K_{h} ( \cdot ) = h^{-1} K ( \cdot )$, where $K ( \cdot )$ is a standard kernel function for continuous covariates. The kernel for discrete covariates is the Aitchison-Aitken kernel, defined as \[ K_{\lambda}^{\mathsf{d}} \left( w_{l}^{\mathsf{d}}, w^{\mathsf{d}} \right) = \left( \frac{\lambda}{r_{l} - 1} \right)^{ \mathbf{1} \left( w_{l}^{\mathsf{d}} \neq w^{\mathsf{d}} \right) } \left( 1 - \lambda \right)^{ \mathbf{1} \left( w_{l}^{\mathsf{d}} = w^{\mathsf{d}} \right) }, \] where $r_{l}$ is the number of distinct atoms in the support of $W_{l}^{\mathsf{d}}$ and $0 \leq \lambda \leq (r_{l} - 1)/r_{l}$. For simplicity, we assume the same bandwidth $h$ and tuning parameter $\lambda$ across all covariates. The use of the Aitchison-Aitken kernel is standard for handling discrete data; see, e.g., Chapter 4 in LiRacine.

Finally, our estimator of $\rho$ is defined as

In Appendix (ref), we develop a highly efficient algorithm to compute $\hat{\rho}$ at a computational cost proportional to $n^2$. We conclude this section with two remarks: one highlighting the computational improvement of our approach over a direct generalization of MaiZou2015, and another comparing our measure with other popular conditional dependence measures in the literature.

Screening and Its Error Bound

We adopt the conditional dependence measure developed in the previous section for our proposed screening procedure. For each candidate covariate $X_j$, where $j = 1, \dots, p$, we define \[ \rho_j = \mathbb{E} \left[ \Lambda_{X_j} \left( Y, \boldsymbol{W} \right) \right], \] as in ((ref)), with its sample counterpart given by ((ref)). We impose the following technical conditions.

Assumption (ref)(1) defines the sets of relevant and irrelevant covariates. The second part of Assumption (ref)(1) imposes a minimum signal strength condition, ensuring that the dependence of each relevant covariate on $Y$ is sufficiently large to be detected in finite samples. Assumption (ref)(2) allows $p$ to grow at any polynomial rate. This can be relaxed to an exponential rate at the cost of a slight sacrifice of power as reflected in the signal strength requirement in Assumption (ref)(1). Assumption (ref)(3) requires a random sample. Assumptions (ref)(4) and (5) guarantee that the denominators in $\hat{\rho}_j$ are uniformly bounded away from zero with high probability. Finally, Assumptions (ref)(6) and (7) are standard kernel regularity and smoothness conditions used to control the bias in nonparametric estimation.

We define the set containing the indices of all relevant variables as $\mathcal{M}^{*} = \left\{ 1,...,s^{*}\right\}$. Its corresponding estimator, $\hat{\mathcal{M}}$, is constructed by collecting the indices of all covariates whose estimated dependence measure, $\hat{\rho}_j$, exceeds a certain threshold $Cn^{-r/\left(2r+q^{\mathsf{c}}+1\right)}\log n$ for a fixed positive constant $C$, i.e.,

Theorem (ref) establishes that, with appropriately chosen thresholds in ((ref)), our screening procedure can select the set of relevant covariates with probability approaching one. This result can also be readily extended to the case where $Y$ is a binary random variable, in which the conditional density of $Y$ reduces to the propensity score. We discuss this extension in detail in Section (ref).

Thresholds, CV Refinement and Post-Selection Estimation

While Theorem (ref) guarantees that our screening procedure works in theory, the threshold defined in ((ref)) depends on an unknown constant $C$ that is difficult to choose in practice. To address this, we move beyond a fixed (subjective) threshold and propose a data-driven second stage to refine the set of selected variables. Our approach utilizes the CV method of HallRacineLi, which eliminates irrelevant covariates through kernel bandwidth selection. The refinement proceeds as follows: first, we use our screening method to select a reasonably large set of $\tilde{p}$ candidate covariates, including the pre-selected $\boldsymbol{W}$ and the top covariates in $\boldsymbol{X}$ ranked by $\hat{\rho}_j$. Then, we apply a CV-based refinement procedure to this set of $\tilde{p}$ variables to determine the final conditioning set. A data-driven method for choosing $\tilde{p}$ is provided in Section (ref).

One important advantage (and added value) of this refinement step, relative to using a fixed threshold in ((ref)), is its ability to further filter out certain “indirectly relevant” variables. Consider the following scenario: suppose there exists an $s^{**}<s^{*}$ such that

and such conditional independence is unambiguous.\footnote{Conditional independence can be ambiguous. For example, it is possible for $Y \perp X_1 \mid X_2$, $Y \perp X_2 \mid X_1$, and $Y \not\perp (X_1, X_2)$ to all hold simultaneously. The notion of “unambiguous” is formalized in Assumption (ref).} In ((ref)), although $X_{s^{**}+1},...,X_{s^{*}}$ may be dependent on $Y$ conditional on $\boldsymbol{W}$ and thus survive the screening step, they may become conditionally independent of $Y$ once we condition on the subset

In this case, the conditional density of $Y$ simplifies as \[ f\left(Y|X_{1},...,X_{s^{*}},\boldsymbol{W}\right) =\frac{f\left(Y,X_{s^{**}+1},...,X_{s^{*}}|\boldsymbol{Z}^{*}\right)}{f\left(X_{s^{**}+1},...,X_{s^{*}}|\boldsymbol{Z}^{*}\right)} = \frac{f\left(Y|\boldsymbol{Z}^{*}\right)f\left(X_{s^{**}+1},...,X_{s^{*}}|\boldsymbol{Z}^{*}\right)}{f\left(X_{s^{**}+1},...,X_{s^{*}}|\boldsymbol{Z}^{*}\right)} =f\left(Y|\boldsymbol{Z}^{*}\right). \] We refer to $X_{s^{**}+1},..., X_{s^{*}}$ as indirectly relevant covariates and to $\boldsymbol{Z}^{*}$ as directly relevant covariates.

To see that ((ref)) is a realistic concern, consider the following illustrative example: \[ Y=g\left(X_{1},\boldsymbol{W},e\right),\text{ }X_{1}=U^{\ast}+U_{1},\text{ }\left(X_{2},...,X_{s}\right)\perp\left(U_{1},\boldsymbol{W},e\right), \] for some $s > 1$, where $g$ is an unknown smooth function. In this setup, $Y\perp\left(X_{2},...,X_{s}\right)|\left(X_{1},\boldsymbol{W}\right)$ holds, and the conditional independence is unambiguous due to the presence of $U_{1}$. However, it is possible that $\left(X_{2},...,X_{s}\right)\not\perp Y|\boldsymbol{W}$ holds through $U^{\ast}$.

For the refinement stage, we retain $\tilde{p}-q^{\mathsf{c}}-q^{\mathsf{d}}$ (recall $\textrm{dim}(\boldsymbol{W})=q^{\mathsf{c}}+q^{\mathsf{d}}$) covariates from $\boldsymbol{X}$ with the largest values of $\hat{\rho}_j$. Together with $\boldsymbol{W}$, this yields a total of $\tilde{p}$ covariates in the conditioning set for us to work with. For the post-selection estimation and the determination of covariates relevance, we do not distinguish between covariates from $\boldsymbol{X}$ and those in $\boldsymbol{W}$ for simplicity.\footnote{In applications, practitioners can exclude certain or all elements of $\boldsymbol{W}$ from this refinement step, as illustrated in the simulations in Section (ref).} Without loss of generality, we denote the surviving covariates after screening (including $\boldsymbol{W}$) as \[ \boldsymbol{\tilde{X}}\equiv\left(X_{1}^{\mathsf{c}},...,X_{\tilde{s}_{1}}^{\mathsf{c}},X_{1}^{\mathsf{d}},...,X_{\tilde{s}_{2}}^{\mathsf{d}}\right)', \] where $\tilde{s}_1$ and $\tilde{s}_2$ denote the number of continuous and discrete covariates, respectively, in $\boldsymbol{\tilde{X}}$.

We then estimate the conditional density of $Y$ given $\boldsymbol{\tilde{X}}$ using the standard Nadaraya–Watson estimator:

where $r_l$ denotes the number of support points of the discrete covariate $X_l^{\mathsf{d}}$, $0\leq\lambda_l\leq\left.\left(r_{l}-1\right)\right/r_{l}$, the kernel functions $K_{h}(\cdot)$, $K_{h_{l}}(\cdot)$, and $K_{\lambda_{l}}^{\mathsf{d}}(\cdot)$ are defined as in ((ref)).

The logic of using CV for variable selection is based on how kernel bandwidths relate to variable relevance. It is easy to see that when the bandwidth $h_l$ for a continuous covariate $X_{l}^{\mathsf{c}}$ is much larger than the range of its support, the kernel $K_{h_{l}}(\cdot)$ assigns nearly the same weight to all observations, meaning $X_{l}^{\mathsf{c}}$ becomes almost uninformative to the estimate $\hat{f}(y|\boldsymbol{\tilde{x}})$. Similarly, when the bandwidth $\lambda_l$ for a discrete covariate $X_{l}^{\mathsf{d}}$ reaches its upper bound of $(r_l-1)/r_l$, $K_{\lambda_{l}}^{\mathsf{d}}\left(x_{l}^{\mathsf{d}},x^{\mathsf{d}}\right)=1/r_{l}$ for all $\left(x_{l}^{\mathsf{d}},x^{\mathsf{d}}\right)$, effectively removing the variable's influence on the density estimate.

In view of this observation, HallRacineLi find that searching for optimal bandwidths to minimize the cross-validated ISE (CV-ISE) will force the bandwidths for irrelevant covariates to these upper extremes. This procedure thus essentially eliminates (“smoothes out”) the influence of irrelevant variables, providing a data-driven method for variable selection. We briefly review the method of HallRacineLi in Appendix (ref).

Formally, the ISE is defined as

The optimal bandwidths are those that minimize the CV-ISE:

and the post-selection conditional density is then obtained using $(\hat{h},\hat{h}_{1},...\hat{h}_{\tilde{s}_{1}},\hat{\lambda}_{1},...,\hat{\lambda}_{\tilde{s}_{2}})$:

For notational simplicity, we abuse the notation to still use $\hat{f}$ to denote the post-selection estimator.

To ensure the validity of $\hat{f}(y|\boldsymbol{\tilde{x}})$ in ((ref)), we impose the following additional assumptions:

Assumption (ref)(1) defines the unambiguity of the conditional independence. $\left(Y,X_{1},...,X_{s^{*}}\right)\perp\left(X_{s^{*}+1},...,X_{p}\right)|\boldsymbol{W}$ is imposed to ensure that the ambiguity does not occur when irrelevant covariates are included. Assumption (ref)(2) is a standard regularity condition required as in HallRacineLi. Assumption (ref)(3) implicitly requires that $s^{*}$ is fixed, so that nonparametric estimation remains feasible. This assumption also guarantees that the initial screening does not mistakenly discard any relevant covariates due to choosing a too small $\tilde{p}$.

The results in HallRacineLi rely on the pure (unconditional) independence to avoid ambiguity. In the following proposition, we (slightly) generalize the results of HallRacineLi to the case in ((ref)) when the conditional independence is not ambiguous.

Proposition (ref) shows that the presence of irrelevant or indirectly relevant variables has no impact on the convergence rate of $\hat{f}\left(y|\boldsymbol{\tilde{x}}\right)$. The post-selection estimator attains the fastest convergence rate achievable as if the estimation were performed using only the conditioning set $\boldsymbol{Z}^{*}$.

A Computationally Efficient Modified CV Procedure

In practice, this procedure tends to be quite time-consuming, and the process can be numerically unstable due to the highly nonlinear nature of the ISE criterion function. Finding the optimal set of bandwidths that minimizes the ISE can be a challenging optimization problem.

We find that an exhaustive search over the full bandwidth space is unnecessary. Instead, we can achieve the same result by evaluating the ISE only at the theoretically optimal bandwidths corresponding to each possible subset (or “selection”) of the $\tilde{p}$ candidate covariates and then taking the one that minimizes the ISE. Although conceptually simple, we need to introduce additional notation to formalize this procedure.

For the $\tilde{p}$ covariates in $\boldsymbol{\tilde{X}}$, there are $2^{\tilde{p}}$ possible selections. Each selection can be represented by a $\tilde{p}\times1$ binary selection vector $\boldsymbol{I}$, with 1 indicating inclusion and 0 indicating exclusion of the corresponding covariate. For example, \[ \boldsymbol{I}=(\underset{\tilde{p}}{\underbrace{0,0,\dots,0}})' \] indicates that none of the covariates are selected, while \[ \boldsymbol{I}^{*} \equiv \left( \underbrace{1, \dots, 1}_{p^{\mathsf{c}*}}, \underbrace{0, \dots, 0}_{\tilde{s}_1 - p^{\mathsf{c}*}}, \underbrace{1, \dots, 1}_{p^{\mathsf{d}*}}, \underbrace{0, \dots, 0}_{\tilde{s}_2 - p^{\mathsf{d}*}} \right)' \] corresponds to the correct selection of all directly relevant covariates $\boldsymbol{Z}^{*}$, where $p^{\mathsf{c}*}$ and $p^{\mathsf{d}*}$ denote the numbers of continuous and discrete variables in $\boldsymbol{Z}^{*}$, respectively. We collect all such candidate selection vectors into the set

We define $\boldsymbol{h}_{\boldsymbol{I}}^{\textrm{opt}}\in\mathbb{R}^{\tilde{p}+1}$ as the vector of theoretically optimal bandwidths corresponding to the selection $\boldsymbol{I} \in \mathcal{S}$. For example, when no covariates are selected, i.e., $\boldsymbol{I}=(0,0,\dots,0)'$, we have \[ \boldsymbol{h}_{\boldsymbol{I}}^{\textrm{opt}}=\left(c_{0}n^{\frac{-1}{1+2r}},\underset{\tilde{s}_{1}}{\underbrace{\infty,\dots,\infty}},\underset{\tilde{s}_{2}}{\underbrace{\frac{r_{1}-1}{r_{1}},\frac{r_{2}-1}{r_{2}},\dots,\frac{r_{\tilde{s}_{2}}-1}{r_{\tilde{s}_{2}}}}}\right)', \] where $c_{0}n^{\frac{-1}{1+2r}}$ is the optimal bandwidth for $Y$ with no conditioning covariates, and the bandwidths for all covariates reach their upper limit. In practice, “$\infty$” for a continuous covariate means that it is assigned uniform weight across all observations, which is analogous to using $(r_{l}-1)/r_{l}$ for a discrete covariate $X_{l}^{\mathsf{d}}$. When only $\boldsymbol{Z}^{*}$ is selected, corresponding to $\boldsymbol{I}^{*}$, we obtain \[ \boldsymbol{h}_{\boldsymbol{I}^{*}}^{\textrm{opt}}=\left(h_{0}^{*},h_{1}^{*},\dots,h_{p^{\mathsf{c}*}}^{*},\infty,\dots,\infty,\lambda_{1}^{*},\dots,\lambda_{p^{\mathsf{d}*}}^{*},\frac{r_{p^{\mathsf{d}*}+1}-1}{r_{p^{\mathsf{d}*}+1}},\dots,\frac{r_{\tilde{s}_{2}}-1}{r_{\tilde{s}_{2}}}\right)', \] where $h_{0}^{*}$, $h_{l}^{*}$, and $\lambda_{l}^{*}$ are the optimal bandwidths as given in Proposition (ref).

The following corollary provides the key theoretical justification for our modified CV procedure. Since the corollary follows directly from Proposition (ref), the proof is omitted. It states that the ISE evaluated at the optimal bandwidths, $\boldsymbol{h}_{\boldsymbol{I}^{*}}^{\textrm{opt}}$, is, with high probability, the minimum among the ISEs evaluated for all other $\boldsymbol{h}_{\boldsymbol{I}}^{\textrm{opt}}$. This result is crucial as it reduces the complex continuous optimization problem to a much faster and more stable discrete search over the $2^{\tilde{p}}$ selection vectors in $\mathcal{S}$.

In practice, we apply the proposed discrete search algorithm to the CV-ISE to obtain an estimate of the optimal selection vector, $\hat{\boldsymbol{I}}$, and its corresponding bandwidths, $\boldsymbol{h}_{\hat{\boldsymbol{I}}}^{\textrm{opt}}$.

Propensity Score Estimation

The propensity score plays an important role in many statistical and econometric methods. This section details how our procedure can be applied to its estimation. We denote the binary response (treatment) variable as $D \in \{0, 1\}$. To identify covariates relevant to the propensity score, we adapt the dependence measure from Section (ref). Specifically, we modify our KS-based distance measure to compare the distribution of a candidate covariate $X$ for the treated group ($D=1$) with that of the control group ($D=0$), conditional on $\boldsymbol{W}$: \[ \Lambda_{X}\left(\boldsymbol{W}=\boldsymbol{w}\right) =\sup_{x}\left|F\left(x|D=1,\boldsymbol{w}\right)-F\left(x|D=0,\boldsymbol{w}\right)\right| \equiv\sup_{x}\left|F_{1}\left(x|\boldsymbol{w}\right)-F_{0}\left(x|\boldsymbol{w}\right)\right|. \] The overall dependence measure and its sample counterpart are then defined analogously to the general case: \[ \rho=\mathbb{E}\left[\Lambda_{X}\left(\boldsymbol{W}\right)\right] \quad \text{and} \quad \hat{\rho}=\frac{1}{n}\sum_{i=1}^{n}\hat{\Lambda}_{X}\left(\boldsymbol{W}=\boldsymbol{w}_{i}\right), \] where \[ \hat{\Lambda}_{X}\left(\boldsymbol{W}=\boldsymbol{w}\right)=\max_{i=1,...,n}\left|\hat{F}_{1,n}\left(x_{i}|\boldsymbol{w}\right)-\hat{F}_{0,n}\left(x_{i}|\boldsymbol{w}\right)\right| \] with \[ \hat{F}_{d,n}\left(x|\boldsymbol{w}\right)=\frac{\sum_{i=1}^{n}\mathbf{1}\left(x_{i}\leq x\right)\mathbf{1}\left(D_{i}=d\right)\left[\Pi_{l=1}^{q^{\mathsf{c}}}K_{h}\left(w_{li}^{\mathsf{c}}-w_{l}^{\mathsf{c}}\right)\right]\left[\Pi_{l=1}^{q^{\mathsf{d}}}K_{\lambda}^{\mathsf{d}}\left(w_{li}^{\mathsf{d}},w_{l}^{\mathsf{d}}\right)\right]}{\sum_{i=1}^{n}\mathbf{1}\left(D_{i}=d\right)\left[\Pi_{l=1}^{q^{\mathsf{c}}}K_{h}\left(w_{li}^{\mathsf{c}}-w_{l}^{\mathsf{c}}\right)\right]\left[\Pi_{l=1}^{q^{\mathsf{d}}}K_{\lambda}^{\mathsf{d}}\left(w_{li}^{\mathsf{d}},w_{l}^{\mathsf{d}}\right)\right]},\textrm{ }d=0,1. \]

One notable difference here from the general case is that because $D$ is binary, no kernel smoothing is required for it.

After screening, the post-selection estimate $\hat{f} \left( D = d \mid \tilde{\boldsymbol{x}} \right)$ can be obtained via the refinement procedure outlined in Section (ref). The final estimator is given by

where $(\hat{h}_{1},...,\hat{h}_{\tilde{s}_{1}},\hat{\lambda}_{1},...,\hat{\lambda}_{\tilde{s}_{2}})$ are the optimal bandwidths that minimize the ISE in the CV-refinement step.

The proofs of Corollaries (ref) and (ref) are straightforward applications of Theorem (ref) and Proposition (ref) and are thus omitted.

Summary of the Proposed Procedure

For ease of reference, we summarize the detailed steps of our procedure proposed in Sections (ref)--(ref) below. Recall that the dimension of the pre-selected covariates $\boldsymbol{W}$ is $q^{\mathsf{c}}+q^{\mathsf{d}}$. We assume that $q^{\mathsf{c}}+q^{\mathsf{d}}\leq3$ without loss of generality.

It is important to note the limitation of this approach: its performance relies on the sparsity assumption, i.e., the number of truly relevant covariates is fixed. If the procedure continues to iterate until the final selected $\tilde{p}$ is too large for reliable nonparametric estimation (e.g., $\tilde{p}=10$ for $n=1000$), our method may not be suitable.

Estimating Average Treatment Effects

In this section, we apply the methods proposed in Section (ref) to estimate and conduct inference on ATE in high-dimensional settings. We use this context to demonstrate how our screening procedure and dimension-reduced density estimator, combined with ML methods, enhance the applicability and performance of established approaches. The identification results presented here may also be of independent interest. Given the extensive literature on ATE, we provide only a brief overview of the framework in Section (ref). For a comprehensive review of related methods, readers may refer to ImbensRubin and ding2024first. Additionally, chernozhukov2024applied discuss recent advances in applying ML and AI techniques to causal inference. It is worth emphasizing that the applications of our proposed methods are not limited to ATE estimation; they are broadly applicable and can complement many of the procedures discussed in the aforementioned works and the references therein.

Reviewing Average Treatment Effects and Doubly Robust Estimator

In this section, we review the potential outcomes framework and the doubly robust estimator for ATE, providing the necessary background for the subsequent discussion.

Suppose we observe a random sample of $n$ individuals indexed by $i=1,...,n$. We consider the classic potential outcomes framework for observational studies under the stable unit treatment value assumption (rubin1980randomization). Let $Y$ denote an outcome variable and $D$ a binary treatment indicator---$D=0$ for the control group and $D=1$ for the treatment group. Define $Y(0)$ and $Y(1)$ as potential outcomes corresponding to $D=0$ and $D=1$, respectively. The observed outcome $Y$ is determined by $Y=Y(0)\left(1-D\right)+Y(1)D$. In addition to $(Y,D)$, we observe a vector of pre-treatment covariates (control variables) $\boldsymbol{Z}\equiv\left(\boldsymbol{X}',\boldsymbol{W}'\right)'$, where we use this notation to align with the setting in Section (ref). We denote the conditional means of the potential outcomes as

The population ATE to estimate is then defined as \[ \psi=\mathbb{E}[g_{1}\left(\boldsymbol{Z}\right)-g_{0}\left(\boldsymbol{Z}\right)]. \]

Let $m\left(\boldsymbol{Z}\right)\equiv\mathbb{E}\left[D|\boldsymbol{Z}\right]=\Pr(D=1|\boldsymbol{Z})$, commonly referred to as the propensity score. The following assumptions are standard in the literature. Assumption (ref)(1) is known as the ignorability condition (rubin1978bayesian and rosenbaum1983central), while Assumption (ref)(2) is referred to as the overlap condition, which ensures the estimability of the ATE (KhanTamer).\footnote{The ignorability assumption is also known as the unconfoundedness or selection on observables assumption.}

Under Assumption (ref), we have $g_{0}(\boldsymbol{Z})=\mathbb{E}[Y|D=0,\boldsymbol{Z}]$, $g_{1}(\boldsymbol{Z})=\mathbb{E}[Y|D=1,\boldsymbol{Z}]$, and the following identification formula for $\psi$:

Expression ((ref)) motivates the following moment estimator for $\psi$:

where $(\hat{g}_{0},\hat{g}_{1})$ and $\hat{m}$ are generic estimators for $(g_{0},g_{1})$ and $m$, respectively. The estimator $\hat{\psi}$ is doubly robust (scharfstein1999adjusting and bang2005doubly) because it is consistent if either $(\hat{g}_{0},\hat{g}_{1})$ are consistent estimators for $(g_{0},g_{1})$ or $\hat{m}$ is a consistent estimator for $m$. Similar doubly robust estimators can be constructed for other treatment effects, such as the average treatment effect on the treated (ATT) and on the control (ATC).

Doubly Robust and De-Biased ATE Estimators with Dimension-Reduced Propensity Score

Assumption (ref)(1) (ignorability) plays a central role in observational studies. It assumes the absence of unobserved confounders---covariates that affect both the treatment and the outcome. However, this assumption is generally not testable and is often justified based on the researcher's domain knowledge. To minimize the risk of violating this assumption, practitioners may control for as many pre-treatment covariates $\boldsymbol{Z}$ as possible (rubin2007design and vander2011new). Nevertheless, there is an intrinsic trade-off between the ignorability assumption and Assumption (ref)(2) (overlap) (d2021overlap). Controlling for a rich set of covariates $\boldsymbol{Z}$ makes the former more plausible but simultaneously reduces the randomness in the propensity score model, making it harder to satisfy the latter as $D$ becomes more predictable. In applications, even if the overlap assumption holds, the estimated propensity scores $\hat{m}(\boldsymbol{z})$ can still be close to 0 or 1, particularly when overfitting occurs due to the inclusion of many covariates irrelevant to $D$. This can lead to poor finite-sample performance of estimators involving inverse propensity weighting, such as ((ref)). Therefore, it is practically important to refine the conditioning set $\boldsymbol{Z}$ to ensure the ignorability condition while mitigating its tension with the overlap condition.\footnote{Adjusting for many covariates also increases the risk of over-adjustment problems (ding2015adjust). Since this paper focuses on estimation and inference, we do not discuss their implications for the identification of the ATE.}

Furthermore, the doubly robust estimator ((ref)) can become inconsistent (or even “doubly fragile”) when both the outcome regression model $(g_{0},g_{1})$ and the propensity score model $m$ are misspecified (KangSchafer2007). To address this issue, recent advancements in DML methods, such as farrell2015robust, ChernozhukovEtal2018, and FarrellEtal2021, propose employing (nonparametric) ML tools to estimate $(g_{0},g_{1})$ and $m$. These approaches enhance the robustness of the estimator ((ref)) by avoiding strong parametric restrictions, such as logistic or linear regressions. However, the drawback is that the estimator may suffer from the curse of dimensionality when $\textrm{dim}(\boldsymbol{Z})$ is large or even diverges as the sample size increases. This challenge highlights another practical importance of refining the conditioning set $\boldsymbol{Z}$--dimension reduction for enabling robust DML estimation.

In what follows, we demonstrate how our proposed methods facilitate such refinement under certain sparsity conditions. The proof of the following result is provided in Appendix (ref).

Proposition (ref) ensures that under the additional sparsity restriction $D\perp\boldsymbol{\underline{Z}}|\boldsymbol{Z}^{*}$, the ignorability condition remains valid when conditioning only on $\boldsymbol{Z}^{*}$. Consequently, the ATE $\psi$ can be consistently estimated using ((ref)) or ((ref)) with the lower-dimensional subvector $\boldsymbol{Z}^{*}$ instead of the full covariate set $\boldsymbol{Z}$, achieving a “dual dimension reduction” for both propensity score and outcome regression models. Moreover, the overlap condition required for these estimators is weaker than Assumption (ref)(2), expanding their applicability when the overlap condition on $m(\boldsymbol{Z})$ fails but holds for $m(\boldsymbol{Z}^{*})$.

This result is particularly relevant when the sparsity condition in Assumptions (ref) and (ref) holds. Specifically, we can apply our approach proposed in Section (ref) to identify $\boldsymbol{Z}^{*}$ with high probability, as ensured by Theorem (ref) (Corollary (ref)). Since $\text{dim}(\boldsymbol{Z}^{*})$ is fixed and much smaller than $\text{dim}(\boldsymbol{Z})$, the overlap condition on $\boldsymbol{Z}^{*}$ is both weaker and more plausible than that on $\boldsymbol{Z}$.

Furthermore, there exists a complementary result, symmetric to Proposition (ref), that relies on a similar sparsity condition on $Y \mid D$. This result is presented in Appendix (ref). Taken together, these results illustrate that our proposed approach can exploit a “double sparsity” structure: when either $D$ or $Y|D$ exhibits a sparsity pattern similar to that in Proposition (ref), our method can recover the corresponding relevant covariates and achieve effective dimension reduction. We provide a more detailed discussion of this point in Appendix (ref).

To illustrate the usefulness of Proposition (ref) in robust causal inference, we replicate the well-established results of DML developed in farrell2015robust, ChernozhukovEtal2018, and FarrellEtal2021. We can show that the ATE estimator ((ref)), when using identified relevant covariates, has the following asymptotic properties:

Theorem (ref) states that under standard conditions, $\hat{\psi}$ has the same asymptotic distribution as it would if $(g_{0},g_{1})$ and $m$ were known. This result follows directly from the existing literature and Corollary (ref), which establishes that our screening procedure can recover $\boldsymbol{Z}^{*}$ with high probability. Thus, we omit the proof. The key implication is that our procedure can achieve dimension reduction not only for $D$ and $m$, but also for $Y\left(d\right)$ and $g_{d}$, $d=0,1$.

We briefly verify conditions (1)--(3) in Theorem (ref). For conditions (1) and (2), we assume that $(\hat{g}_{0},\hat{g}_{1})$ are certain estimators from the literature that satisfy $\{ n^{-1}\sum_{i=1}^{n}\left[\hat{g}_{d}\left(\boldsymbol{z}_{i}^{*}\right)-g_{d}\left(\boldsymbol{z}_{i}^{*}\right)\right]^{2}\} {}^{1/2}=O_{P}(n^{-1/4})$ for $d=0,1$, under standard regularity conditions. We show in Lemma (ref) that $\hat{m}\left(\boldsymbol{z}_{i}^{\hat{*}}\right)$ also satisfies the required convergence rate, provided that $r/\left(p^{\mathsf{c}\ast}+2r\right)>1/4$, aligning with results in FarrellEtal2021. Condition (3) can be satisfied through a “sample splitting” procedure, as described in ChernozhukovEtal2018.

We end this section with the following remarks.

Simulations

This section assesses the finite-sample performance of the procedure proposed in Section (ref) through Monte Carlo experiments. All simulations presented in this section are conducted using the R programming language, based on 1,000 independent replications. In Section (ref), we examine dimension reduction in conditional density estimation with a continuous response variable \( Y \). In Section (ref), we study the performance of our procedure in estimating propensity scores for the estimation and inference of ATE.

Simulations with Continuous Response

We explore three distinct simulation designs, each differing in how the response variable \( Y \) depends on covariates $\boldsymbol{X}$, based on the following linear varying coefficient model:

Throughout this section, we assume the error term \( \epsilon \sim N(0,1) \) and consider the total number of covariates \( p \in \{20, 50, 100\} \), with \( \boldsymbol{X}^{\mathsf{c}} \) and \( \boldsymbol{X}^{\mathsf{d}} \) each having dimension \( p/2 \). We assume there is no \( \boldsymbol{W}^{\mathsf{d}} \) and only a scalar \( W^{\mathsf{c}} \), i.e., \( \boldsymbol{W}= W^{\mathsf{c}}\). The true parameter functions for the continuous and discrete covariates are specified as follows: \[ \underset{(p/2) \times 1}{\boldsymbol{\beta}^{\mathsf{c}}(W^{\mathsf{c}})} = \left( 0.4 W^{\mathsf{c}} + 0.5, \sin(2\pi W^{\mathsf{c}}) + 0.5, 0, \dots, 0 \right)', \] and \[ \underset{(p/2) \times 1}{\boldsymbol{\beta}^{\mathsf{d}}(W^{\mathsf{c}})} = \left( 0.5 W^{\mathsf{c}} + 0.7, 0, \dots, 0 \right)'. \] Thus, \(\boldsymbol{X}(\mathcal{M}^{*}) \equiv (X_{1}^{\mathsf{c}}, X_{2}^{\mathsf{c}}, X_{1}^{\mathsf{d}})'\) are directly relevant to \( Y \), while the remaining covariates are not.

\noindentDesign 1: We generate a \(\left(p+1\right) \times 1\) vector \(\boldsymbol{e} = (e_1, \dots, e_{p+1})'\) independently from \( N(0, \boldsymbol{I}_{p+1}) \), where \(\boldsymbol{I}_{p+1}\) denotes the \(\left(p+1\right) \times \left(p+1\right)\) identity matrix. We set \( W^{\mathsf{c}} = e_{p+1} \) and, for \( j = 1, \dots, p/2 \), define $X_{j}^{\mathsf{c}} = e_{j}$, and \[ X_{j}^{\mathsf{d}} = (-1) \cdot \boldsymbol{1} \left( \Phi(e_{j+p/2}) \leq 1/3 \right) + \boldsymbol{1} \left( \Phi(e_{j+p/2})>2/3 \right), \] where \(\Phi(\cdot)\) denotes the CDF of the standard normal distribution. Thus, in this design, \[ \left(Y, W^{\mathsf{c}}, \boldsymbol{X}(\mathcal{M}^{*})\right) \perp \left(X_{3}^{\mathsf{c}}, \dots, X_{p/2}^{\mathsf{c}}, X_{2}^{\mathsf{d}}, \dots, X_{p/2}^{\mathsf{d}}\right). \]

\noindentDesign 2: Let \(\boldsymbol{\Sigma}_{k}\) denote a \(k \times k\) covariance matrix with all diagonal elements equal to 1 and all off-diagonal elements equal to 0.25. We generate a \((p-2) \times 1\) vector \(\boldsymbol{e}_{1} = (e_{1,1}, \dots, e_{1,p-2})' \sim N(0, \boldsymbol{\Sigma}_{p-2})\) and a \(3 \times 1\) vector \(\boldsymbol{e}_{2} = (e_{2,1}, e_{2,2}, e_{2,3})' \sim N(0, \boldsymbol{\Sigma}_{3})\) independently. We set \( W^{\mathsf{c}} = e_{1,p-2} \), $X_{j+2}^{\mathsf{c}} = e_{1,j}$, for $j = 1, \dots, p/2 - 2$, and \[ X_{j+1}^{\mathsf{d}} = (-1) \cdot \boldsymbol{1} \left( \Phi(e_{1,j+p/2-2}) \leq 1/3 \right) + \boldsymbol{1} \left( \Phi(e_{1,j+p/2-2}) > 2/3 \right), \] for \( j = 1, \dots, p/2 - 1 \). For the relevant covariates, we define $X_{j}^{\mathsf{c}} = e_{2,j}$, for $j = 1,2$, and \[ X_{1}^{\mathsf{d}} = (-1) \cdot \boldsymbol{1} \left( \Phi(e_{2,3}) \leq 1/3 \right) + \boldsymbol{1} \left( \Phi(e_{2,3}) > 2/3 \right). \] Thus, in this design, \[ Y \not\perp \left(X_{3}^{\mathsf{c}}, \dots, X_{p/2}^{\mathsf{c}}, X_{2}^{\mathsf{d}}, \dots, X_{p/2}^{\mathsf{d}}\right), \textrm{ but } Y \perp \left.\left(X_{3}^{\mathsf{c}}, \dots, X_{p/2}^{\mathsf{c}}, X_{2}^{\mathsf{d}}, \dots, X_{p/2}^{\mathsf{d}}\right) \right| W^{\mathsf{c}}. \]

\noindentDesign 3: We generate a \((p+1) \times 1\) vector \(\boldsymbol{e} = (e_1, \dots, e_{p+1})' \sim N(0, \boldsymbol{\Sigma}_{p+1})\) and define \(\boldsymbol{X}^{\mathsf{c}}\), \(\boldsymbol{X}^{\mathsf{d}}\), and \(W^{\mathsf{c}}\) following the same formulae based on \(\boldsymbol{e}\) as in Design 1. Thus, in this design, \[ Y \not\perp \left.\left(X_{3}^{\mathsf{c}}, \dots, X_{p/2}^{\mathsf{c}}, X_{2}^{\mathsf{d}}, \dots, X_{p/2}^{\mathsf{d}}\right)\right| W^{\mathsf{c}}, \textrm{ but } Y \perp \left.\left(X_{3}^{\mathsf{c}}, \dots, X_{p/2}^{\mathsf{c}}, X_{2}^{\mathsf{d}}, \dots, X_{p/2}^{\mathsf{d}}\right)\right| \left(W^{\mathsf{c}}, \boldsymbol{X}(\mathcal{M}^{*})\right). \]

Design 1 is the benchmark design. Design 2 examines the effectiveness of our approach in handling the dependence between \( Y \) and \( (X_{3}^{\mathsf{c}}, \dots, X_{p/2}^{\mathsf{c}}, X_{2}^{\mathsf{d}}, \dots, X_{p/2}^{\mathsf{d}} )'\) through their shared dependence on \(\boldsymbol{W}\). Design 3 assesses the performance of our procedure under a more general dependence structure, as described in ((ref)), where the dependence between \( Y \) and \(( X_{3}^{\mathsf{c}}, \dots, X_{p/2}^{\mathsf{c}}, X_{2}^{\mathsf{d}}, \dots, X_{p/2}^{\mathsf{d}} )'\) arises from two sources: their common dependence on \(\boldsymbol{W}\) and their relationship with the set \(\boldsymbol{X}(\mathcal{M}^{*})\) of covariates directly relevant to \( Y \).

For each simulation design, we implement the dimension reduction procedure described in Section (ref) with sample sizes \( n \in \{250, 500, 1000, 2000\} \) and set \(\tilde{p} = 5\) to identify the set $\boldsymbol{X}(\mathcal{M}^{*})$ in the data-generating process (DGP) ((ref)). Let \( \boldsymbol{X}(\hat{\mathcal{M}}^{(1)}) \) denote the set of covariates selected after the first two steps of our procedure (screening), and \(\boldsymbol{X}(\hat{\mathcal{M}}^{(2)}) \) denote the set of covariates not assigned extreme bandwidths in the CV refining steps 3 and 4. In Appendix (ref), we report two sets of simulation results:

The two variable screening approaches based on our proposed conditional dependence measure can identify relevant covariates with high probability across Designs 1--3, even in high-dimensional settings, when $n \geq 500$. The “\(\textrm{Quantile-}\rho\)” approach, which uses the statistic defined in ((ref)) with multiple reference values, generally outperforms the “$\rho$” approach that relies solely on the sample median, particularly when screening discrete covariates in small samples. Supplementary simulations suggest that incorporating additional quantiles, such as the $(0.1, \ldots, 0.9)$ quantiles of $Y$, can further enhance the small-sample performance of the “\(\textrm{Quantile-}\rho\)” method.

In small samples, our approaches exhibit substantially higher CRRs in variable screening than FOCI and slightly lower CRRs than CDCSIS. As noted in Remark (ref), FOCI is the most computationally efficient, with a computational cost of $O(n\log n)$, but our simulation results indicate it requires relatively large sample sizes to attain a desirable CRR. Applying CDCSIS, which has a computational cost of $O(n^3)$, to simulations with $n\geq 2000$ can be quite time-consuming. Thus, we omit the corresponding results in Table (ref)--(ref), given that CDCSIS already achieved 100% CRR with $n=1000$. Our two approaches, however, have a computational cost of $O(n^2)$, and we observe that when $n\geq 1000$, they have almost identical CRR to CDCSIS. In other words, our approaches offer significantly more computationally efficient alternatives to CDCSIS for moderately large datasets.

Table (ref) reports results from the CV refinement step, showing that our “Modified-CV” algorithm achieves exact recovery with probability approaching one. It consistently outperforms the original “CV” procedure of HallRacineLi across all designs and sample sizes. As discussed in Section (ref), the “CV” procedure minimizes a CV-ISE criterion, which can be numerically unstable due to its high nonlinearity and multiple local minima. Our simulations demonstrate that the “Modified-CV” algorithm offers a more computationally efficient and stable post-screening refinement strategy. We omit the results of the original CV for $n=2000$ due to its high computation cost.

Simulations with Discrete Response

This section examines the performance of our procedure in propensity score estimation and its application to dimension-reduced doubly robust estimation of ATE, as proposed in Section (ref).

We consider simulation designs adapted from FarrellEtal2021. Specifically, we assume a binary treatment $D$ with the propensity score given by the varying coefficient logistic model:

where $\boldsymbol{Z}=(\boldsymbol{X}',\boldsymbol{W}')'$ and \(\mathcal{S}(\cdot)\) denotes the standardization operator such that \(\mathcal{S}(X) = [X -\mathbb{E}(X)]/\text{sd}(X)\) for any random variable $X$. As in Section (ref), the dimension of \( \boldsymbol{X} \) is set to \( p \in \{20, 50, 100\} \), where \( \boldsymbol{X}^{\mathsf{c}} \) and \( \boldsymbol{X}^{\mathsf{d}} \) each have dimension \( p/2 \), and \( \boldsymbol{W}=W^{\mathsf{c}} \). The true parameter functions in ((ref)) are specified as: \[ \boldsymbol{\alpha}^{\mathsf{c}}(W^{\mathsf{c}}) = \left( 0.4 W^{\mathsf{c}} + 0.5, \sin\left(2\pi W^{\mathsf{c}}\right)+0.5, 0, \dots, 0 \right)', \] and \[ \boldsymbol{\alpha}^{\mathsf{d}}(W^{\mathsf{c}}) = \left( 0.5W^{\mathsf{c}} + 0.7, 0, \dots, 0 \right)'. \] Thus, only \( \boldsymbol{X}(\mathcal{M}^{*}) = (X_{1}^{\mathsf{c}}, X_{2}^{\mathsf{c}}, X_{1}^{\mathsf{d}})'\) are directly relevant to \( D \), and $\boldsymbol{Z}^*=(\boldsymbol{X}(\mathcal{M}^{*})',W^{\mathsf{c}})'$.

We adopt the same DGPs for \( \boldsymbol{X}\) and \( W^{\mathsf{c}}\) as in Designs 1–3. These DGPs, when applied to the propensity score model specified in equation ((ref)), establish Designs 4–6, which feature the following dependence structures: \noindentDesign 4: This serves as our benchmark design, where \[ \left(D, W^{\mathsf{c}}, \boldsymbol{X}(\mathcal{M}^{*})\right) \perp \left(X_{3}^{\mathsf{c}}, \dots, X_{p/2}^{\mathsf{c}}, X_{2}^{\mathsf{d}}, \dots, X_{p/2}^{\mathsf{d}}\right). \] \noindentDesign 5: \( D \) and \( (X_{3}^{\mathsf{c}}, \dots, X_{p/2}^{\mathsf{c}}, X_{2}^{\mathsf{d}}, \dots, X_{p/2}^{\mathsf{d}} )'\) are dependent solely through their shared dependence on \( W^{\mathsf{c}} \), i.e., \[ D \not\perp \left(X_{3}^{\mathsf{c}}, \dots, X_{p/2}^{\mathsf{c}}, X_{2}^{\mathsf{d}}, \dots, X_{p/2}^{\mathsf{d}}\right), \text{ but } D \perp \left.\left(X_{3}^{\mathsf{c}}, \dots, X_{p/2}^{\mathsf{c}}, X_{2}^{\mathsf{d}}, \dots, X_{p/2}^{\mathsf{d}}\right)\right| W^{\mathsf{c}}. \] \noindentDesign 6: \( D \) depends on \( (X_{3}^{\mathsf{c}}, \dots, X_{p/2}^{\mathsf{c}}, X_{2}^{\mathsf{d}}, \dots, X_{p/2}^{\mathsf{d}} )'\) not only through \( W^{\mathsf{c}} \) but also through the covariates \(\boldsymbol{X}(\mathcal{M}^{*})\), which are directly relevant to \( D \), i.e., \[ D \not\perp \left.\left(X_{3}^{\mathsf{c}}, \dots, X_{p/2}^{\mathsf{c}}, X_{2}^{\mathsf{d}}, \dots, X_{p/2}^{\mathsf{d}}\right)\right| W^{\mathsf{c}}, \text{ but } D \perp \left.\left(X_{3}^{\mathsf{c}}, \dots, X_{p/2}^{\mathsf{c}}, X_{2}^{\mathsf{d}}, \dots, X_{p/2}^{\mathsf{d}}\right)\right| \left(W^{\mathsf{c}}, \boldsymbol{X}(\mathcal{M}^{*})\right). \]

Given the covariates and treatment variable, we specify the following DGP for the outcome \( Y \):

where $\boldsymbol{X}(\mathcal{G}^{*})\equiv (X_{1}^{\mathsf{c}}, X_{3}^{\mathsf{c}}, X_{1}^{\mathsf{d}})'$. Notably, the subsets of covariates directly relevant to the outcome and the treatment differ, with $\boldsymbol{X}(\mathcal{G}^{*})= (X_{1}^{\mathsf{c}}, X_{3}^{\mathsf{c}}, X_{1}^{\mathsf{d}})'$ and $\boldsymbol{X}(\mathcal{M}^{*}) = (X_{1}^{\mathsf{c}}, X_{2}^{\mathsf{c}}, X_{1}^{\mathsf{d}})'$. The untreated potential outcome \( g_{0}(\boldsymbol{X}(\mathcal{G}^{*}), W^{\mathsf{c}})\) and the conditional treatment effect \( \psi(\boldsymbol{X}(\mathcal{G}^{*}), W^{\mathsf{c}}) \) are defined as: \[ g_{0}(\boldsymbol{X}(\mathcal{G}^{*}), W^{\mathsf{c}}) = \varphi(\boldsymbol{X}(\mathcal{G}^{*}))'\beta_{g}(W^{\mathsf{c}}) \] and \[ \psi(\boldsymbol{X}(\mathcal{G}^{*}), W^{\mathsf{c}}) = \varphi(\boldsymbol{X}(\mathcal{G}^{*}))'\beta_{\psi}(W^{\mathsf{c}}), \] where \( \varphi(\boldsymbol{X}(\mathcal{G}^{*})) \) denotes the vector of all 9 second-degree polynomial terms of \( \boldsymbol{X}(\mathcal{G}^{*}) \), including all pairwise interactions. The coefficient functions are given by: \[ \beta_{g,k}(W^{\mathsf{c}}) = W^{\mathsf{c}} \cdot U_{g,k} \textrm{ with } U_{g,k} \sim N(0.3, 0.7) \] and \[ \beta_{\psi,k}(W^{\mathsf{c}}) = W^{\mathsf{c}} \cdot U_{\psi,k} \textrm{ with } U_{\psi,k} \sim U(0.1, 0.22), \] where the subscript \( k =1,...,9 \) indexes each component of \( \varphi(\boldsymbol{X}(\mathcal{G}^{*})) \). We integrate the outcome regression model (ref) into Designs 4–6 to estimate the ATE, defined as \( \psi = \mathbb{E}[\psi(\boldsymbol{X}(\mathcal{G}^{*}), W^{\mathsf{c}})] \).

We consider sample sizes $n\in\{500, 1000, 2000, 4000\}$ and split each sample into two folds $(I_1,I_2)$, of equal size. Note that the sample splitting is required for condition (3) in Theorem (ref). For implementation, we first apply the procedure proposed in Section (ref) with $\tilde{p}=5$ to sample $I_1$ (with sample size $n_1\in\{250, 500, 1000, 2000\}$) to obtain the dimension-reduced estimator $\hat{m}(\boldsymbol{Z}^{\hat{*}})$ of the propensity score \( m(\boldsymbol{Z})\), where \(\boldsymbol{Z}^{\hat{*}} \equiv (\boldsymbol{X}(\hat{\mathcal{M}}^{(2)}), W^{\mathsf{c}}) \), following the notation in Sections (ref) and (ref). This step aims to identify the components of \( \boldsymbol{X} \) that are directly relevant to the treatment \( D \), namely \( \boldsymbol{X}(\mathcal{M}^{*})\), with high probability. It is worth noting that, by Proposition (ref), consistent estimation of $\psi$ does not require recovering $\boldsymbol{X}(\mathcal{G}^{*})$.

Tables (ref)--(ref) and Table (ref) in Appendix (ref) present the dimension reduction results for Designs 4--6. These tables parallel Tables (ref)--(ref) and Table (ref) from Designs 1--3, respectively, reporting the same set of metrics. The only exception is the absence of the “$\textrm{Quantile-}\rho$” results, as $D$ is a binary variable in Designs 4--6. In these designs, our proposed screening method and CDCSIS demonstrate finite-sample performance similar to that observed in Designs 1--3 in almost all aspects. However, it appears that $n_1=2000$ is not always sufficient for FOCI to attain a desirable CRR in some designs. Taken together with our earlier discussion on computational efficiency, these findings underscore the advantage of our approach in nonparametric estimation involving discrete response variables, such as propensity score estimation. The CV refinement results in Table (ref) reveal similar patterns to those in Table (ref). In particular, while our “Modified-CV” algorithm continues to outperform the “CV” procedure of HallRacineLi, the performance gap narrows in these designs.

Following YangDing2018, we focus on the subpopulation whose propensity scores lie within the interval \(\left[0.1,0.9\right]\). This practice helps avoid issues with extreme estimated propensity scores by trimming the sample. To estimate the ATE \( \psi \), we apply the doubly robust estimator defined in ((ref)) and implement the multi-layer perceptron (MLP) procedure employed in FarrellEtal2021. Specifically, we define \(\hat{\omega} = \mathbf{1}\left(0.1\leq \hat{m}\left(\cdot\right)\leq 0.9\right)\) and consider the following variants:

When implementing MLP to estimate the nuisance functions, the network architecture consists of two hidden layers with 64 and 32 neurons, respectively, each followed by a ReLU activation function. Training is performed using the Adam optimizer with a learning rate of $10^{-3}$, a mini-batch size of 32, and 100 training epochs. We use the mean squared error loss function for outcome regression and the binary cross-entropy loss function for propensity score estimation. All computations are conducted using the PyTorch framework.

For all these estimators, we report their mean bias (Bias), root mean squared error (RMSE), and the length (IL) and coverage probability (CP) of the 95% confidence interval. We summarize these results in Tables (ref)--(ref) in Appendix (ref), corresponding to Designs 4--6, respectively. Estimators $\hat{\psi}_3$ and $\hat{\psi}_4$ exhibit similarly satisfactory performance across all designs. As $n$ increases, both estimators display vanishing biases, RMSEs that decline at approximately $\sqrt{n}$ rate, and shrinking confidence intervals with coverage probabilities approaching the target value of 0.95. Their performance, as expected, slightly deteriorates as $p$ increases. Between the two, $\hat{\psi}_4$ tends to have slightly higher RMSE and wider confidence interval than $\hat{\psi}_3$. In contrast, $\hat{\psi}_1$ performs poorly, yielding larger bias, higher RMSE, and wider confidence interval compared to the other estimators. Estimator $\hat{\psi}_2$ also suffers from relatively large bias and, more critically, produces narrow confidence intervals when $p$ is large, resulting in poor coverage. Overall, we recommend that practitioners employ $\hat{\psi}_3$ or $\hat{\psi}_4$ in empirical applications, as we will do in the next section.

Empirical Illustration

There is considerable debate in empirical studies on whether 401(k) eligibility increases (“crowd-in”) or decreases (“crowd-out”) the accumulation of other assets (savings). A key challenge in identifying this effect is that 401(k) eligibility may be correlated with an individual's unobserved savings preferences, potentially leading to selection bias. For example, individuals with strong savings preferences may be more likely to seek employment at firms that offer 401(k) plans. benjamin2003does investigates this question using regression within propensity score subclasses, with the validity of his analysis relying on the absence of unobserved confounding factors. gelber2011401 aims to address this challenge by employing a difference-in-differences (DID) approach using matched observations based on stratified propensity scores.\footnote{gelber2011401's DID approach estimates the ATT of 401(k) eligibility.} However, both the linear DID regression and the Probit propensity score model employed in his analysis may suffer from model misspecification. In this section, we re-analyze the effect of 401(k) eligibility on savings by applying our proposed dimension-reduced propensity score estimator (Section (ref)) and the doubly robust ATE estimator leveraging Proposition (ref) (Section (ref)) to the rich longitudinal dataset used by gelber2011401, attempting to address the methodological limitations in existing studies.

Broadly, evaluating the effect of tax-advantaged retirement savings plans, such as the 401(k) in the US or pension systems in other countries, on private saving or debt, along with its influence channels, has attracted long-standing interest. Related topics have been extensively studied in recent years using various empirical strategies (e.g., DID, regression discontinuity (kink) designs), as illustrated by works such as chetty2014active, andersen2018tax, messacar2018crowd, goodman2020catching, beshears2022borrowing, and garcia2024public. Identifying and estimating the corresponding causal effects also contributes to a better understanding of individual responses to policy incentives for private retirement savings (chan2022income) and the broader macroeconomic implications of pension policy for the structure of the financial system (scharfstein2018presidential). The related literature is vast. We refer interested readers to the cited studies and the references therein for a more comprehensive review.

Data

We use data from Waves 3, 6, 7, 9, and 12 of the 1996 SIPP.\footnote{The raw data and accompanying do-files used for processing were obtained from gelber2011401's replication package, which is publicly available at \url{https://www.openicpsr.org/openicpsr/project/116534/version/V1/view}.} 401(k) eligibility is observed in Wave 7, and asset (or liability) information is collected in the other waves.\footnote{For simplicity, we treat all liabilities as assets and represent them in positive terms, so that a larger value indicates a larger liability.} Following gelber2011401, we define the period from Wave 3 to Wave 6 as “Year 0”, the period from Wave 6 to Wave 9 as “Year 1”, and the period from Wave 9 to Wave 12 as “Year 2”, which approximately correspond to calendar years 1997, 1998, and 1999, respectively. The data timeline is illustrated in Figure (ref), sourced from gelber2011401. Following prior studies, our sample is restricted to reference persons in each household who are 22–64 years old and work at a for-profit firm in Year 1.

As in benjamin2003does and gelber2011401, our treatment variable $D$ equals 1 if the individual is eligible for a 401(k) plan and 0 otherwise. In the literature, eligible individuals with positive 401(k) balances are defined as 401(k) participants, while those with a zero balance are considered non-participants. It is well documented in the literature that firms decide whether to offer employees 401(k) plans, whereas participation is the employee's choice. Therefore, eligibility, rather than participation, is a more plausible treatment variable that is conditionally independent of individual-level potential outcomes, given observed individual, household, and firm characteristics, satisfying Assumption (ref)(1) (ignorability).

We follow gelber2011401 in examining the following asset categories: 401(k) assets, IRA assets, other financial assets, secured debt, unsecured debt, and car value. Let $A_{t}^{k}$ denote the balance of asset category $k$ in Wave $t$. Our outcome variables are the first differences of asset balances, defined as $\Delta A_{2}^{k} \equiv A_{12}^{k} - A_{9}^{k}$, measuring the change in the balance of asset $k$ in Year 2. This transformation helps control for time-invariant unobserved heterogeneity, such as individual saving preferences. To mitigate the influence of large outliers commonly observed in asset data, we follow gelber2011401 and winsorize outcome variables at the 5th and 95th percentiles. Summary statistics for these outcome variables are presented in Panel A of Table (ref). The reported $p$-values, from regressions of each $\Delta A_{2}^{k}$ on $D$ without other covariates, indicate that eligible and ineligible individuals differ significantly in their accumulation of 401(k) and IRA assets. The latter suggests a “crowd-in” effect, a finding consistent with that of gelber2011401.

To strengthen the plausibility of Assumption (ref)(1), we include a comprehensive set of covariates, encompassing all those used by benjamin2003does and gelber2011401, along with many additional variables.\footnote{These include the number of children under 18, indicators for home ownership, gender, race, metropolitan residence, state of residence, private and public health insurance coverage, financial aid receipt, union or employee association membership, lifetime armed forces service, among others.} Leveraging the longitudinal nature of the data, we also control for changes in 401(k) and other asset balances in Year 0. In total, we include 93 covariates in our analysis.\footnote{These are all stand-alone variables, not including any transformation such as polynomial or interaction terms.} All monetary variables are measured in 1996 US dollars (in units of \$1,000). After excluding observations with missing values, our working data comprises $n = 6741$ observations, of whom 59.98% are 401(k) eligible ($D=1$) and 40.02% are ineligible ($D=0$).

Empirical Model

Given a set of covariates $\boldsymbol{Z}$, we consider the following empirical model for a representative individual's (potential) asset holdings in wave $t \in \{3, 6, 9, 12\}$:

where $A_{t}^{k}(d)$ denotes an individual's holdings of asset $k$ in Wave $t$, given the 401(k) eligibility status $D = d \in \{0, 1\}$. The term $\alpha_{d}^{k}$ captures unobserved time-invariant confounding factors, $u_{t,d}^{k}$ is an idiosyncratic error, and $G_{t,d}^{k}(\cdot)$ is an unknown smooth function. benjamin2003does estimates model ((ref)) using cross-sectional SIPP data, assuming a linear specification for $G_{t,d}^{k}(\boldsymbol{Z})$ and excluding the unobserved confounder $\alpha_{d}^{k}$. Taking first differences between Waves 9 and 12 removes $\alpha_{d}^{k}$ and yields

which defines the potential outcome for our analysis. Then, the observed outcome can be expressed as $\Delta A_{2}^{k} = \Delta A_{2}^{k}(1) D + \Delta A_{2}^{k}(0) (1 - D)$, following the standard potential outcomes framework.

Suppose Assumption (ref) holds in this application. We apply the approaches proposed in Sections (ref) and (ref) to the observed data $(\Delta A_{2}^{k}, D, \boldsymbol{Z})$ to estimate the ATE:

where $m(\boldsymbol{Z}) = \Pr(D | \boldsymbol{Z})$ and $g_{d}^{k}(\boldsymbol{Z}) = \mathbb{E}[\Delta A_{2}^{k}(d) | \boldsymbol{Z}]$ for $d\in\{0,1\}$, consistent with the notation in previous sections.

Variable Selection

We begin by applying the variable selection procedure from Section (ref) to evaluate the dependence between each pre-treatment covariate and $D$. For the procedure's implementation, we set $\tilde{p}=8$ and use the change in an individual's 401(k) balance in Year 0, i.e., $\Delta A_{0}^{\textrm{401(k)}}\equiv A_{6}^{\textrm{401(k)}}-A_{3}^{\textrm{401(k)}}$, as $\boldsymbol{W}$. This variable serves as a proxy for 401(k) participation in Year 0 and is therefore expected to be highly relevant for predicting subsequent eligibility. By capturing prior participation behavior, this variable serves as a sufficient surrogate for many observed and unobserved characteristics that influence 401(k) eligibility, thereby simplifying the dimension reduction task. To minimize subjective judgment, no other covariates are included in $\boldsymbol{W}$.

To ensure a robust selection, we apply our procedure to the full sample as well as to two equal-sized (randomly) split samples ($I_1, I_2$), which yields slightly different sets of selected covariates conditional on $\Delta A_{0}^{\textrm{401(k)}}$. In accordance with the principle of conservatism, our final set of relevant covariates for subsequent ATE estimation, denoted by $\boldsymbol{Z}^{\hat{*}}$, is the union of the variables selected from all three analyses. Definitions and summary statistics for these covariates are provided in Panel B of Table (ref). The reported $p$-values show significant imbalances between the eligible and ineligible groups for most of these characteristics. Notably, our choice of $\tilde{p}=8$ proves to be sufficiently large, as there are always some covariates “smoothed out” in the CV-refining step of our procedure across all three analyses.

Our procedure selects eight covariates as relevant to one’s 401(k) eligibility. As expected, firm characteristics are important. Larger firms and those in certain industries\footnote{According to the 1996 SIPP data documentation, $\texttt{Industry4}$ corresponds to firms engaged in transportation, communication, and utilities.} are more likely to offer 401(k) plans. Among individual-level variables, prior changes in the balances of 401(k) accounts, other financial assets, and unsecured debt may serve as plausible proxies for unobserved saving preferences. Gender and income are also selected, likely because they are important determinants of occupation and employment status. These selected covariates are also used in benjamin2003does and gelber2011401 for propensity score estimation. Notably, our method identifies private health insurance as a strong predictor of 401(k) eligibility, a factor that has been overlooked in both prior studies. This finding highlights the importance of data-driven variable selection in uncovering economic relationships and enhancing model specification.

ATE Estimation and Balance Tests

We first estimate $\psi^{k}$ by applying the estimators $\hat{\psi}_{3}$ and $\hat{\psi}_{4}$, as defined in Section (ref), using the full sample. For all procedures, we standardize all continuous covariates and use the trimming rule $\hat{\omega} = \mathbf{1}(0.05 \leq \hat{m}(\boldsymbol{Z}^{\hat{*}}) \leq 0.95)$ to handle extreme propensity scores. For the MLP implementation, we use a neural network with three hidden layers containing $(64, 32, 16)$ neurons, respectively. The training is conducted with a learning rate of $10^{-5}$ and a batch size of 64. All other settings follow those in Section (ref). These estimation results are reported in Table (ref). In addition, we consider two alternative estimators. First, we implement the H{\'a}jek form of the doubly robust estimator proposed in robins2007comment. Second, we apply the DML procedure of ChernozhukovEtal2018 using a two-fold data partition $(I_{1}, I_{2})$. These results are presented in Tables (ref) and (ref), respectively, for comparison.

A correctly specified propensity score, $m(\boldsymbol{Z}^{*})$, should balance the distribution of covariates between the treatment and control groups, i.e., $D\perp \boldsymbol{Z}^{*}|m(\boldsymbol{Z}^{*})$. Accordingly, a reliable estimator should yield approximately balanced samples within subclasses defined by the estimated propensity scores. Table (ref) reports the $p$-values from two-sample $t$-tests used to assess covariate balance within given subclasses for both the kernel and MLP propensity score estimators. Subclasses corresponding to intervals $[0, 0.1]$ and $[0.9, 1.0]$ are omitted due to trimming. An entry of “NA” indicates that the subclass contains too few treated observations to perform a meaningful test. The results show that the kernel estimator passes the balance test for most covariates and subclasses, and substantially outperforms the MLP estimator in this regard. In view of this, our subsequent interpretation of the empirical results focuses on the ATE estimates obtained using the kernel-based propensity score, as presented in Panel A of Tables (ref), (ref), and (ref).

Empirical Results and Concluding Remarks

Our ATE estimates across different estimators yield largely similar results. As expected, 401(k) eligibility significantly increases contributions to 401(k) accounts. In addition, we find weak evidence suggesting that eligibility leads to increases in IRA savings, other financial assets, and secured debt, which indicates potential “crowd-in” effects. The estimated ATEs for these outcomes are only marginally significant. However, due to large standard errors and consequently wide confidence intervals, we cannot rule out the possibility of economically meaningful effects in these categories. Finally, we find no evidence that 401(k) eligibility significantly affects unsecured debt or car values in either a statistical or an economic sense. A practical limitation of our analysis is the imprecision of the ATE estimates, which may stem from the large noise in asset (liability) data, as reflected by the large standard deviations reported in Table (ref).

The approach we illustrate here allows for unobserved heterogeneity and flexibly accommodates pre-treatment covariates to account for observed heterogeneity. This offers a more robust framework for estimating treatment effects compared to the common practice of using linear or logistic regressions with homogeneous coefficients. The results presented in Panel A of Table (ref) provide a ready example. The reported $p$-values essentially summarize estimation results from two-way fixed effects DID regressions without covariates. These regressions do not provide suggestive evidence for the potential “crowd-in” effects of 401(k) eligibility on other financial assets and secured debt that are detected by our proposed methods.\footnote{A direct comparison should be made with caution, as DID estimates the ATT, while our analysis focuses on the ATE. A more direct comparison would involve integrating our methods into doubly robust ATT estimation (e.g., shinozaki2015brief) or modern doubly robust DID frameworks (e.g., sant2020doubly). This is beyond the scope of this paper, and we leave it for future research.} This highlights the value of adopting a more robust, data-driven approach for uncovering important causal effects.

As a final remark, while our procedure is data-driven, it should be applied under the supervision of researchers with domain knowledge and a degree of conservatism. Once the assumed sparsity structure (like in Proposition (ref)) in the data can be justified, as in our application, nonparametric estimation in fact becomes more straightforward to implement than many commonly used parametric approaches. Prior studies, such as benjamin2003does and gelber2011401, expend significant effort determining and justifying the best specification for their empirical models. This challenge is doubled in causal inference, which often requires specifying both a propensity score and an outcome regression model. A key advantage of our approach is therefore clear: by focusing on identifying the relevant variables, our method enables us to sidestep the need to argue for a specific functional form, thereby simplifying the empirical analysis while maintaining robustness.

Conclusion

This paper introduces a novel method for nonparametric conditional density estimation in high-dimensional settings, addressing the curse of dimensionality through a computationally efficient approach. Rather than imposing strong distributional or functional form assumptions, our dimension reduction approach leverages a sparsity assumption that only a small subset of covariates is truly relevant. The core of our methodology is a screening procedure based on a new measure of conditional dependence, which achieves a favorable balance between small-sample statistical power and computational efficiency. Simulation studies show that it offers a substantial improvement over existing approaches in terms of recovery rates and processing time. Our method is particularly well-suited for applications with moderate sample sizes, where alternative methods may be either computationally costly (e.g., CDSIS by WangEtal2015) or lack desirable small sample power (e.g., FOCI by azadkia2021simple).

We further propose to refine the screening process using a modified CV algorithm that provides a data-driven approach to effectively eliminate irrelevant covariates through bandwidth selection, thereby refining the conditioning set for subsequent density estimation. Simulations demonstrate that our algorithm consistently outperforms the original CV procedure by HallRacineLi in terms of exact recovery probability and numerical stability. Combined with the screening step as a whole, our proposed variable selection procedure provides a practical and theoretically grounded solution to the curse of dimensionality, which often hinders the application of nonparametric methods in complex, high-dimensional settings.

We demonstrate the practical utility of our dimension-reduced density estimator by applying it to the estimation of ATE with doubly robust methods. We show that our dimension reduction procedure offers a dual benefit in this context: it not only refines the set of variables for estimating the propensity score, mitigating the practical challenges posed by the trade-off between ignorability and overlap conditions, but also, through a refined ignorability condition, effectively reduces the dimensionality of the outcome regression models. Simulations confirm our method's potential to enhance the robustness and precision of causal inference in the presence of many covariates.

We further illustrate the real-world applicability of our approach through an empirical analysis of the effect of 401(k) eligibility on savings. Aiming to address unobserved confounders and potential model misspecification issues in prior literature, this application underscores how our proposed dimension-reduced propensity score estimator and doubly robust ATE estimator can lead to more credible empirical findings.

In conclusion, this paper provides a robust and efficient framework for nonparametric conditional density estimation in high-dimensional contexts. Beyond improving the performance and applicability of established causal inference methods, our approach can enhance practitioners’ ability to uncover data structures and address research questions. While we focus on ATE estimation, the proposed procedure offers a general-purpose toolkit for advancing a broad class of econometric and statistical models where understanding conditional distributions is essential.

We are grateful for the valuable feedback and discussions provided by seminar participants at the University of Queensland, as well as conference attendees at the 19th International Symposium on Econometric Theory and Applications (SETA 2025). All errors are our responsibility.