Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
84,610 characters · 16 sections · 0 citation commands
Multivariate kernel regression in vector and product metric spaces
This paper extends nonparametric kernel regression to more general regressor settings than those considered in the literature. The general regression model
is free from the difficulty of choosing a parametric specification. We focus on the Nadaraya-Watson (NW) estimator, introduced by Nadaraya (1965) and Watson (1964), which recognizes that a continuous regression function can be estimated pointwise by a weighted average that attaches higher weights to close-by observations.
Here we emphasize the fact that the regressor $X$ can be a vector in $ \mathbb{R}^{q}$, or alternatively $X$ may belong to a function space, a more general metric space, or comprise of components from several such metric spaces. In fact, the components of $X$ do not necessarily have to belong to spaces of vectors or functions but could be intervals, graphs, or networks, as long as a metric (or even a semi-metric) can be defined for each space. $Y$ represents a scalar dependent variable, $u$ denotes an unobserved error, and the conditional mean function $m$ satisfies some smoothness assumptions.
The NW estimator has been used extensively with $X\in \mathbb{R}^{q}$ (see e.g. the textbook Li and Racine, 2007, for discussion and examples) and has recently been introduced to functional regression by Ferraty and Vieu (2004). Well known limit distributional results for the NW estimator were derived in $\mathbb{R}^{q}$ under restrictions requiring existence and smoothness of the density. For functional regression (where there is no density) the limit distributional results were derived in a univariate context only.
In this paper we establish asymptotic normality of the NW estimator for regression in the presence of a general regressor $X$ that could have a multivariate {singular distribution} in $\mathbb{R}^{q}$ or is comprised of any number of functional and vector regressors. A singular distribution does not admit a density function that integrates to it.
In settings where data has both discrete and continuous components Racine and Li (2007) obtained asymptotic normality of the NW estimator without having to deal with the singularity by treating the discrete and (absolutely) continuous components separately. However, sometimes the distinction between discrete and continuous variables is not straightforward; continuous variables could be discretized with different levels of discretization. When data with both discrete and continuous components is viewed as a vector in a Euclidean space, $X\in \mathbb{R}^{q},$ the distribution of $X$ is singular. In our simulations we demonstrate that there may be no gain from avoiding the singularity by considering the discrete regressors separately.
The presence of latent factors in the continuous regressors, common in macroeconomic and finance models (e.g., Bai and Ng, 2006, for portfolio, stock returns and macroeconomic data) could also imply a singular distribution. For example, if the regressor is a \thinspace $q\times 1$ vector $X\sim N\left( 0,\Sigma \right) $ with $\Sigma $ a singular matrix of rank $r<q$, the distribution is singular. Similarly, non-linear common factors, such as in Hotelling's (1929) spatial model of horizontal differentiation which assumes that each consumer has an `ideal' variety identified by his location on the unit circle (see also Desmet et al., 2010), imply a smaller effective dimension for the regressor space, resulting in a singular distribution over $\mathbb{R}^{q}.$ This also is true when there exists a functional relation between the regressors (e.g. with exact collinearity that can arise in production functions, Ackerberg et al., 2015).
Singularity also originates from a fractal structure in the data; examples in economics include the daily prices in the cotton market (Mandelbrot, 1997), financial markets, and networks (see Takayasu et al., 2009). Fractals are common to many geographic features, including coastlines, river networks and landforms, and have been used in urban growth studies (e.g., Shen, 2002) and spatial econometrics in general. Furthermore, singularities also result when continuously distributed variables exhibit mass points (e.g. Arulampalam et al., 2017, for neonatal mortality and Olson, 1998, for weekly hours worked).
We demonstrate the benefit of extending the NW estimator to regression with singular data by applying it to the data from a randomized experiment in a job training program evaluation study by LaLonde (1986). Following the work by Rosenbaum and Rubin (1983), Dehejia and Wahba (1999, 2002) applied propensity score methods to the LaLonde data for estimation of causal treatment effects in an attempt to generalize the experimental results to nonexperimental data. The propensity score matching was used instead of a multivariate nonparametric model with matching on individual characteristics which was deemed impractical because of the high dimensionality of the regressors. The benefits of kernel regression were analyzed in Heckman et al. (1997, 1998) (without allowing for singularity). The LaLonde data and methodologies were discussed by Angrist and Pischke (2009) and in a recent review by Imbens and Xu (2024). As shown here the discreteness of most of the regressors implies reduced dimension of the support of the joint distribution; for continuous variables existence of continuous density is not imposed and mass at zero in income is accounted for. Our asymptotic results provide the validity of the NW estimator for this singular distribution. The kernel estimators we employ give new insights into the heterogeneous effects of the program, based on a variety of individual characteristics and compare quite well with random forest estimates of the conditional average treatment effect on the treated (CATT) (Wager and Athey, 2018).
Our results also extend to functional regression (see, e.g. Ramsey and Silverman, 2005) where estimation and inference techniques have been developed by Ferraty and Vieu (2004) and pointwise asymptotic normality was established in regression for a Banach or metric space by Ferraty et al. (2007), Ferraty and Vieu (2006), and Geenens (2015) in the i.i.d. case. Masry (2005) derived the limit distribution for a strongly mixing process. Recently Kurisu et al. (2025) made a case for extending the univariate set-up of functional regression by considering jointly a random vector and a function to obtain an estimate for the propensity score used in evaluating the average treatment effect. We establish asymptotic normality in multivariate functional regression with regressors in any number of heterogeneous metric spaces. This provides a basis for simultaneously evaluating the impact of the different predictors rather than comparing their performance in distinct models, as in Caldeira et al. (2020) and Ferraty and Nagy (2022).\footnote{ E.g., Caldeira et al. (2020) compares the model forecasting aggregate stock market excess return on a function representing the history of returns with regression models based on traditional predictors. Ferraty and Nagy (2022) compare the performance of separate models for predicting adult height with functional regressors (one being growth velocity profiles from ages 1-10 and the other for 5-8).} Our simulations show that using multivariate rather than univariate functional regression can improve the fit of the kernel estimator.
We derive asymptotic normality results for a random regressor $X$ supported on some domain in a vector space, $\mathbb{R}^{q},$ or metric, semi-metric space, $\Xi ^{\left[ 1\right] },$ or a product of such spaces $\Xi ^{\lbrack q]}\equiv \Xi _{1}^{[1]}\times \cdots \times \Xi _{q}^{[1]}$. The metrics on $\Xi _{l}^{[1]}$, $\left\Vert .\right\Vert _{l}$, may differ for each of the $q$ components of function spaces, thus as in Kurisu et al. (2025) one may be the $\mathbb{R}^{1}$ space and the other one a function space. A key ingredient in our technical derivations is small cube probability, which characterizes local properties of $X$ in the general multivariate case in place of the density. We introduce this concept here.
In the univariate metric space, $\Xi =\Xi ^{\lbrack 1]}$, the probability measure is characterized by the small ball probability (e.g., see Ferraty and Vieu, 2006): for the ball $B\left( x,h\right) =\left\{ X:\left\Vert x-X\right\Vert \leq h\right\} $ {centered at $x$} in $\Xi ^{\left[ 1\right] } $ {the probability measure} is denoted $P_{X}\left( B\left( x,h\right) \right) $. Characterizing the measure locally via a ball is insufficient when we wish to examine heterogeneous regressors in $\mathbb{R}^{q}$ {or, in general, in product metric spaces $\Xi ^{\left[ q\right] }$.}
Kankanala and Zinde-Walsh (2024) introduced small cube probability for a cuboid. A cuboid $C\left( x,h\right) $ centered around $x=\left( x^{1},...,x^{q}\right) \in \mathbb{R}^{q}$ for a vector $h=\left( h^{1},...h^{q}\right) ^{\prime }$ with positive finite components is defined as the set $C\left( x,h\right) =$ $\left\{ X\in \mathbb{R}^{q}:\left\vert X^{l}-x^{l}\right\vert \leq h^{l},l=1,\cdots ,q\right\} .$ With the distribution function of $X$ given by $F_{X}$ the corresponding probability measure is $P_{X}\left( C\left( x,h\right) \right) =\int_{C\left( x,h\right) }dF_{X}.$ The small cube probability permits us to extend the regression on univariate metric spaces to $\Xi ^{\lbrack q]}$ where the probability measure $P_{X}$ is defined. For the cuboid
the corresponding small cube probability is also denoted $P_{X}\left( C\left( x,h\right) \right) .$
One of our contributions is the derivation of auxiliary technical results that express moments for the multivariate kernels and related functions in terms of the small cube probabilities without appealing to differentiability on which previous multivariate derivations relied. The moments and moment bounds are derived for general multivariate local functions under arbitrary distributions over $\mathbb{R}^{q}$ or probability measures over $\Xi ^{ \left[ q\right] }$. Bounds on a moment functional expressed via power of small cube probability pinpoint the rate of growth of the functional. These results generalize the derivations for the univariate kernel used in functional regression to the multivariate setting. The full details of these auxiliary technical results are presented in the supplemental material (Appendix A). The moment expressions could find use in other contexts, for instance for local linear and local polynomial estimation in $\mathbb{R}^{q}$ or in products of suitable metric spaces, $\Xi ^{\left[ q\right] },$ kernel estimation of distribution functions and conditional distributions in $ \mathbb{R}^{q}$ as well as to kernel regression of objects in metric spaces on objects in products of spaces.
Implementation of the NW estimator relies on a tuning bandwidth parameter. We show that in $\mathbb{R}^{q}$ a popular cross-validation method of choosing a bandwidth with properties that were worked out for the absolutely continuous (a.c.) case, has similar properties in some empirically relevant classes of singular distributions with dimension-reducing singularity. We also examine adaptive bandwidth selection for regressor distributions that are represented by a mixture of a continuous distribution with some mass points.
We provide simulation evidence on some important features of the behavior of the NW estimator under possible singularity of the distribution of regressors, $F_X$, in $\mathbb{R}^{q}$, in particular on the pointwise rate of convergence and specific impact of mass points. We examine the behavior of the NW estimator for models with dependence on both a functional object in $\Xi ^{\left[ 1\right] }$ and a random variable.
The structure of the paper is as follows. Section 2 provides the set-up suitable for the multivariate vector and functional regression highlighting the probability measure for the regressor. Section 3 gives the asymptotic normality results under the most general distributional assumptions. Section 4 discusses implementation, in particular, bandwidth selection. Section 5 provides a sketch of the simulation results and Section 6 is devoted to the empirical study. The supplementary material collects various auxiliary results and the proofs as well as the details of the Monte Carlo simulations and the empirical study.
This section provides the formula for the Nadaraya-Watson (NW) kernel estimator over $\Xi ^{\left[ q\right] }$, introduces some useful notation and gives formal assumptions. The distributional assumptions are very general in that they do not restrict the distribution over $\mathbb{R}^{q}$ to have absolutely continuous components, and apply to the probability measure over the multivariate metric space $\Xi ^{\left[ q\right] }$ for an arbitrary $q.$
In the product space $\Xi ^{\left[ q\right] }$ we define a vector $w$ as $ \left( w^{1},...,w^{q}\right) ^{T}$ where each component is in the corresponding space, thus for $\Xi =\mathbb{R}^{q},$ $w$ is a $q$ -dimensional vector of reals, in $\Xi =\Xi ^{\lbrack q]}$ each $w^{l}\in \Xi _{l}^{[1]},\ l=1,...,q.$ The bandwidth vector is $h=\left( h^{1},...,h^{q}\right) \in \mathbb{R}^{q}$ with $0<\underline{h}=\min \left\{ h^{1},...,h^{q}\right\} >0$ and $\bar{h}=\max \left\{ h^{1},...,h^{q}\right\} $. We use the same notation $\left\Vert \cdot \right\Vert $ for the absolute value of a scalar in $\mathbb{R}^{1}$, the Euclidean norm for a vector in $\mathbb{R}^{q}$ or norm for a function in $ \Xi =\Xi ^{\lbrack 1]},$ with $\Xi ^{\lbrack 1]}$ a Banach space, or metric (semi-metric) in a metric space $\Xi ^{\lbrack 1]}$; where the meaning is not clear from the context we shall specify.
The NW estimator, $\widehat{m}(x)$, for a sample $\left\{ \left( Y_{i},X_{i}\right) \right\} _{i=1}^{n}$ generated by ((ref)) is defined below. Generically the argument of the kernel function is
The NW estimator is given by
where $K(W_{i}(x))=K(W_{X_{i}}\left(x\right) )$ is a multivariate (non-negative) kernel function and $h$ usually depends on $n$; $x$ such that at least for some $i$ we have that $K\left( W_{i}(x)\right) >0$. The kernel function $K$ and bandwidth vector $h$ determine the properties for the NW estimator. In the metric space $\Xi =\Xi ^{\lbrack 1]}$ the kernel function $ K$ is defined for a univariate non-negative argument; in the case $\Xi =\Xi ^{\lbrack q]}$ with $q>1$ different bandwidths could appear for the different components $W_{X}^{l}\left( x\right) ,l=1,..,q.$ With a symmetric kernel on $\mathbb{R}^{q}$ we can just write $W_{X}\left( x\right) =h^{-1}\left\Vert x-X\right\Vert $ for any $\Xi .$
We restrict the multivariate kernel functions on $\mathbb{R}^{q}$ to have bounded support and be suitably differentiable in the interior.
Let $I_{{\Greekmath 0118} }$ denote any subset of the set $\left\{ 1,...,q\right\} $ $\ $ of consecutive non-negative integers; there are $2^{q}$ such subsets including the empty set $\varnothing ;$ denote by $q\left( {\Greekmath 0118} \right) $ the cardinality of the set $I_{{\Greekmath 0118} }=\left\{ j_{1},...,j_{q({{\Greekmath 0118} })}\right\} $ with $j_{1}<...<j_{q({{\Greekmath 0118} })}.$ We use $\prod_{j\in I_{{\Greekmath 0118} }}^{{}}\left( \partial _{j}\right) $ to denote an operator that, when applied to a differentiable function $g\left( z\right) =g\left( z^{1},...,z^{q}\right) $ at $z$, maps it to its partial derivative for $j_{1}<...<j_{q({{\Greekmath 0118} })}$, that is
We call a function $g(z)$ \textquotedblleft sufficiently differentiable\textquotedblright\ if for any set $I_{{\Greekmath 0118} }$ the derivative $ \left( \prod_{j\in I_{{\Greekmath 0118} }}^{{}}\left( \partial _{j}\right) \right) g\left( z\right) $ exists and is continuous at any point on the interior of its support.
The following assumption is made on the kernel function.
Assumptions (ref)(a--c) are satisfied by the commonly employed product kernels of Epanechnikov or quartic kernels. Assumption (ref) (d) is not usual for kernel regression on $\mathbb{R}^{q};$ in the context of univariate functional regression\ it is satisfied by a Type I kernel defined in Ferraty and Vieu (2006) as \thinspace $K:C_{1}I_{\left[ 0,1\right] }\leq K\leq C_{2}I_{\left[ 0,1\right] }$ with some $0<C_{1}\leq C_{2}<\infty .$ Condition (d) in conjunction with (a-c) provides the same type of univariate kernel. Extended to a multivariate setting it can be said that a kernel that satisfies Assumption (ref) (a-d) is a type I kernel. The uniform kernel is an example. The functional regression literature demonstrates that with kernels of type I asymptotic normality can be established in more general univariate settings. As commonly used in $ \mathbb{R}^{q}$ kernels are not of type I, the asymptotic normality results are given separately to apply under Assumption (ref)(a,b,c) and under the full Assumption (ref).
Consider the process $\left\{ \left( X_{i},Y_{i}\right) \right\} _{i\in \mathbb{N}}.$ An i.i.d sequence would provide the simplest characterization, but strong mixing makes it possible to extend the results to time series data. Denote by $\mathcal{F}_{a}^{b}$ the sigma algebra generated by $ \left\{ \left( X_{i},Y_{i}\right) \right\} _{i=a}^{b}.$ Define
Recall that the process is strong mixing if ${\Greekmath 010B} \left( l\right) \rightarrow 0$ as $l\rightarrow \infty .$
The assumption requires a polynomial bound on the rate of decline of the mixing coefficient with a link to the moment of $Y;$ it is similar to those in Masry (2005) and Hong and Linton (2020).
Assumption (ref) requires Holder continuity of $m\left( x\right) ;$ it would follow from differentiability or Lipschitz continuity in $\mathbb{R} ^{q}$ with ${\Greekmath 010E} =1$. In the above assumptions, and below, $L$ and $M$ denote lower and upper bounds of functions where the subscript typically denotes the function whose bounds are provided. The bounds could depend on the point $x.$
For the probability measure $P_{X}$ on a generic space $\Xi ,$ that could coincide with $\mathbb{R}^{q},$ $\Xi ^{\left[ 1\right] },$ or $\Xi ^{\left[ q \right] },$ any point $x\in \Xi $ is a point of support if for $\underline{h} >0$ {the measure $P_{X}\left( C\left( x,\underline{h}\right) \right) >0.$}
By the Lebesgue decomposition, the distribution $F_{X}$ on $\mathbb{R}^{q}$ can be represented as a mixture of an absolutely continuous distribution, $ F^{a.c.},$ a singular distribution (the distribution function is continuous but there is no function that integrates to it), $F^{s},$ and a discrete distribution, $F^{d}:$
In a multivariate setting as soon as at least one variable is continuously distributed, mass points do not arise and the joint distribution is a continuous function, but with some discrete components or mass points in some of the continuous components the distribution can no longer be absolutely continuous and is singular. In many applications at least one of the variables is assumed continuous and in a semiparametric regression often an index model is assumed (single index in Ichimura, 1993; multiple index in Donkers and Schafgans, 2008) to avoid singularity as well as to reduce dimensionality of the model. In a general multivariate distribution the presence of singularity achieves reduction of dimension (see, e.g. examples 2-4 in Kankanala and Zinde-Walsh, 2024) that will have a similar beneficial effect on the convergence of the kernel estimator.
The discussion in this section applies to the space $\mathbb{R}^{q}$ as a special case. Particular classes of probability measures considered in univariate functional regression (e.g. Ferraty et al., 2007, Ferraty and Vieu, 2006) have a small ball probability centered at a point $x$ of support either with a polynomial (fractal) rate of decline $P_{X}(B(x,h))\sim C\left( x\right) h^{{\Greekmath 011C} }>0$ (${\Greekmath 011C} >0),$ or with an exponential type rate of decline $P_{X}(B(x,h))\sim C\left( x\right) \exp (-h^{-{\Greekmath 011C} _{1}}\log h^{-{\Greekmath 011C} _{2}})$ (${\Greekmath 011C} _{1}>0,{\Greekmath 011C} _{2}>0)$ as $h\rightarrow 0.$ This characterization can be applied to the multivariate setting by replacing $ B(x,h)$ with the cuboid and the univariate bandwidth in the rate with $\bar{h }.$ The exponential rate of decay of the small ball probability requires a type I kernel and leads to slow convergence for the estimators (curse of dimensionality). There are ways to mitigate the curse of dimensionality arising from such exponential decay. It is common to apply finite dimensional approximation of these functionals as suggested in Gasser et al. (1998). Indeed, the case where functional data can be accurately approximated in a finite dimensional space is not rare (corresponds to observation of smooth curves with common shape) as noted by Ferraty and Nagy (2022).
Kernels of type I play an important role in establishing pointwise asymptotic normality in the absence of any restrictions on the decline of the small cube measure. For kernels that may not be of type I sufficient conditions on the shrinkage of the probability measure as $h\rightarrow 0$ were proposed in Assumption $H_{3}$ in Ferraty et al. (2007), Ferraty and Vieu (2006) and were referred to in various subsequent papers on functional regression, e.g. Hong and Linton (2020). The assumption below generalizes these conditions to apply to $C\left( x,h\right) $ on $\Xi ^{\left[ q\right] };$ the assumption is both necessary for the conditions to hold (see the supplementary material, Appendix B) and at the same time sufficient for the convergence results.
A polynomial decay condition places a measure into class $\mathcal{D}.$ Indeed if the small cube probability satisfies
{where for some $c$, $H(x),$ $1\leq c<\infty ,$ $0<H\left(x\right) <\infty $ and $\bar{h}=c\underline{h}<H\left( x\right) ,$ $0\leq s\left( x\right) \leq 1$ and $M_{F}\left( x\right) /L_{F}\left( x\right)<B<\infty $ at all points of support $x$, ((ref)) holds with $C_{{\Greekmath 0122} }=B\left( c/{\Greekmath 0122} \right) ^{q}.$}
Condition ((ref)) applies quite widely and holds for many distributions of regressors used in econometric models. In $\mathbb{R}^{q}$ it is satisfied by any absolutely continuous distribution with\ a positive bounded density function $f_X\left( x\right) $ where $M_{P}\left( x\right) \geq \underset{\tilde{x}\in C\left( x,H\right) }{\sup }f_X\left( \tilde{x} \right) ;$ $L_{F}\left( x\right) =\underset{\tilde{x}\in C\left( x,H/c\right) }{\inf }f_X\left( \tilde{x}\right) $ and $s\left( x\right) =1.$ If $x$ is an isolated mass point then ((ref)) applies with $s=0.$ If $X$ has a linear structure with $r$ common factors, the probability measure is singular and satisfies ((ref)) with $s\left( x\right) =s=\frac{q}{r}.$ For a fractal distribution that is singular with constant $ s,$ $0<s<1$, the bounds also apply.
Condition ((ref)) is satisfied by the general class of Ahlfors (1966) regular (A-r) distributions common in statistics, where for this class $s\left( x\right) =s,$ and $L_{P}\left( x\right) =L$ and $M_{P}\left( x\right) =M$ are constants, as well as by a finite mixture of such distributions (as proved in the supplementary material, Appendix B). Thus an absolutely continuous distribution ($s=1)$ or, more generally, a measure given by a continuous possibly singular distribution function that satisfies ((ref)) contaminated with some mass points ($s=0$) is in $\mathcal{D}$ $;$ this applies to the empirical example examined here, ensuring the pointwise asymptotic normality of the NW estimator with standard kernels.
Consider the product space $\Xi ^{\left[ 2q\right] }=\Xi ^{\left[ q\right] }\times \Xi ^{\left[ q\right] };$ the measure on this product space has marginals $P_{X}$ on each $\Xi ^{\left[ q\right] }$ (see, e.g., Pollard, 2001). {The joint measure $P_{s,t}\left( C\left( x,h\right) \times C\left( x,h\right) \right) ,$ defined as $\Pr \left( X_{t}\in C\left( x,h\right) ,X_{s}\in C\left( x,h\right) \right) ,$ is a product of the measures of the cuboid in the case of independency.} With dependence an additional assumption is made on how the joint measure relates to the small cuboid measure. {We provide the same assumption as in e.g. Masry (2005) and Hong and Linton (2020) for the small cube probability.}
Consider the NW estimator as given by ((ref)), ((ref)). As the sample size increases the bandwidths {are assumed} go to zero. For $\Xi = \mathbb{R}^{q}$ the denominator, $B_{n}(x)$ is proportional to the usual kernel density estimator, given by $h^{-q}B_{n}(x),$ at point $x.$ {When the density, $f_{X}\left( x\right) $, exists and is continuous}, the estimator $ h^{-q}B_{n}(x)$ consistently estimates $f_{X}\left( x\right) ,$ but if the density does not exist, $h^{-q}B_{n}\left( x\right) $ diverges to infinity. Consistency of the NW estimator $\widehat{m}\left( x\right) $ over a univariate metric space was established in Gy\"{o}rfi et al. (2002), the limit distribution in Ferraty et al. (2007), Masry (2005) and Geenens (2015).
The key to the asymptotic normality result is the derivation of the moments for multivariate functions of the form $g(X)K^{m}\left( h^{-1}\left\Vert x-X\right\Vert \right) $ for general probability measures and establishing lower and upper bounds (derivations in the supplementary material, Appendix B). The bounds provide expressions in terms of the small cube probability:
with constants $L_{EgK^{m}}\left( x\right) $ and $M_{EgK^{m}}\left( x\right) $ at $x.$ Most important, ((ref)) provides a lower bound on $ EB_{n}\left( x\right) =EK\left( K^{m}\left( h^{-1}\left\Vert x-X\right\Vert \right) \right) $, given by $L_{EK}P_{X}\left( C\left( x,h\right) \right) $, with appropriate conditions for $L_{EK}$ to be strictly positive to ensure that the denominator of the NW estimator is such that it exists and the limit does not blow up. Type I kernel automatically entails that $L_{EK}>0$, but for kernels such as Epanechnikov the bound requires Assumption {(ref). }With $g\left( X\right) $ that is continuous at $x$
These moment expressions for distributions over $\mathbb{R}^{q}$ hold under the standard assumptions of existence and continuity of (bounded) density $ f_X,$ and the function $g$, {where}
with more details about the $o\left( 1\right) $ term under smoothness of $f_X $ (see, e.g. derivations in Li and Racine, 2007). Once the moments and the bounds are derived, the proofs of asymptotic normality proceed along similar lines to those in Masry (2005).
The point-wise limit normality is provided in the next theorem under two alternative types of conditions: (i) with type I kernel without imposing further constraints on $F_X$, and (ii) not imposing the type I kernel but with the distributional Assumption (ref). Denote the bias of the estimator given $x$, $E\left( \widehat{m}\left( x\right) \right) -m\left( x\right) ,$ by $bias\left( \widehat{m}\left( x\right) \right) .$ The difference $\widehat{m}\left( x\right) -m\left( x\right) $ is delivered by $ \frac{A_{n}^{c}(x)}{B_{n}(x)}$ with the \textquotedblleft centered\textquotedblright\ $A_{n}^{c}(x)=A_{n}(x)-m\left( x\right) B_{n}(x). $
\setcounter{theorem}{0}
Remarks.
The convergence rate in (c) of Theorem (ref) is $O\left( (nP_{X}\left( C\left( x,h\right) \right) ^{-1/2}\right).$\footnote{ This convergence rate obtains under $h\rightarrow 0.$ In the presence of an irrelevant regressor, say $x^{(2)}$, such that $m\left( x\right) =m\left( x^{\left( 1\right) }\right) $ for all $x=\left( x^{\left( 1\right) },x^{\left( 2\right) }\right) ,$ this requirement can be restricted to the function $m\left( x^{\left( 1\right) }\right) $ with the irrelevant $ x^{\left( 2\right) }$ eliminated. For the estimator this elimination can be achieved by setting the bandwidth on components of $x^{\left( 2\right) }$ to be larger than the range of those variables, possibly infinite.} Existence of a limit variance ${\Greekmath 011B} _{\widehat{m}\left( x\right) }^{2}$ requires that $\left( nP_{X}\left( C\left( x,h\right) \right) \right)\frac{\left[ EK\left( h^{-1}\left\Vert x-X\right\Vert \right) \right]^2}{{{\Greekmath 0116} _{2}\left( x\right) E\left[ K^{2}\left( h^{-1}\left\Vert x-X\right\Vert\right)\right] }} $ converges. Without additional assumptions it is possible that the ratio does not converge; see example in the supplementary material (Appendix B) that provides a case when convergence does not hold; this happens when the small cube probability declines very rapidly and the kernel is not uniform. Suitable additional assumptions on the distribution, such as $H_{3}$ in Ferraty et al. (2007) and Condition 3(i) in Masry (2005) and similar ones in subsequent papers provide restrictions on the probability measure on $\Xi ^{\left[ 1 \right] }$ that are sufficient for the convergence. Generally, one needs to ensure that the limits given below on the expectation of the kernel function and its square hold.\footnote{ This implies that the extra condition is also required for the Corollary 1 of Hong and Linton (2020).}
This assumption holds quite widely. From the moment expressions it can easily be shown that it holds for the uniform kernel without any additional distributional assumptions. In the case of continuous density it holds by virtue of ((ref)) with
{Suppose that singularity arises, because of combining discrete and continuous variables in $\mathbb{R}^{q}$ or functional dependence between the regressors, that restrict the support of the distribution to be in some subspace of dimension $r<q$, $V\left(r\right) \subset $ $\mathbb{R}^{q}.$} If the distribution on $V\left( r\right) $ is absolutely continuous with a continuous density, then derivations provide similar limits to ((ref)) with integration over $V\left( r\right) $ and density restricted to $ V\left( r\right) .$
Define now
This limit extends the results that were obtained in the literature on kernel estimation in $\mathbb{R}^{q}$ under smoothness assumptions on the distribution $F_{X}.$ {For functional regression our assumptions are comparable to those of Ferraty et al. (2007), Masry (2005), and subsequent papers while they make the extension to multivariate functional regression possible. }
Estimation of $m\left( x\right) $ requires a selection of the kernel, $K,$ and bandwidth, $h.$ As may be clear from the results here and the literature, type I kernel (such as the uniform) is preferred but other kernels can also deliver asymptotic rates provided the small cube probability does not decline exponentially fast. Aside from the estimator of the conditional mean, estimators of variance and mean squared error are needed to evaluate the performance of the estimator. While in the literature on kernel regression on $\mathbb{R}^{q},$ the leading term of the limit variance is expressed via the density function, often in the actual implementation the corresponding estimators do not make use of plug-in expressions, instead estimating the variance directly from the data and possibly with bootstrap (see Hall and Horowitz, 2013).
Cross-validation procedures in popular statistical packages (such as R) provide a single bandwidth (vector) that was shown to be consistent for the \textquotedblleft optimal\textquotedblright\ bandwidth: minimizer of weighted integrated mean squared error, WIMSE, (e.g. Li and Racine, 2007). The proofs of consistency relied on absolute continuity of the regressors. The consistency results extend to some classes of singular distributions.
WIMSE is defined for an absolutely continuous distribution with density function $f_X\left( x\right) $ as
with some weighting function $M\left( x\right) $ chosen to mitigate boundary effects. The expression can be written with $dF_X$ replacing $f_X\left( x\right) dx$ (valid in the case of singularity):
This function depends on the bandwidth vector $h$ used in the estimator (see the review of bandwidth selection methods, including cross-validation and plug-in in K\"{o}hler et al., 2014). The \textquotedblleft optimal\textquotedblright\ bandwidth vector $h^{0}$ is a minimizer of the WIMSE criterion function based on a trade-off between the variance and bias of the NW estimator.
In the cross-validation procedure the finite sample analogue of WIMSE replaces the expectation by
employing the leave-one-out kernel estimator, $\widehat{m}_{-i},$ and provides the bandwidth vector $h_{cv}$ by minimizing the CV criterion.
Hall et al. (2007) gave a general result about consistency of the cross-validated bandwidth for regression over $\mathbb{R}^{q}$ with discrete and continuous regressors, with some of the regressors possibly being irrelevant. Their general result in Theorem 2.1 was obtained under a set of assumptions that required independent identically distributed observations, restrictions on the support of the probability measure, two continuous derivatives for density, the regression function, and the conditional variance of the error; in addition, for the $d$ continuous relevant regressors $h^{o}=n^{-\frac{1}{4+rd}}a^{o}$ holds with the vector $a^{o}$ having unique, positive and finite components. This result was extended to weakly dependent data by Li et al. (2009) under assumptions that replaced the i.i.d. assumption by requiring strict stationarity and ${\Greekmath 010C} -$mixing in the process for $\left\{ x,y\right\} $ and martingale difference error, with suitable restrictions on the mixing parameters.
The result on the cross-validated bandwidth applies more widely. For instance, consider a singular distribution of $X\in \mathbb{R}^{q}\,$\ where there is a functional dependence among the continuous variables in the presence of possibly some discrete covariates such that the support of the distribution is restricted to a subspace $V\left( r\right) \subset \mathbb{R}^{q}$ of dimension $r<q$ represented by a union of affine subspaces. If, restricted to $V\left( r\right) ,$ the distribution function is such that the conditions of Theorem 2.1 of Hall et al. (2007) or Theorem 1 of Li et al. (2009) are satisfied (Assumption CV) then the conclusions of those theorems are valid and the consistency of the bandwidth and automatic dimension reduction by smoothing out irrelevant regressors hold for this singular distribution. More details are provided in the supplementary material (Appendix B).
Importantly, no knowledge of $V\left( r\right) $ or $r$ is required. This implies that for functionally dependent continuous regressors the knowledge of the number of factors is not required for the consistency of the cross-validated bandwidth or the automatic dimension reduction. We conjecture that in many other cases with possible singularity the cross-validation procedure will facilitate dimension reduction by smoothing out irrelevant variables.
Bandwidth selection could benefit from adaptation to different types of singularity. The treatment of adaptive bandwidth selection in the literature (Fan and Gijbels, 1996, Sain, 1994, Demir et al., 2010) typically focuses on adjusting the smoothing parameter to accommodate the varying data density, but not dealing with singularity or mass points. Adaptive bandwidths can provide a better fit of the criterion function by increasing the number of observations used to estimate the function at a point of sparsity.\footnote{ Given some initial bandwidth $\tilde{h}$ and density estimate at this bandwidth, $\hat{f}_{X},$ an adaptive bandwidth is defined for each point as $h\left( X_{i}\right) =\tilde{h}\left( \frac{\hat{f}_{X}\left( X_{i}\right) }{G}\right) ^{-{\Greekmath 010B} }$ where $G=\left( \prod \hat{f}_{X}\left( X_{j}\right) \right) ^{1/n}$ is the geometric mean of the densities and $ {\Greekmath 010B} $ is typically selected to be $1/2.$ {One could construct $\hat{f} _{X}(x)$ $\ $\ with a uniform kernel in which case it is identical to an estimate of $P(C(x,\tilde{h}))$ by the proportion of observations in the $ \tilde{h}$ cuboid around $x.$}} Such bandwidths can similarly be constructed for cases of singular distributions. But these adaptation procedures still need to be investigated in the case of general mixtures of singular distributions. However, singularity adaptation simplifies considerably for the empirically important case of a mixture of an absolutely continuous distribution with mass points, where the two levels of singularity can be separated. The approach is detailed in the supplementary material (Appendix B).
This section provides the highlights of various simulations that show features of the finite sample performance of the NW estimator under singularity. Additional details and features are in the supplementary material (Appendix C).
In this example we consider the regression distribution with mass points. Alongside we examine the trinormal mixture considered in Kotlyarova et al. (2016), an a.c. distribution which represents features (high density derivatives) that makes it comparable to a singular distribution.
The distribution with mass points, following Jun and Song (2019), is given by
where $F^{d}$ is the discrete uniform distribution function with $ D=\{-1,0,1\}$ the set of mass points; $\Phi $ is the standard Gaussian distribution function.
We simulated 500 random samples $\{(Y_{i},X_{i})\}_{i=1}^{n}$ using the model
for different sample sizes. The error $\{{\Greekmath 0122} _{i}\}_{i=1}^{n}$ is drawn independently of the regressor and has a standard Gaussian distribution; ${\Greekmath 011B} $ is selected to yield a given signal to noise ratio, $ snr$, here selected to equal one. We use the Epanechnikov kernel $K(u)=\frac{ 3}{4}(1-{u^{2}})1(u^{2}\leq 1)$ and obtain the leave-one-out cross-validated bandwidth.
We analyze the pointwise RMSE at a coarse grid of points across samples of size $n$ equal to 50, 100, 200, 400, 800, 1600, 3200 based on 500 replications from the above DGP. To obtain empirical rates of convergence we regress $\log (RMSE)$ on $\log (n)$ and a constant. The coefficient on $\log (n)$ is the \textquotedblleft realized\textquotedblright\ rate of convergence; for example if $RMSE\propto n^{-2/5}$ (univariate kernel regression with smooth density and second order kernel) then $\log (RMSE)={\Greekmath 010B} _{0}+{\Greekmath 010B} _{1}\log (n)$ and ${\Greekmath 010B} _{1}$ should be close to -0.4.\footnote{ The authors thank Jeff Racine for suggesting this insightful exercise. See also Hall and Racine (2015).}
In Table 1, illustrative results are provided for the regressor distribution with mass points and the trinormal distribution on a set of support points.
For the distribution with mass points, the NW estimator with cross-validated bandwidth performs remarkably well at points sufficiently far from our mass points (faster than the expected rate of -0.4). The empirical rate at mass points is close to -0.5 when the bandwidth is set equal to zero. The empirical convergence rate is slow for points close to the mass points (within the small ball probability measure under cross validated bandwidth) due to the boundary weight associated with mass in the neighborhood. { Bandwidth adaptive to masspoints improves the rate.} The convergence rates for the trinormal distribution, are reflective of usual smooth nonparametric regression although are somewhat faster at points with high derivatives.
We consider a model where $m(X)=\log (X_{1})+\log (X_{2})$ with regressors $ X_{1}$ and $X_{2}$ satisfying $X_{1}+X_{2}=d(k),$ with fixed $d(k)$ corresponding to $k=1,2,3$.\footnote{ An example could be where $X_{1}$ and $X_{2}$ represent earnings of the husband and wife and, for tax purposes, their combined income is set at some $d(k)$.} This is equivalent to a model with one continuous and one discrete regressor $m(X)=\log (X_{1})+\log (D-X_{1}),$ with $D=d(k)$.
We simulated 500 random samples $\{(Y_{i},X_{1i},X_{2i}\}_{i=1}^{n}$ using the model
for different sample sizes with the additive error chosen as in the previous simulation. The probability of an observation belonging to a sub-population with $k=1,2,3$ is set equal to $0.5$, $0.3$, and $0.2$ respectively and $ d(1)=4,d(2)=6,d(3)=7$; $X_{1}$ is drawn from the uniform distribution: $ U[1,3]$.
We implement the NW estimator first using $X_{1}$ and $X_{2}$ as regressors (NW.c) and second using $X_{1}$ and $D$ as regressors (NW.d) and obtain the leave-one-out cross-validated bandwidths. For the discrete regressor $D$ we use special discrete kernel weights proposed by Wang and van Ryzin (1981) in accordance with Racine and Li (2004).\
In Table 2 we provide illustrative results comparing the empirical rate of convergence of the NW.c and NW.d at a grid of points.
The reduced dimensionality is reflected in the estimates of the pointwise rate of convergence which are around $-0.40$ rather than the slower rate of $ -0.33$ the presence of two continuous regressors would suggest ($q=2$). The estimate of the empirical rate for NW.c is slightly faster than NW.d, moreover, indicating that there is no gain from separate treatment of discrete regressors. With the reduced dimension structure here therefore one gets the rate corresponding to the Hausdorf dimension of the regressor space automatically without the need to recognize that it is possible to transform the regressors to one discrete, and one continuous variable.
Here we examine a functional regressor in a multivariate setting. Consider a bivariate conditional mean function $m(X)=m(X_{1},X_{2})$, where $X_{1}$ is a functional regressor and $X_{2}\in \mathbb{R}$ \ may be correlated with some $m_{1}(X_{1})$. Let
Following Ferraty et al. (2007), the functional regressor is defined as
with $a_i$ and $b_i$ drawn from $U(-1,1)$, $w_i$ drawn from $U(-{\Greekmath 0119} ,{\Greekmath 0119} )$ and
For $X_{2}$ we consider two possibilities: (a) a N(0,1) random variable independent of $X_{1}$; (b) $X_{2}=m_{1}(Z)$ where $Z(t)$ is a functional regressor similar to $X_{1}(t)$ with $(a_i,b_i,w_i)$ replaced by $ (a_i^{\prime },b_i^{\prime },w_i^{\prime })$ where the correlation between $ (a_i^{\prime },b_i^{\prime },w_i^{\prime })$ and $(a_i,b_i,w_i)$ is given by ${\Greekmath 011A} $ (and set equal to either $0$ or $0.8$).
For the functional regressor $X_{1}$ we use the same metric as in Ferraty et al. (2007), that is $\left\Vert x_{1}-X_{1}\right\Vert _{1}=\sqrt{ \int_{-1}^{1}\left( x_{1}^{\prime }(t)-X_{1}^{\prime }(t)\right) ^{2}dt}.$ We use a product kernel with kernel $K(u)=1-u^{2}$ defined on $[0,1]$ for the functional regressor and the Epanechnikov kernel defined on $[-1,1]$ for $X_{2}$.
Table 3 shows RMSE of the NW estimator at the cross-validated bandwidths as well as RMSE where either the functional or scalar regressor is dropped. The loss from misspecifying the functional regression as univariate can be substantial.
The causal inference literature has made extensive use of the LaLonde (1986) data on the National Supported Work Demonstration (NSW) program following the release of that data by Dehejia and Wahba (1999, 2002). Their finding, that propensity score-based methods provide a way to generalize the experimental results on the impact of training to nonexperimental data, was influential and led to significant methodological advances and practical changes as discussed in the review by Imbens and Xu (2024). Here, we consider the experimental sample to analyze potential heterogeneous treatment effects using multivariate kernel estimation. Kernel-based matching on individual characteristics, advocated in Heckman et al. (1997, 1998), was not considered due to the claimed high dimensionality of the regressors. We show that it is both feasible and insightful for this data due to the dimension reduction implied by the presence of several discrete, discretized and categorical regressors. The mass point of the regressor on pre-treatment earnings at zero further contributes to the regressor singularity and we do not require continuity or indeed existence of density over positive values, thus kinks or mass at positive values are not excluded. The kernel estimator is applicable to such singular distributions.
We focus here on the full LaLonde NSW male sample which contains 297 treated individuals and 425 controls where the pre-intervention variables are well-matched.\footnote{ In the sub-sample with 1974 earnings data in Dehejia and Wahba (1999) the distribution of 1975 earnings exhibits a significantly different mass at zero between the treated (68%) and untreated (60%).}
With $Y$ denoting the post-treatment outcome, $T$ the treatment and $X$ the individual pre-treatment characteristic(s), we use the nonparametric regression model
to evaluate the heterogeneous effects of the treatment as
The heterogeneous effect of treatment on the treated, also known as the conditional average treatment effect, CATT, is given by
Focusing on the latter, we use the NW estimates to evaluate
for all treated individuals $n_T$ (i.e., we use both the actual and the counterfactual treatment for our estimates).
First, we consider a bivariate kernel regression model where we only use the pre-treatment earnings (re75) as regressor $X$. Following that, we estimate the multivariate model with the full set of variables $X$, where in addition to the pre-treatment earnings we include years of education, high school “no degree” status, race, age, marital status, and pre-treatment unemployment status, u75. It is not unreasonable to attempt nonparametric estimation for this problem where the only truly continuous regressor is earnings (and possibly age and education) as singularity provides dimension reduction.
As with our simulations, we use the np package in R for the nonparametric estimation where we consider the Epanechnikov (e), Uniform (u) and discrete (d) kernel.\footnote{ We use the discrete kernel proposed by Aitchison and Aitkin, 1976, where $ K((d-d_{i})/h)=1-h$ if $d=d_{i}$, else $h$ where $h\in \lbrack 0,1/2]$.} Bandwidth selection is based on cross validation and we consider the adaptive bandwidth selection approach that accounts for the masspoint. As was shown in our simulations the rate improvement associated with singularities does not require special attention to discrete variables to benefit from it.
Estimation results are reported in detail in the supplementary material (Appendix D). Below the main findings are summarized.
For the bivariate regression model, the cross validated bandwidths confirm that we should not smooth across treated and untreated observations and that local heterogeneous treatment effects as related to pre-treatment earnings are present. Figure 1, displays estimates of the conditional expectation using the Epanechnikov kernel by treatment status and pre-treatment earnings together with the bootstrapped confidence bounds. It suggests that treatment for individuals at low levels of pre-treatment earnings, in particular, is beneficial.
The adaptive bandwidth results in a slightly better in-sample correlation between the post-treatment outcome, $re78$, and its fit (increasing from 0.2098 to 0.2116 (for OLS the correlation is 0.1697)); bandwidths obtained using non-masspoint-observations only are quite similar to those obtained when including the masspoints in this case. The NW estimates with the adaptive bandwidth provide values of CATT that on average equal $\$920$ (76), \$920 (76), and \$906 (80) (standard error in brackets) for the (e,e), (d,e), (d,u) kernels on $(T,X)$, respectively.\footnote{ As discussed in the supplemental material (Appendix D), we denote the kernel with two arguments: the first argument denotes the kernel applied to all binary regressors (treat, u75, nodegree, black, hispanic, and married) and the second argument denotes the kernel applied to the other regressors (re75, educ, and age).} For comparison, the local linear kernel based estimates on average equal \$822 (49) with the (e,e) kernel, while the average of the CATT estimates based on random forest (RF) equal \$848 (52). The CATT results of the kernel regression based approach are more variable than those obtained using the random forest approach. For observations at mass points, CATT estimates using adaptive bandwidth are closer to those obtained using the random forest based approach.
For the multivariate model the cross-validated bandwidths provide important insights. Firstly, even though pre-treatment earnings is still relevant, the bandwidth is much larger than in the baseline model for all kernels, suggesting a reduction of the heterogeneous impact with individual's pre-treatment earnings. The bandwidths selected for nodegree, hispanic and married are large, signaling that these variables are not relevant (these regressors are automatically smoothed out from the regression function). At the same time, the bandwidths for education and age imply a heterogeneous impact associated with those characteristics, although the size of the bandwidth for age is fairly large.
The inclusion of additional controls yields an improvement in the in-sample correlation between the post-treatment outcome and its fit. For the (e,e) kernel we see an increase in correlation from 0.210 in the bivariate model to 0.338 (for comparison, for OLS the correlation equals 0.209 when age squared is included as well); the results for the (d,e) and (d,u) kernel are comparable.
To highlight the heterogeneity of the treatment effect of education and its interplay with race, we display in Figure 2 estimates of the conditional expectation by treatment status, years of education, and race for an individual with median age and pre-treatment earnings together with the bootstrapped confidence bounds.
The graph reflects a heterogeneity of the impact of treatment whereby the more educated individuals identified as black appear to benefit more from treatment than their nonblack counterparts. Gains of treatment arise where the confidence band around the estimated nonparametric fit $\widehat{m}(x,1)$ lies above that of $\widehat{m}(x,0)$; for non-black individuals this is at the middle range of education, for black individuals this starts around 10 years of education and is rising over that range. These results are further supported when evaluating the average CATT for black individuals across different levels of education (see supplemental material, Appendix D).
Box-plots of the CATT estimates for the multivariate model using the NW regression estimate and the RF estimates are presented in Figure 3. The limit distributional results of Wager and Athey (2018) do not apply here as many components of $X$ are not continuously distributed.
The kernel based regression CATT results remain more variable than those provided by the random forest approach, but their interquartile range is comparable. The NW kernel based estimates of the CATT on average exceed the RF based estimates: $\$1,045$ (107), \$1,018 (108), and \$1,019 (104) for the (e,e), (d,e) and (d,u) kernel on $(T,X)$ against \$794 (54) based on the random forest.
The NW based results are stable across kernel, give interpretable insights and with cross-validation make it possible to detect irrelevant regressors.
Acknowledgements: The authors thank the participants at the Econometric Study Group conference in Bristol, the Canadian Econometric Study Group conference, and Saraswata Chaudhuri for their comments. We thank Jeffrey Racine for his discussion and valuable suggestions at the CESG 2023 and Sid Kankanala for insightful comments on earlier versions of the paper. We thank the Associate Editor and three anonymous referees for their careful reading of the paper and very helpful comments and suggestions.
Funding: Victoria Zinde-Walsh gratefully acknowledges financial support from the Natural Sciences and Engineering Research Council of Canada (NSERC) grant 253139.