Nonparametric estimation of conditional densities by generalized random forests

Nonparametric Estimation of Conditional Densities by Generalized Random Forests

\def\spacingset#1{ {#1}} \spacingset{1}

\if00 {

I did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors. I report that there are no competing interests to declare. }\\ Deparment of Economics,\\ University of Nebraska--Lincoln} } \fi

\if10 {

} \fi

{\it Keywords:} conditional density, information projection, generalized random forests, high-dimensional infinite-order $U$-statistic, asymptotic properties.

\spacingset{1.8}

Introduction

The conditional density plays a key role in econometrics, as it provides a comprehensive description of how a random variable behaves when a specific value of another variable is given. Counterfactual exercises and policy recommendations are usually based on this description. To have a convincing assessment in this regard, it is crucial to rely on a conditional density estimator with good finite-sample performance and theoretical guarantees that rely on realistic assumptions.

Consider a continuous random vector $(Y , X)$ taking values on $[0,1] \times \mathscr{X}$, where $ \mathscr{X} \subset \mathbb{R}^{d}$ is a compact subset with nonempty interior and $d \in \mathbb{N}$, and let $f(\cdot | x)$ denote the conditional density of $Y$ given $X=x$. In this setting, I propose a nonparametric estimator of $f(\cdot | x)$. Such an estimator is built by combining atw19's forest-based design with the exponential-series approach adopted by bs91 and wu10 to estimate unconditional densities. The combination of these approaches can be briefly described as follows. Let $\boldsymbol{\phi}$ be a $J$-dimensional vector of Legendre basis functions on $[0,1]$. Consider the conditional mean $\boldsymbol{\mu}_x : = \mathbb{E} [ \boldsymbol{\phi} (Y) | X = x ] $ together with the information projection of $f(\cdot |x)$ onto the corresponding exponential family, which is defined by $\tilde{f} (y ; \boldsymbol{\theta}_x) : = {\exp \left[ \boldsymbol{\theta}_x^\tau \boldsymbol{\phi} ( y ) \right]} / { \int_0^1 \exp \left[ {\boldsymbol{\theta}}_x^\tau \boldsymbol{\phi} (t) \right] d t }$ for $y \in [0, 1]$ and with ${\boldsymbol{\theta}}_x \in \mathbb{R}^J$ satisfying $ \int_0^1 [ \boldsymbol{\phi} (y) - \boldsymbol{\mu}_x ] \exp [ {\boldsymbol{\theta}}_x^\tau \boldsymbol{\phi} (y) ] d y = 0 $, i.e., matching the conditional moments. It follows from existing results that $\tilde{f}( \cdot ; \boldsymbol{\theta}_x )$ approximates $f(\cdot | x)$ as $J \rightarrow \infty$; hence, $\boldsymbol{\theta}_x$ can provide a good finite-dimensional summary about the heterogeniety of $f(\cdot | x)$ across $x$.

The proposed estimator $\hat{f}(\cdot | x)$ arises naturally from the precedent discussion and can be computed in two steps as follows. Let $\{ (Y_1, X_1) ,\dots, (Y_n , X_n) \}$ be a random sample from $(Y , X)$. In the first step, we estimate $\boldsymbol{\mu}_x$ by $\hat{ \boldsymbol{\mu}}_x = \sum_{i=1}^n \omega_i (x) \boldsymbol{\phi} ( Y_i ) $, where $\omega_i (x)\in \mathbb{R} $ are weights generated by atw19's random forest algorithm. Specifically, each $\omega_i (x)$ is defined as the fraction of trees in which $X_i$ appears in the same leaf as $x$. So, the weights are adaptive and generated via recursive partitioning on subsamples, where in each split we focus on maximizing the heterogeneity of $\boldsymbol{\theta}_x$ across $x$.\footnote{I refer to has09, as well as ha22, for precise definitions of machine learning concepts such as tree, leaf, and random forest.} In the second step, we set $\hat f(y|x) = {\exp \left[ \hat{\boldsymbol{\theta}}_x^\tau \boldsymbol{\phi} (y) \right]} / { \int_0^1 \exp \left[ \hat{\boldsymbol{\theta}}_x^\tau \boldsymbol{\phi} (t) \right] d t }$, where $ \hat{\boldsymbol{\theta}}_x \in \mathbb{R}^{J}$ is obtained by solving the nonlinear system $ \int_0^1 [ \boldsymbol{\phi} (y) - \hat{\boldsymbol{\mu}}_x ] \exp [ \boldsymbol{\theta}^\tau \boldsymbol{\phi} (y) ] d y = 0 $ with respect to $\boldsymbol{\theta} \in \mathbb{R}^J$.

The distinguishing feature of this approach is the use of adaptive weights based on a random forest design.\footnote{Indeed, with the aim of estimating an auction model with risk-averse bidders, zin18 has suggested using kernel weights in the first step instead of $\omega_i (x)$.} Such a choice is motivated by the increasing success of machine learning techniques and, in particular, random forest algorithms in empirical applications. Among other advantages, these algorithms have succeeded in dealing with the curse of dimensionality and screening out irrelevant covariates wa18. They have also shown good predictive power in supervised learning schemes. In econometrics, machine learning techniques has been applied to a wide range of models that involve, e.g., demand estimation ba15machine, as well as estimating heterogeneous treatment effects in regression discontinuity designs as16class.

Regarding the asymptotic properties, I show that $\hat{f}(\cdot | x ) $ is uniformly consistent on $[ 0 , 1]$ and pointwise asymptotically normal as $J$ increases with $n \rightarrow \infty$. I also provide a ratio-consistent variance estimator to assess the precision of $\hat{f}(y | x ) $ and to build asymptotically valid confidence intervals for ${f}(y | x )$. To derive these results, first, I characterize $\hat{ \boldsymbol{\mu}}_x $ as a high-dimensional infinite-order $U$-statistic with a random kernel. Second, I adapt the arguments in bs91 to conditional densities and extend the results of wa18 and atw19 to allow $J \rightarrow \infty$, relying on inequalities from vit92 for the variance of $U$-statistics.

The proposed estimator $\hat{f}(\cdot|x)$ aims to help the applied researcher in estimating not only $f(\cdot | x)$, but also infinite-dimensional parameters in nonparemetric models such as the ones studied in mat03,mat07,mat08id,mat15. Empirical auction models constitute another example where the proposed estimator $\hat{f}(\cdot|x)$ can be useful, as the parameters of interest can be usually expressed as a functional of the conditional density of bids. Among others, ah07 and pv21 provide comprehensive surveys on this topic. In this paper, I apply the proposed estimator to a real-world dataset of timber auctions from the U.S.\ Forest Service and, in the Supplementary Material, I provide nonparametric estimates of the conditional density of bids given certain characteristics of the auctioned timber track.

There is an extensive body of literature on conditional density estimation. Kernel-based methods have constituted the traditional approach; however, their poor finite-sample performance under many covariates has triggered the developments of bandwidth selection methods and alternative techniques. To name a few, harali04 develop a cross-validation bandwidth selection method that is able to distinguish relevant from irrelevant covariates. efro07 studies adaptive optimal rates using orthogonal polynomials for one-dimensional continuous variables. shen2016 studies the frequentist properties of nonparametric Bayesian models, focusing on adaptive density regression in high-dimensional settings; they develop priors based on orthogonal polynomials and splines, achieving adaptive optimal contraction rates. In a similar setting, nor2017 and nor2022 use priors based on mixtures of regressions, addressing the issue of irrelevant covariates. In addition, kernel-based estimators also suffer from boundary bias problems and, in this regard, series and local polynomial techniques constitute effective solutions: see, e.g., ca22 and the references cited therein.

This paper contributes to this literature but differs from the aforementioned references in that it does not examine the effect on the convergence rate of the smoothness of $f(y |x)$ with respect to $x$. This limitation commonly arises in the study of the asymptotic properties of forest-based estimators, as in wa18 and atw19, although there are recent notable exceptions such as mou20 and cat23inf. Thus, this paper is closely related to recent studies that incorporate machine learning techniques into the development of conditional density estimators such as dal20 and gh22. The distinguishing aspect of my paper, compared to these articles, is the use of atw19's forest design and the establishment of desired asymptotic properties. Additionally, the use of a series method ensures that the uniform consistency result applies to the entire support, preventing the proposed estimator from experiencing boundary bias problems.

The rest of this paper is organized as follows. Section (ref) introduces the proposed estimator together with the algorithm to compute its weights. Section (ref) presents the asymptotic properties and a ratio consistent estimator of the asymptotic variance. Section (ref) reports the results of Monte Carlo experiments. Section (ref) concludes with a discussion of possible extensions. The proofs of the lemmas, theorems, and corollaries are relegated to the Supplementary Material, where I also provide additional results from Monte Carlo experiments and an empirical illustration. The following notation will be employed hereafter.

\paragraph{Notation.} An array with its elements separated by commas --such as $v = (v_1,\dots, v_m)$\sloppy -- is always considered a column vector. The super-script $^\tau$ denotes the transpose, $\mathbf{1}_m$ is an $m$-dimensional column vector of ones, $\mathbf{I}_m$ stands for the $m \times m$ identity matrix, and $\mathbb{N}_0 = \mathbb{N} \cup \{ 0 \}$. The Euclidean and sup- norms of a real vector $v$ are denoted by $\| v \|_2$ and $\| v \|_\infty$, respectively. With a slight abuse of notation, for a matrix $w$, $\| w \|_2$ will denote the matrix norm induced by the Euclidean norm, i.e., $\| w \|_2 = \sup_{\| v \|_2 = 1} \| w v\|_2$. For a real-valued function $\varphi$ defined on a set $T$, the notation is $\| \varphi \|_{T,2} = ( \int_{T} \varphi (t)^2 dt )^{1/2}$ and $\| \varphi \|_{T,\infty} = \sup\{ | \varphi(t) | : t \in T \}$. The interior of $T$ is denoted by $\mathrm{int}(T)$. For an arbitrary countable set $T$, $| T| $ stands for the number of elements. The lexicographical order is our default order relation. So, e.g., the first two elements of $\{ (4,1) , (1,2) , (3,0) , (5,6) \}$ refer to $(1,2)$ and $(3,0)$. Given an ordered set $T = \{ t_1 < t_2 < \dots \}$, we write $( v_t )_{t \in T} = ( v_{t_1}, v_{t_2},\dots )$. Whenever possible, the dependence of certain terms --such as $s$, $J$, $N$, and $\sigma$-- on the sample size $n$ will be omitted from the notation. All asymptotic results are derived as $n \rightarrow \infty$, unless otherwise stated. The symbol $\overset{d}{\rightarrow}$ stands for convergence in distribution and w.p.a.1 abbreviates with probability approaching one. Moreover, $\mathcal{N}(0,1)$ and $z_a$ denote a standard normal distribution and its $a$-quantile, respectively. With a slight abuse of notation, $\mathbb{V}$ denotes either the variance of a random variable or the variance-covariance matrix of a random vector. The convention $ 0 / 0 = 0$ is adopted.

The estimator

Let $Z = (Y,X)$ be a continuous random vector supported on $[0,1] \times\mathscr{X}$, where $\mathscr{X} \subset \mathbb{R}^{d}$ is a compact subset and $d \in \mathbb{N}$. Assume that $\mathrm{int}( \mathscr{X} )$ is nonempty and that the marginal p.d.f.\ of $X$ is continuous and bounded away from zero on $\mathscr{X}$. Throughout this paper, consider a fixed $x \in \mathrm{int}( \mathscr{X} )$ and let $ f(\cdot | x)$ denote the conditional density of $Y$ given $X=x$.

The objective of this section is to build a forest-based estimator of $ f(\cdot | x)$ from a random sample $\{ Z_1 ,\dots, Z_n \}$ of $Z$, being $Z_i = (Y_i , X_i)$. I start with an assumption that imposes smoothness conditions on $ f(\cdot | x)$, similar to those in wu10. Let ${\operatorname{\text{\usefont{U}{BOONDOX-cal}{m}{n}m}}} \geq 2$ be an integer.

Two remarks are noteworthy. First, this assumption introduces different degrees of smoothness for $f(y | x)$ along the $y$- and $x$- directions: for a fixed $x$, we assume that $f(\cdot | x)$ has $m$ continuous derivatives, while $f(y | x) $ only needs to satisfy Lipschitz continuity when it is considered a function of both $( y , x)$. This asymmetric treatment arises because the subsequent analysis will not focus on the effect of the smoothness of $f( \cdot | \cdot)$ on the convergence rate of the proposed estimator, nor aim to achieve the optimal convergence rate. Consequently, I impose the minimal regularity condition necessary to ensure that such an estimator remains asymptotically unbiased.\footnote{ This approach might be particularly convenient, for instance, when estimating a structural model where exogenous covariates are present in the data, but they are excluded from the theoretical analysis and the empirical content does not provide guidance on the smoothness of certain equilibrium outputs with respect to these covariates. }

Second, the support condition $(Y,X) \in [0,1] \times \mathscr{X}$ is standard in the context of series-based estimators and it can be replaced with a more general statement of the form $(Y,X) \in\ \{ (y, x) \in \mathbb{R}^{1 + d} : \ushort{{\operatorname{\text{\usefont{U}{BOONDOX-cal}{m}{n}y}}}} (x) \leq y \leq \bar{{\operatorname{\text{\usefont{U}{BOONDOX-cal}{m}{n}y}}}}(x) , \ x \in \mathscr{X}\}$, where $\ushort{{\operatorname{\text{\usefont{U}{BOONDOX-cal}{m}{n}y}}}} (\cdot) < \bar{{\operatorname{\text{\usefont{U}{BOONDOX-cal}{m}{n}y}}}}(\cdot)$ are continuously differentiable real-valued functions on $\mathrm{int}( \mathscr{X} )$. In this context, if $\ushort{{\operatorname{\text{\usefont{U}{BOONDOX-cal}{m}{n}y}}}} $ and $\bar{{\operatorname{\text{\usefont{U}{BOONDOX-cal}{m}{n}y}}}}$ were known, we can work with the transformation

Thus, throughout this section and the next one, I work under the assumption $(Y,X) \in [0,1] \times \mathscr{X}$ to simplify the exposition. In Section \ref*{sec:ei} of the Supplementary Material, I discuss strategies for handling real-world data in situations where both $\ushort{{\operatorname{\text{\usefont{U}{BOONDOX-cal}{m}{n}y}}}} $ and $\bar{{\operatorname{\text{\usefont{U}{BOONDOX-cal}{m}{n}y}}}}$ are unknown.

Now let $\phi_{\ell}$ be an orthonormal Legendre polynomial on $[0,1]$ defined as follows:

For $J \in \mathbb N$, consider a density $\tilde f$ from the exponential family on $[0,1]$ of the form

where $\boldsymbol{\theta} = ( \theta_1 , \dots , \theta_J) \in \mathbb R^{J}$ and $\boldsymbol{\phi} (y) = \left( \phi_1 (y), \dots, \phi_J (y) \right)$. If we allow $J$ to grow to infinity with $n$ and if $\boldsymbol{\theta}$ is appropriately chosen, then $\tilde f ( \cdot ; \boldsymbol{\theta} ) $ can approximate $f(\cdot | x)$. Specifically, let $\boldsymbol{\theta}_x \in \mathbb R^{J} $ be the vector of the information-projection coefficients that satisfies

The next lemma establishes that $\boldsymbol{\theta}_x$ is indeed well-defined when $J$ is sufficiently large and provides the convergence rate of $\tilde f \left( \cdot ; \boldsymbol{\theta}_x \right)$ towards $f(\cdot | x)$ in terms of ${\operatorname{\text{\usefont{U}{BOONDOX-cal}{m}{n}m}}}$.

The proof of this lemma is provided in the Supplementary Material and it follows essentially by adapting the results of bs91 to conditional densities. I remark that Legendre polynomials are employed as basis functions because of their computational convenience. However, other orthonormal basis can be used instead and, in such a case, the existence and approximation results of Lemma (ref) can be modified accordingly; see, e.g., bs91 and ch07.

From Eqs.\ ((ref)) and ((ref)), I propose estimating $f (\cdot| x)$ by

where $\hat{\boldsymbol{\theta}}_x$ is a forest-based estimator of ${\boldsymbol{\theta}}_x $ that can be computed in two steps as follows. In the first, $\boldsymbol{\mu}_x$ is estimated by generalized random forests. In the second, letting $\hat{\boldsymbol{\mu}}_x$ be the estimator from the previous step, $\hat{\boldsymbol{\theta}}_x$ is obtained by solving the nonlinear system of equations

The rest of this section focuses on the first step as the second is computationally straightforward. In this regard, I remark that existence of $\hat{\boldsymbol{\theta}}_x$ w.p.a.1 is established in Lemma (ref).3 below and that, even though there is no analytical solution for $\hat{\boldsymbol{\theta}}_x$, a numerical solution can be obtained using Newton's method: see wu10 for further discussion.

Consider a positive integer $s <n$, define the set $\mathscr I_{n,s} = \{ ( \iota_1, \dots, \iota_s ) \in \mathbb{N}^s : \iota_1 < \iota_2 < \dots < \iota_s \leq n \}$, and let $\mathscr I_{n,s}^\ast$ be a bootstrap sample from $\mathscr I_{n,s}$ of size $N$: note that $\mathscr I_{n,s}^\ast$ can be obtained by non-replacement subsampling, i.e., by drawing $N$ subsamples of size $s$ from $\{1,\dots,n\}$ without replacement.\footnote{See chk19 for further discussion on this interpretation of the resampling process.} The estimator of ${\boldsymbol{\mu}}_x$ is then defined by

where $\{ \omega_i (x) \in \mathbb{R} : \ i=1,\dots,n \}$ are similarity weights that measure the relevance of the $i$th observation to fitting $\boldsymbol{\theta}_x$, i.e.,

and $\mathcal{L}_x^\ast \left( Z_{\iota_1^\ast} , \dots, Z_{\iota_s^\ast} \right) $ is a subsample of $\{ X_{\iota_1} , \dots, X_{\iota_s} \}$ that contains the observations falling in the same leaf as $x$. Based on bre01's algorithm, a detailed procedure for constructing $\mathcal{L}_x^\ast $ and $\omega_i (x)$ is provided in Subsection (ref) below, together with two suggestions for the splitting schemes.

Choosing the weights to maximize heterogeneity

This subsection provides the algorithm to compute $\mathcal{L}_x^\ast ( z_{\iota_1} , \dots, z_{\iota_s} ) $ via recursive partitioning, for a given training subsample $\{ z_{\iota_1} , \dots, z_{\iota_s} \}$ in which each $z_\iota := (y_\iota , x_{\iota} ) $ can be interpreted as a realization, or data point, of $Z_\iota = (Y_\iota , X_\iota)$. I consider two splitting schemes: one that prioritizes the heterogeneity of $\boldsymbol{\theta}_x$, which is the recommended approach in this paper, and another that targets the heterogeneity of $\boldsymbol{\mu}_x$, offering a low computational cost alternative. I remark that none of these schemes are based on large-sample optimality considerations.

For given data points $\{z_1, \dots, z_n \}$ and $J \in \mathbb N$, every split starts with a rectangular parent node $P \subset \mathrm{int}( \mathscr{X} ) $ containing $x$, and a subsample of indexes $\mathcal{I} = \{ \iota_1 < \iota_2 < \dots < \iota_s \}$. Then we seek to solve an optimization problem of the form

where $ (C_1 , C_2)$ are two axis-aligned children of $P$ and $\mathrm{err}(C_1 , C_2 ; P , \mathcal{I})$ is an error function. In this setting, I consider two splitting schemes whose error functions are as follows. Denote $n ( \mathcal{I} ,C) = | \{ i \in \mathcal{I} : x_i \in C \} |$.

Essentially, each splitting scheme involves a finite-dimensional parameter that summarizes certain aspects of the conditional distribution of $Y$ given $X=x$. In the first scheme, the error function is based on atw19's approach, taking $\boldsymbol{\theta}_x$ as the parameter of interest. The purpose of using the $\arg\min$ set in Eq.\ ((ref)) is to guarantee the existence of a solution to this optimization problem, in a context where $n ( \mathcal{I} ,P)$ does not necessarily grow to infinity as $n \rightarrow \infty$. In the second scheme, we target $\boldsymbol{\mu}_x$ due to its computational simplicity: unlike the first, this approach bypasses the nonlinear optimization problem in Eq.\ ((ref)). It is worth noting that none of these functions guarantee uniqueness of the solution of the optimization problem ((ref)); however, we are just seeking for one solution to this problem, which we know that exists by standard arguments.

Solving optimization problem ((ref)) is often unfeasible in practice. This becomes apparent in the splitting schemes (ref) and (ref) due to the unknown distributions of $X$ and $\boldsymbol{\theta}_X $, along with the extremely high computational cost. Hence, we consider solving an approximated version of problem ((ref)). For that purpose, I introduce an approximating function $\tilde{\Delta}$ that satisfies

where ${\operatorname{\text{\usefont{U}{BOONDOX-cal}{m}{n}t}}} ( P , \mathcal{I} )$ is an estimate of the targeted parameter of the splitting scheme, computed using only the data points $\{ z_i : \ x_i \in P , \ i \in \mathcal{I} \}$. Thus, solving optimization problem ((ref)) becomes approximately equivalent to maximizing $\tilde{\Delta} [ C_1 , C_2 ; \mathcal{I} , {\operatorname{\text{\usefont{U}{BOONDOX-cal}{m}{n}t}}} ( P , \mathcal{I} ) ] $, with respect to $(C_1 , C_2)$.

Before proceeding, I provide the functional forms of $\tilde{\Delta} $ and ${\operatorname{\text{\usefont{U}{BOONDOX-cal}{m}{n}t}}}$ for the Splitting Schemes (ref) and (ref). Since the approximation result in Eq.\ ((ref)) is not needed in the next section for deriving the asymptotic properties of $\hat{f} (\cdot | x)$, I just provide an heuristic discussion in this regard.\footnote{I refer to Sections 2.2 and 2.3 in atw19 for a detailed discussion and formal results related to Eq.\ ((ref)), noting that such results can be applied by considering a fixed $J$: exploring the approximation error as $J \rightarrow \infty$ is beyond the scope of this paper and left for future research.}

To compute $\mathcal{L}_x^\ast ( z_{\iota_1} , \dots, z_{\iota_s} ) $, next I provide a recursive partitioning algorithm that solves $\max_{C_1, C_2} \tilde{\Delta} [ C_1 , C_2 ; \mathcal{I} , {\operatorname{\text{\usefont{U}{BOONDOX-cal}{m}{n}t}}} ( P , \mathcal{I} ) ] $ in each splitting step. Before doing so, I introduce an optimization problem along with input parameters and two randomization devices. Given an integer $\ushort{k} \geq 2$, a constant $\ushort{\alpha} \in (0,1/2)$, a nonempty subset ${\operatorname{\text{\usefont{U}{BOONDOX-cal}{m}{n}d}}} \subseteq \{1,\dots,d\} $, and a parent node $P$ satisfying $| \{ i \in \mathcal{I} : x_i \in P \} | \geq 2 \ushort{k}$, consider the optimization problem

subject to the following constraints:

Note that constraint (c3) implies that a split must put at least a fraction $\ushort{\alpha}$ of the observations of the parent node into the resulting children, $C$ and $P \backslash C$, and that both must have at least $\ushort k$ observations each: see Section 3 and Definition 4 in wa18 for further discussion on this restriction.

Consider also the following randomization devices. Let $\mathcal{W}$ be discrete random vector taking values on $\{ (w_1,\dots, w_s) \in \{ 0, 1\}^s: \ \sum_{m=1}^s w_m = \lfloor s/2 \rfloor \}$ and uniformly distributed over this set. Let $\mathcal{D}$ be a set-valued random element taking values on the power set of $\{1 ,\dots, d \}$ and having a distribution (chosen by the researcher) that satisfies $\mathbb{P} ( \mathcal{D} = \emptyset ) = 0$ and $ \min_{m = 1,\dots, d} \mathbb{P} ( \mathcal{D} = \{ m \} ) > \pi / d$ for some $\pi \in ( 0, 1)$.

The steps for computing $\mathcal{L}_x^\ast ( z_{\iota_1} , \dots, z_{\iota_s} ) $ are provided in Algorithm 1.

For given data points $\{ z_1 , \dots, z_n \}$, now the weights $\omega_i (x)$ can be computed in three steps as follows. First, draw $N$ subsamples of size $s$ from $\{1,\dots, n \}$ without replacement. Second, apply Algorithm 1 to each subsample $\{ \iota_1^\ast < \iota_2^\ast < \dots < \iota_s^\ast \} $ and obtain $\mathcal{L}_x^\ast ( z_{\iota_1^\ast} , \dots, z_{\iota_s^\ast} ) $.\footnote{If the same subsample is drawn twice, or multiple times, the second step must be implemented only one time per subsample. In other words, the described procedure must be applied to each subsample, so we may not need to repeat the second step $N$ times.} Third, compute $\omega_i (x)$ by applying Eq.\ ((ref)) to each $i = 1,\dots,n$.\footnote{If $ x_k \notin \mathcal{L}_x^\ast ( z_{\iota_1^\ast }, \dots, z_{\iota_s^\ast }^\ast ) \ \forall k = 1,\dots, n$ for some subsample $\{ \iota_1^\ast < \iota_2^\ast < \dots < \iota_s^\ast \} $, then follow the convention 0/0=0.}

I highlight that Algorithm 1 does not grow the entire tree, rather it only grows the branch where $x$ resides because we focus on a fixed $x$. Such a branch is grown with a double-sample procedure as the training data is divided into two parts, $\mathcal{I}_0 $ and $\mathcal{I}_1$. The resulting branch is honest because it does not use observations in $\mathcal{I}_0 $ to make the splits in Lines 5-10 of the algorithm. In other words, the branch is grown using one sub-subsample, $\mathcal{I}_1$, while the predictions at the leaf are estimated in another. In addition, the resulting leaf comes from a random-split tree in the sense that the probability that the next split occurs along the $m$th feature of $X$ is bounded from below by $\pi /d$, regardless of the relevance of such a feature and the choice of the splitting scheme. All these properties are crucial to make Assumption (ref) in the next section feasible.

To compute $\mathcal{L}_x^\ast \left( z_{\iota_1} , \dots, z_{\iota_s} \right) $ over a grid of $x$ points, Line 8 of Algorithm 1 should be modified to store each resulting child that contains at least one grid element. Since this would increase the computational cost, we may consider using Splitting Scheme (ref) or a thicker grid in Line 8. However, it is worth noting that using Splitting Scheme (ref) still imposes the same computational burden as the one in Algorithms 1 and 2 of atw19.

To conclude this section, I emphasize that other splitting schemes can be used in Line 8 rather than the ones presented above. In other words, we can chose different approximating function $\tilde\Delta$ for the optimization problem ((ref)). So, in practice and keeping into considerations the computational resources, one can adapt the splitting schemes to target certain characteristics of the parameter of interest, which often can be written as a functional of $f(\cdot|x)$.

Asymptotics

In this section, I establish uniform consistency and (pointwise) asymptotic normality of $\hat f (\cdot | x)$. I also provide a computationally tractable standard error formula to build asymptotically valid confidence intervals.

As a starting point, I introduce a high-level assumption on the tuning parameters, specifically focusing on the diameter of $\mathcal{L}_x^\ast ( Z_{1} , \dots, Z_{s} ) \cup \{ x\}$. This diameter, denoted as $\bar{\mathrm{d}}( \mathcal{L}_x^\ast)$, is defined as the length of the longest segment parallel to one of the axes; namely,

with $\mathrm{d}_m( \mathcal{L}_x^\ast) = \sup \left\{ | x_m^{\prime\prime} - x_m^{\prime} | : \ \{ x^{\prime\prime} , x^{\prime} \} \subset \mathcal{L}_x^\ast ( Z_{1} , \dots, Z_{s} ) \cup \{ x\} \right\}$.

In words, this assumption requires that the terminal leaf has vanishing diameter. This will be used later to derive the asymptotic properties $\hat{f}(\cdot | x)$ since it produces a bound in the bias of $\hat{\boldsymbol{\mu}}_x$, as in Lemma 1 and Theorem 3.2 in wa18. From these results, combined with Algorithm 1, it also follows that Assumption (ref) is automatically satisfied with

To have an idea about the magnitude of ${\operatorname{\text{\usefont{U}{BOONDOX-cal}{m}{n}c}}}_b > 0$, e.g., setting $\ushort{\alpha} = 1/5$ yields ${\operatorname{\text{\usefont{U}{BOONDOX-cal}{m}{n}c}}}_b \approx 0.07 \pi /d$. Concrete suggestions for choosing these tuning parameters will be provided below in Section (ref) (Monte Carlo experiments).

Eq.\ ((ref)) holds essentially for any error function $\mathrm{err}(\cdot , \cdot)$ and any approximating function $\tilde\Delta$. This is made possible by the randomization device introduced in Line 7 of Algorithm 1 that ensures that there will be a split in any direction with positive probability, regardless of the criterion function used in Line 8 and irrespective of the relevance of the covariate. So, under proper choices of $( \ushort{\alpha} , \ushort{k} , \pi )$, Assumption (ref) remains valid even with irrelevant covariates and under practically any choice of $\tilde\Delta$. In comparison to other conditional density estimators, this flexibility comes at the price of having a bias that vanishes at a slow rate.

I also remark that the tuning parameters and the initial parent node can be chosen in any way as long as the condition $\mathbb{P} [ \bar{\mathrm{d}}( \mathcal{L}_x^\ast) \geq s^{-{\operatorname{\text{\usefont{U}{BOONDOX-cal}{m}{n}c}}}_b} ] = \mathcal{O} (s^{-{\operatorname{\text{\usefont{U}{BOONDOX-cal}{m}{n}c}}}_b} ) $ is satisfied. This vanishing diameter condition has been introduced as a high-level assumption for two reasons. First, to maintain generality, acknowledging that different choices of these tuning parameters might yield larger (better) values of $ {\operatorname{\text{\usefont{U}{BOONDOX-cal}{m}{n}c}}}_b $ than the one derived in Eq.\ ((ref)). Second, to simplify the presentation of the asymptotic properties later, as the obtained rates of convergence and the additional required assumptions will be directly expressed in terms of $ {\operatorname{\text{\usefont{U}{BOONDOX-cal}{m}{n}c}}}_b $ without referring to the tuning parameters.

Following mh16's approach, the next step consists in characterizing $\hat{\boldsymbol \mu}_x$ as an infinite-order $U$-statistic with a random kernel. With this aim, I introduce two random vectors indexed by the set $\mathscr{I}_{n,s}$ that will play the role of randomization devices: $(V_{\boldsymbol{\iota}})_{\boldsymbol{\iota} \in \mathscr{I}_{n,s}}$ and $(W_{\boldsymbol{\iota}})_{\boldsymbol{\iota} \in \mathscr{I}_{n,s}}$. Both are assumed to independent between each other and from $(Z_1,\dots,Z_n)$. The role of $(V_{\boldsymbol{\iota}})_{\boldsymbol{\iota} \in \mathscr{I}_{n,s}}$ is to capture the randomness in generating the bootstrap sample $\mathscr{I}_{n,s}^\ast$ or, equivalently, in drawing $N$ subsamples of size $s$ from $\{1,\dots, n\}$ without replacement. Thus, it is assumed that $(V_{\boldsymbol{\iota}})_{\boldsymbol{\iota} \in \mathscr{I}_{n,s}}$ has a multinomial distribution with $N$ trials and all cell probabilities equal to $1/\binom{n}{s}$. The second random vector, $( W_{\boldsymbol{\iota}} )_{\boldsymbol{\iota} \in \mathscr{I}_{n,s}}$, captures the randomness of $( \mathcal{W} , \mathcal{D} )$ in Algorithm 1 for generating the sub-subsample and producing the splits.

Given these randomization devices, I can write

for $\boldsymbol{\iota} = ({\iota}_1, \dots, {\iota}_s) \in \mathscr I_{n,s} $ and

where

I further consider the $U$-statistic

noting that $\mathbf{h} ( z_1,\dots, z_s )$ is a symmetric nonrandom kernel, and I denote by $\tilde{\boldsymbol{\theta}}_{x}$ the solution of the nonlinear system of equations $ \int \boldsymbol{\phi} (y) \tilde f ( y ; \boldsymbol{\theta} ) dy = \mathbb{E} ( \tilde{\boldsymbol \mu}_x )$ with respect to $\boldsymbol{\theta} \in \mathbb R^{J}$, whenever it exists. Next, I can determine the convergence rates of $ \hat{\boldsymbol \mu}_x $, $\tilde{\boldsymbol{\theta}}_x$, and $\hat{\boldsymbol{\theta}}_x$ by introducing the next assumption that determines the rate at which the tuning parameters $(s,J,N)$ should grow towards infinity.

Now I can state the first lemma.

The first part of this lemma establishes that the difference $\hat{\boldsymbol \mu}_x - \tilde{\boldsymbol \mu}_x $ is asymptotically negligible and converges to zero at a very fast rate, which is a consequence of Assumption (ref).(iii). This implies that $\hat{\boldsymbol{\theta}}_x$ can be treated asymptotically as the solution of the system $\int \boldsymbol{\phi} (y) \tilde f \left( y ; \boldsymbol{\theta} \right) dy = \tilde{\boldsymbol \mu}_x $. Then, the results in Eq.\ ((ref)) follows by extending the arguments in wa18 to allow $J \rightarrow \infty$. In particular, the vanishing bias result is a direct consequence of the Lipchitz continuity condition of Assumption (ref) combined with the vanishing diameter condition of Assumption (ref). Existence of both $\tilde{\boldsymbol{\theta}}_x$ and $\hat{\boldsymbol{\theta}}_x $ is obtained by adapting the arguments of Lemma 5 in bs91 to conditional densities, while their approximation rates can be obtained by standard arguments.

As a corollary of Lemmas (ref) and (ref), it follows that $\hat{f} ( \cdot | x ) $ is a uniformly consistent estimator of $f ( \cdot | x)$ on $[0,1]$.

From this corollary, $ \hat{f} ( \cdot | x ) - f ( \cdot | x)$ attains the fastest uniform rate of convergence by setting $\beta = ( 1 + 2 {\operatorname{\text{\usefont{U}{BOONDOX-cal}{m}{n}c}}}_b )^{-1}$ and $\gamma = {\operatorname{\text{\usefont{U}{BOONDOX-cal}{m}{n}c}}}_b[ m ( 1 + 2 {\operatorname{\text{\usefont{U}{BOONDOX-cal}{m}{n}c}}}_b ) ]^{-1}$. This leads a uniform rate of $n^{-( {\operatorname{\text{\usefont{U}{BOONDOX-cal}{m}{n}m}}} - 1) {\operatorname{\text{\usefont{U}{BOONDOX-cal}{m}{n}c}}}_b / [ m ( 1 + 2 {\operatorname{\text{\usefont{U}{BOONDOX-cal}{m}{n}c}}}_b ) ] }$. To have an idea about its magnitude, e.g., if we set $\ushort \alpha = 1/5$ in Eq.\ ((ref)) so that ${\operatorname{\text{\usefont{U}{BOONDOX-cal}{m}{n}c}}}_b \approx 0.07 \pi /d$, the resulting rate becomes $n^{-0.07( {\operatorname{\text{\usefont{U}{BOONDOX-cal}{m}{n}m}}}-1)/[{\operatorname{\text{\usefont{U}{BOONDOX-cal}{m}{n}m}}}(d/\pi +0.14)]}$, which indicates the optimal convergence is not achieved.

With the aim of deriving the asymptotic distribution of $\hat{f} ( \cdot | x )$, hereafter, I consider a fixed $ y \in [0,1]$ such that $\phi_ {{\ell}} ( y) \neq \int \phi_ {{\ell}} (t) f (t | x) dt$ for some ${\ell} \in \mathbb{N} $.\footnote{The purpose of introducing this condition is to avoid super-consistency issues; see, e.g., Theorem 8 in wu10.} Then, the next lemma establishes that the difference $ \hat f ( y | x) - \tilde f (y ; \tilde{\boldsymbol{\theta}}_x ) $ can be approximated by a zero-mean infinite-order $U$-statistic.

Now consider the H\'ajek projection of $\mathcal{U}$, which is given by

Denote its variance by $\sigma^2 = s^2 \mathbb{E} [ g_y (Z_1)^2 ] /n$, noting that its dependence on $y$ is omitted from the notation to simplify the exposition.

The next lemma is the building block for obtaining the desired asymptotic distribution.

The first part of this lemma implies that

while the second establishes essentially a Lyapunov condition, from which the asymptotic distribution of $\mathcal{U}^\circ / \sigma $ can be automatically derived. The proof of these parts are mainly based on the proofs of Theorems 3.3 and 3.4 in wa18. The key technical challenges here consist in adapting wa18's arguments to allow $J \rightarrow \infty$ and to deal with the fact that the Lipschitz constant of the mapping $x^\prime \mapsto \mathbb{E} [ \mathcal{T}_x (y) \boldsymbol{\phi} (Y) | X = x^\prime ]$ can increase with $n$. Part 3.\ of Lemma (ref) establishes that the distributions of $\mathcal{U}$ and $\mathcal{U}^\circ$ must be very similar. The proof of this result relies on $U$-statistic theory vit92 and projections vdVaart98.

Now we can derive the desired asymptotic distribution.

I remark that the reason for introducing extra conditions, $\beta > (1 + 2 {\operatorname{\text{\usefont{U}{BOONDOX-cal}{m}{n}c}}}_b)^{-1}$ and $\gamma < \beta(1 + 2 {\operatorname{\text{\usefont{U}{BOONDOX-cal}{m}{n}c}}}_b) - 1$, is to eliminate the remainder $\mathcal{O}_p $-term from Lemma (ref). The first condition is standard as it is also required for estimating heterogeneous treatment effects wa18 and parameters of interest identified via a local moment equations atw19. Combined with Assumption (ref).(ii), it automatically implies that $(1-\beta)/3 < 2 {\operatorname{\text{\usefont{U}{BOONDOX-cal}{m}{n}c}}}_b /3 $. The second condition, $\gamma < \beta(1 + 2 {\operatorname{\text{\usefont{U}{BOONDOX-cal}{m}{n}c}}}_b) - 1$, just prevents $J$ from increasing too fast. Despite the fact that $J \rightarrow \infty$, the rate at which $\hat f ( y | x)$ approximates $\tilde f (y ; \tilde{\boldsymbol{\theta}}_x )$ can be compared with the one obtained in Theorem 5 of atw19; essentially, they only differ by the term $\| \mathcal{T}_x (y) \|_2 = \mathcal{O}( \sqrt{J})$.

The discussed extra conditions also guarantee that $[ \tilde{f} ( y ; \tilde{\boldsymbol{\theta}}_x ) - \tilde f ( y ; \boldsymbol{\theta}_x ) ] / \sigma = {\operatorname{\text{\usefont{U}{BOONDOX-cal}{m}{n}o}}}(1)$. However, further requirements, namely ${\operatorname{\text{\usefont{U}{BOONDOX-cal}{m}{n}m}}} \geq 3$ and $\gamma > (1- \beta )/[ 2 ( {\operatorname{\text{\usefont{U}{BOONDOX-cal}{m}{n}m}}} - 1 )]$, are still needed to center the asymptotic distribution of $[ \hat f ( y | x) - f ( y | x) ]/\sigma $ at zero by making the approximation error $[ \tilde{f} ( y ; \boldsymbol{\theta}_x ) - f ( y | x ) ]/ \sigma$ negligible.

The applicability of Theorem (ref) stems from the ability to construct confidence intervals using the quantiles of the standard normal as critical values, provided that a valid standard error for $\hat{f}(y|x)$ is available. Specifically, for given a $\alpha \in (0,1)$ and a valid standard error $\mathrm{se} [ \hat{f}(y|x) ] $, the confidence interval

ensures proper asymptotic coverage. The next corollary formalizes this statement.

The next subsection provides an example of a ratio consistent estimator of $\sigma$.

Ratio consistent standard error formula

This subsection provides a ratio-consistent estimator $\hat\sigma^2$ of $\sigma^2$, allowing $\hat\sigma := \sqrt{\hat\sigma^2}$ to serve as a valid standard error for $\hat{f} (y | x)$. With this estimator, we can assess the precision of $\hat{f} (y | x)$ and build asymptotically valid confidence intervals for $f (y |x)$.

Following closely atw19 and combining their approach with the delete-$D$ Jackknife method sw89, I begin by introducing an ideal estimator of $\sigma^2$:

where $D_\sigma \in \mathbb{N}$ is a sequence satisfying $ D_\sigma < n - s$, $\hat{\mathcal{T}}_x (y)$ is the plug-in estimator of ${\mathcal{T}}_x (y)$, i.e.,

and $\tilde{\boldsymbol \mu}_{x,\boldsymbol{\iota}} $ denotes the estimator $\tilde{\boldsymbol \mu}_x$ computed using the subsample $\{ Z_{\iota_1},\dots, Z_{\iota_{n-D_{\sigma}}}\}$ whose indexes are taken from $\boldsymbol{\iota} = (\iota_1 , \dots, \iota_{n-D_{\sigma}}) $, after deleting $D_\sigma$ observations from the sample. To be specific,

with $\check{\mathscr{I}}_{\boldsymbol{\iota},s} = \{ ( \tilde\iota_1, \dots, \tilde\iota_s ) \in \{ \iota_1, \dots, \iota_{n-D_\sigma} \}^s : \tilde\iota_1 < \tilde\iota_2 < \dots < \tilde\iota_s \}$.

The idea behind $\hat\sigma^2$ is to directly approximate the variance of $\mathcal{U}$ via Jackknife methods for $U$-statistics, for which I recall that $E ( \mathcal{U}^2) \approx \sigma^2$ by Lemma (ref).3, and the next theorem formalizes this intuition.

As a consequence of Corollary (ref) and Theorem (ref), for a given $\alpha \in (0,1)$, we have that the confidence interval defined by

provides proper asymptotic coverage, i.e., $\mathbb{P} [ f(y|x) \in \mathrm{CI}_{1-\alpha} ] \rightarrow 1-\alpha$.

I emphasize that allowing the researcher to choose $D_\sigma$ such that $n/D_\sigma \rightarrow c_\sigma > 1$ offers more flexibility in comparison with the half-sampling method adopted in the previous literature, which essentially consists in selecting $D_\sigma = \lfloor n/2 \rfloor$ and $c_\sigma = 2$. Having an additional input parameter to choose can be helpful to improve the finite-sample performance.

Since $\hat\sigma$ is computationally unfeasible, I conclude this section by proposing a modification of the subsample scheme to make it feasible. Specifically, I suggest modifying the subsample scheme as follows:\footnote{See sl09se for further discussion about this computational procedure.}

Then, the estimators $\hat{\boldsymbol\mu}_{x}$ and $\hat{f} ( \cdot | x)$ can be computed by applying the procedure described in Section (ref) to the subsamples $\mathscr{S}_{1} , \dots , \mathscr{S}_{N}$, while a feasible version of $\hat\sigma$ can be computed as follows:

where $\hat{\boldsymbol \mu}_{x, -l} $ denotes the estimator of $\boldsymbol\mu_x$, calculated solely from the subsamples $\mathscr{S}_{1} , \dots , \mathscr{S}_{N}$ that do not contain $\mathscr{S}_{\sigma ,l}$, and ${\boldsymbol \mu}_{x,av} = (1/N_\sigma) \sum_{l} \hat{\boldsymbol \mu}_{x, -l} $. Note that Lines 3-4 of Algorithm 2 ensure a minimum of two such subsamples.

Monte Carlo experiments

This section presents the results of Monte Carlo experiments to evaluate the finite-sample performance of the proposed estimator $\hat{f}( y | x)$ when $d = 4$. I consider the the bias, standard deviation, and the Mean Integrated Squared Error (MISE) as evaluation criteria. In addition, I analyze the performance of confidence intervals that use $\hat\sigma_{\mathrm{fe}} $ as standard error.

The design of the experiments is as follows. The vector of covariates $X \in \mathbb{R}^4$ has a multivariate normal distribution with mean $(1/2) \mathbf{1}_4$ and variance-covariance matrix $(1/8)\mathbf{I}_4$, and truncated on $[0,1]^4$. The following designs are considered for the conditional distribution of $Y$ given $X$, noting that $X_4$ is an irrelevant covariate.

I set $x = (1/2) \mathbf{1}_4$ and the conditional density $f(y| x )$ corresponding to each design can be visualized in Figure (ref) below. The design points for the values of $y$ are provided in the second column of Table (ref) below. Two sample sizes are considered: $n =500 , 1000 $.

I have performed 500 replications for each design and sample size.\footnote{These replications were conducted using the Holland Computing Center’s Swan system at the University of Nebraska, a high-performance computing resource with 56 cores and 256GB of memory.} In each replication, first, I generated a random sample from the corresponding scenario. Then, I computed the estimator $\hat{f} ( \cdot | x)$ and the standard error $\hat\sigma_{\mathrm{fe}}$ using the subsample scheme suggested in Algorithm 2. The tuning parameters of Algorithm 1 were chosen as follows: $J=8$, $N= 2240$, $P = [1/4, 3/4]^4$, $\ushort{k}= 10$, and $\ushort{\alpha} = 0.05$. In addition, I used $s \in \{ 125 , 250 \}$ and $s \in \{ 200, 400 \}$ as subsample sizes for $n=500$ and $n = 1000$, respectively. For choosing the covariates in each splitting step (Line 7 of Algorithm 1), I randomly chose $\min\{ \max\{ \mathrm{Poisson}(5) , 1 \} , 4\}$ integers from $\{ 1, 2, 3, 4\}$ without replacement, so $\min_{m = 1,\dots, d} \mathbb{P} ( \mathcal{D} = \{ m \} ) \approx 0.0101$ and therefore $\pi \approx 0.04$. Line 8 of Algorithm 1 was implemented using a simple grid search. The tuning parameters of Algorithm 2 were $N_\sigma = 560$ and $D_\sigma = n /20$.

I also computed two alternative estimators for comparison purposes. One of them is an oracle kernel estimator of the form $\hat{f}^{ker}( y | x) = \hat{f}_{yx}^{ker}( y , x) / \hat{f}_x^{ker}( x) $ that correctly ignores $X_4$. As tuning parameters for this estimation, I have used a tri-weight kernel and bandwidths of the form $h_y = 1.06\hat{\sigma}_y n^{-1/8}$ and $h_{x_m} = 1.06\hat{\sigma}_{x_m} n^{-1/8} $ in the numerator and $h_{x_m} = 1.06\hat{\sigma}_{x_m} n^{-1/7}$ in the denominator for $m=1,2,3$. The other estimator is a conditional Maximum Likelihood Estimator (MLE) based on a over-fitted model that has six parameters and includes $X_4$. To estimate D1, I used the specifications $b_1 + b_2 X_1 + b_2 X_2 + b_4 X_4$ and $b_5 + b_6 X_3$ for the parameters $\alpha$ and $\beta$, respectively. To estimate D2 and D3, I employed $b_1 + b_2 X_1 + b_3 X_2 + b_4 X_4 $ and $b_5 + b_6(X_3 - 1/2)^2$ for $\mu$ and $\sigma^2$, respectively. The maximization problem was then carried out with respect to $(b_1,\dots,b_6)$ over a compact set, which included the true values of the parameters.

The results of the Monte Carlo simulations are presented in Tables (ref)-(ref) and Figure (ref) at the end of this section. Table (ref) reports the bias and the standard deviation (in parentheses) for the design points of $y$, as well as for the MISE over $[0.15 , 0.85]$. The column `GRF' displays the results linked to the proposed estimator $\hat{f}( y| x)$, while the columns `Kernel' and `MLE' provide the results from the alternative estimators. As can be noted, the bias of $\hat{f}( y| x)$ is relatively small except at certain values of $y$ in Design D3; namely, $y = 0.250, 0.500, 0.750$. Increasing $s$ seems not to have a significant impact on the bias. The standard deviation of $\hat{f}( y| x)$ is extremely small when compared to that of the oracle kernel estimator and, as expected, it is larger than the one obtained from the MLE procedure. In terms of the MISE, the proposed estimator exhibits a very good performance and clearly dominates the oracle kernel estimator. Both the standard deviation and the MISE of $\hat{f}( y | x)$ decreases when $n$ increases.

The performance of $\hat{f}( y | x)$ when $n = 1,000$ can be visualized in Figure (ref). In this figure, the solid line represents the true conditional density for each scenario. The circles and squares correspond to the mean and median of the estimator, respectively, obtained in the simulations for each value of $y$ provided in the second column of Table (ref). As suggested by the symmetric limiting distribution, there is an overlap between the mean and median, and both are close to the solid line: this effect is particularly noticeable in Design D1. The black up- and down-pointing triangles are the $10^\text{th}$ and $90^\text{th}$ percentiles of $\hat{f}( \cdot| x)$, respectively. It is noteworthy that the difference between these two percentiles tends to be relatively small, as one may expect from the values of the standard deviations in Table (ref).

Additional experimental results are provided in the Supplementary Material. Specifically, Table (ref) and Figure (ref) were replicated using subsample sizes of $s \in \{ 75, 200 \}$ and $s \in \{ 125, 350 \}$ for $n=500$ and $n = 1000$, respectively, to conduct robustness checks. The results remain largely consistent with the main findings, indicating that the conclusions are not sensitive to these variations in subsample sizes.

Finally, Table (ref) presents the results to evaluate the performance of the 95% confidence intervals $[ \hat{f}(y | x) \pm \hat\sigma_{\mathrm{fe}} \times z_{0.975}] $ across the design points of $y$. Given the nominal confidence level of 0.95, an ideal confidence interval should achieve coverage close to this value. Among the three designs, D1 and D2 show relatively stable coverage, with some deviations, particularly for smaller and intermediate design points. On the other hand, Design D3 exhibits more variability, with noticeably lower coverage at 0.500 that can be attributed to the bias. Table (ref) also provides the average standard errors (shown in parentheses), which can be compared to the values in parentheses under the ‘GRF’ columns in Table (ref). Overall, $\hat\sigma_{\mathrm{fe}}$ effectively approximates the standard deviation of $\hat{f}(y | x)$, though it tends to overestimate it at higher values of $y$.

Conclusion

I have proposed a nonparametric estimator for the conditional density of $Y$ given $X=x$. The proposed estimator $\hat{f}(\cdot |x)$ is based on atw19's generalized random forest design, targeting the heterogeneity of $\boldsymbol{\theta}_x$ in each splitting step. I have shown that $\hat{f}(\cdot |x)$ is uniformly consistent and asymptotically normal. I also have provided a standard error formula that allows the researcher to build asymptotically valid confidence intervals, using the percentiles of a standard normal as critical values.

I conclude with suggestions on potential avenues for future research. First, the obtained results can be generalized to scenarios where $Y$ is a random vector, i.e., $Y \in [0,1]^{d^\prime}$ with $d^\prime \in \mathbb{N}$. To do so, one can consider the multivariate orthonormal Legendre polynomials on $[0,1]^{d^\prime}$, defined by

together the following multivariate conditional density from the exponential family:

where here $\boldsymbol{\theta} = ( \theta_1 , \dots , \theta_J) \in \mathbb R^{J}$, $\boldsymbol{\phi} (y) = \left( \phi_{\boldsymbol{\ell}_1} (y), \dots, \phi_{\boldsymbol{\ell}_J} (y) \right)$, and $\{\boldsymbol{\ell}_1 < \boldsymbol{\ell}_2 < \dots < \boldsymbol{\ell}_J\}$ are the first $J$ elements of the set $ \mathbb{N}_0^{d^\prime} \backslash \{ (0,\dots,0) \}$. Then, the multivariate conditional density of $Y$ given $X=x$ can be estimated by $\hat f ( y | x) : = \tilde f ( y ; \hat{\boldsymbol{\theta}}_x ) $, where $\hat{\boldsymbol{\theta}}_x$ satisfies

and $\{ \omega_i (x) \in \mathbb{R} : \ i=1,\dots,n \}$ can be obtained by adapting the arguments of Section (ref). The asymptotic findings of Section (ref) can be extended to encompass this multivariate scenario by combining together the results of this paper with those from wu10, which considers a multivariate unconditional density.

Second, in Section (ref), the heterogeneity of $\boldsymbol{\theta}_x$ is targeted by maximizing Eq.\ ((ref)) and using the pseudo-outcomes $\boldsymbol{\rho}_{i} (\boldsymbol{\theta} )$, $i=1,\dots,n$. Since $\boldsymbol{\rho}_{i} (\boldsymbol{\theta} )$ can be treated as a multivariate output, an alternative approach here might be using the methods described in sch23tree that are based on modified splitting or stopping rules for multi-output regressions. Third, an alternative variance estimator can be constructed by combining the approximation of Lemma (ref) with the recent developments on variance estimation for random forests with infinite-order $U$-statistics xu23var. Fourth, it may be helpful for practical purposes to develop a selection rule for $J$, taking into account the computational challenges that this may involve.