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.
69,900 characters · 14 sections · 12 citation commands
lpdensity: Local Polynomial Density Estimation and Inference
Nonparametric estimation of a probability density function (PDF), as well as its associated cumulative distribution function (CDF) or higher-order derivatives thereof, plays an important role in empirical work across many disciplines. Sometimes these quantities are the main objects of interest, while in other cases they are useful ingredients in forming more complex nonparametric or semiparametric statistical procedures. See Wand-Jones_1995_Book and Fan-Gijbels_1996_Book for classical textbook introductions to kernel-based density and local polynomial methods.
This article discusses the main methodological and numerical features of the software package \pkg{lpdensity}, available in both \proglang{R} and \proglang{Stata}, which implements the local polynomial smoothing approach proposed and studied in \citet*{Cattaneo-Jansson-Ma_2020_JASA,Cattaneo-Jansson-Ma_2021_LocRegDistribution} for estimation of and inference on a smooth CDF, PDF, and derivatives thereof. In a nutshell, the idea underlying this estimation approach is to first approximate the discontinuous empirical CDF using local polynomial methods, and then employ that smoothed approximation to construct estimators of the distribution function, density function, and higher-order derivatives.
The resulting local polynomial density estimator is intuitive and easy to implement, and exhibits several interesting theoretical and practical features. For example, it does not require pre-binning or any other complicated pre-processing of the data, and enjoys all of the celebrated features associated with local polynomial regression estimation Fan-Gijbels_1996_Book. In particular, it automatically adapts to the (possibly unknown) boundaries of the density's support, a feature that is unavailable for most other density estimators in the literature. See Karunamuni-Albert_2005_SM for a review on this topic. Two exceptions are the local polynomial density estimators of Cheng-Fan-Marron_1997_AoS and Zhang-Karunamuni_1998_JSPI, which require pre-binning of the data or, more generally, pre-estimation of the density near the boundary, thereby introducing additional tuning parameters that need to be chosen for implementation. In contrast, the density estimator implemented in the \pkg{lpdensity} package requires choosing only one tuning parameter: the bandwidth entering the local polynomial approximation, for which the package also offers data-driven selectors. Furthermore, following the results in \citet*{Calonico-Cattaneo-Farrell_2018_JASA,Calonico-Cattaneo-Farrell_2020_CEOptimal}, robust bias-corrected inference methods are also implemented, which allow using mean squared error (MSE) optimal or the integrated mean squared error (IMSE) optimal bandwidth choices when forming confidence intervals or conducting hypothesis testing.
The software implementation covers smooth estimation of the distribution and density function, and derivatives thereof, for any polynomial order at both interior and boundary points. \citet*{Cattaneo-Jansson-Ma_2020_JASA,Cattaneo-Jansson-Ma_2021_LocRegDistribution} give formal large-sample statistical results for these estimators, including (i) asymptotic expansions of the leading bias and variance, (ii) asymptotic pointwise and uniform normal approximations, (iii) consistent standard error estimators, (iv) consistent data-driven bandwidth selection based on asymptotic MSE expansions of the point estimators, and (v) asymptotically valid uniform inference and confidence bands. Importantly, all these results apply to both interior and boundary points simultaneously. We briefly summarize these results in the upcoming sections, and illustrate them numerically, including a comparison with other methods available in \proglang{R}.
In the remaining of this article, we focus on the \proglang{R} implementation of the software package \pkg{lpdensity} (\url{https://CRAN.R-project.org/package=lpdensity}), but all functionalities are also available in \proglang{Stata}. See Appendix (ref) for more details. The \proglang{R} package includes the following two main functions.
In addition, the methods \code{coef()}, \code{confint()}, \code{plot()}, \code{print()}, \code{summary()} and \code{vcov()} are supported for objects returned by \code{lpdensity()}, and the methods \code{coef()}, \code{print()} and \code{summary()} are supported for objects returned by \code{lpbwdensity()}. In particular, the function \code{plot()}, building on the \pkg{ggplot2} package in \proglang{R}, can be used to plot the estimated CDF, PDF, or higher-order derivatives for graphical illustration. This function takes the output from \code{lpdensity()}, and plots both point estimates and confidence intervals/bands for a collection of grid points.
There are several other packages and functions available for kernel-based density estimation in \proglang{R}. Table (ref) gives a summary of their functionalities. As shown in that table, the package \pkg{lpdensity} is the first to offer consistent estimation of the CDF, PDF and density derivatives for both interior and boundary points, higher-order bias reduction, and valid inference both pointwise (confidence intervals) and uniformly (confidence bands). Section (ref) compares the numerical performance of these packages in a simulation study.
This article continues as follows. Section (ref) provides a brief, self-contained overview of the main ideas underlying the local polynomial estimators implemented in the package \pkg{lpdensity}. Section (ref) illustrates the main features of our package. Section (ref) showcases its finite-sample performance and compares it with other \proglang{R} packages implementing kernel-based density estimators. Section (ref) concludes. We also include two appendices. Appendix (ref) discusses in more detail our data-driven bandwidth selectors, and Appendix (ref) illustrates the \proglang{Stata} implementation of the \pkg{lpdensity} package. Installation details, scripts replicating the numerical results reported herein, links to software repositories, and other companion information, can be found in the package's website: \url{https://nppackages.github.io/lpdensity/}.
This section offers a brief overview of the main methods implemented in the \proglang{R} and \proglang{Stata} package \pkg{lpdensity}. For formal results, including assumptions, proofs and any other technical details see \citet*[hereafter CJM]{Cattaneo-Jansson-Ma_2020_JASA,Cattaneo-Jansson-Ma_2021_LocRegDistribution}.
We assume that $X_1,X_2,\cdots,X_n$ is a random sample from the random variable $X\in\mathcal{X}$, where $F(x)$ denotes its smooth CDF, $f(x)$ denotes its smooth PDF, and $\mathcal{X}\subseteq\mathbb{R}$ denotes its (possibly restricted) support, which can be bounded or unbounded. As it is well known, conventional kernel density estimators will be biased at or near boundary points, and other density estimators must be used if the goal is to estimate a density function on a compact support Karunamuni-Albert_2005_SM. The package \pkg{lpdensity} implements a simple, easy-to-interpret and boundary adaptive density estimator based on local polynomial methods. As a by-product, the package also offers a smooth local polynomial estimate of the CDF as well as density derivatives. To cover all cases in an unified way, we employ the notation $g^{(\nu)}(x)=\partial^{\nu}g(x)/\partial x^{\nu}$ and $g(x)=g^{(0)}(x)$ for any smooth function $g(\cdot)$, and define $F(x)=F^{(0)}(x)$, $f(x)=F^{(1)}(x)$, and derivatives of the density function as $f^{(\nu-1)}(x)=F^{(\nu)}(x)$ with $\nu=1,2,\cdots$, with $f(x)=f^{(0)}(x)$.
To describe the estimators implemented in the package \pkg{lpdensity}, consider first the weighted empirical distribution function \[ \hat{F}(x)=\frac{1}{n}\sum_{j=1}^{n}W_j\mathds{1}(X_{j}\leq x), \] where $\mathds{1}(\cdot)$ is the indicator function, and the weights $W_j$ are introduced for empirical applications such as missing data or counterfactual comparison (see CJM for examples). We assume these weights are normalized so that $\sum_{j=1}^n W_j/n=1$. The package \pkg{lpdensity} allows for a possibly estimated weighting scheme embedded in $\hat{F}(x)$, although for simplicity we will assume that each $W_i=1$ throughout this article. That is, $\hat{F}(x)$ is taken to be the standard root-$n$ consistent empirical distribution function estimator of $F(x)$.
As an alternative to conventional kernel density estimators, consider an estimator that first smooths out $\hat{F}(x)$ using local polynomials, and then constructs an estimator of $f(x)$ (and its derivatives). For $x\in\mathcal{X}$, the estimator implemented in the \pkg{lpdensity} package is
where $\mathbf{r}_{p}(u)=(1,u,u^{2},\cdots ,u^{p})^\mathsf{T}$ is the $p$-th order polynomial expansion, $K(\cdot)$ is a kernel function such as the uniform or triangular kernel, and $h$ is a positive vanishing bandwidth sequence. The estimator approximates the discontinuous empirical CDF $\hat{F}(x)$ by a smooth local polynomial expansion using the weighting scheme implied by the kernel function, localized around the evaluation point $x$ according to the bandwidth $h$. CJM showed that \[\widehat{\boldsymbol{\beta}}_{p}(x) \overset{\mathbb{P}}{\to} \boldsymbol{\beta}_p(x) =
^\mathsf{T}, \] as $h\to0$ and $nh^{2p-1}\to\infty$, where $\overset{\mathbb{P}}{\to}$ denotes convergence in probability. This implies that the least squares coefficients $\widehat{\boldsymbol{\beta}}_{p}(x)$ are consistent estimators of the CDF, PDF, and derivatives thereof at the evaluation point $x$.
Therefore, the generic local polynomial distribution estimator takes the form: \[\hat{F}_{p}^{(\nu)}(x) = \hat{f}_{p}^{(\nu-1)}(x) = \nu! \mathbf{e}_{\nu}^\mathsf{T}\widehat{\boldsymbol{\beta}}_{p}(x), \qquad 0\leq\nu\leq p, \] where $\mathbf{e}_{\nu}$ denotes the conformable unit vector that extracts the $(\nu +1)$-th estimated coefficient. This estimator is implemented in the function \code{lpdensity()}, given a choice of evaluation point $x$, polynomial degree $p$, derivative order $\nu$, kernel function $K$, and bandwidth $h$. In particular, the local polynomial density estimator is $\hat{f}_{p}(x)=\hat{F}_{p}^{(1)}(x)=\mathbf{e}_1^\mathsf{T}\widehat{\boldsymbol{\beta}}_{p}(x)$, which is implemented via the default \code{lpdensity(...,p=2,v=1)}, employing a quadratic approximation to the empirical distribution function to construct the density estimator $\hat{f}_{2}(x)$. Similarly, higher-order derivatives of the CDF can be estimated through $\hat{f}^{(\nu-1)}_{p}(x)=\hat{F}^{(\nu)}_{p}(x)=\nu! \mathbf{e}_{\nu}^\mathsf{T}\widehat{\boldsymbol{\beta}}_{p}(x)$ for $2\leq\nu\leq p$. We recommend using a local polynomial that is one order higher than the derivative to be estimated, $p=\nu+1$, or more generally to set $p-\nu$ odd. Of course, it is possible to achieve more bias reduction by increasing the local polynomial order $p$.
The \code{lpdensity()} function employs the triangular kernel by default. Other available options include the uniform kernel and the Epanechnikov kernel. Generally speaking, the uniform kernel delivers a smaller asymptotic variance but a larger asymptotic bias for the resulting point estimator. It is possible reduce its asymptotic bias by using a kernel that is more concentrated around the origin, such as the triangular kernel, at the cost of increasing the asymptotic variance. The choice of the kernel, however, does not affect the orders of the bias and the variance, and hence this is less of a concern compared to bandwidth selection, which we discuss in more detail below. In addition, Cattaneo-Jansson-Ma_2021_LocRegDistribution show that more efficiency gains can be achieved by first including a higher-order polynomial term and then partialling out this term using a minimum distance approach.
The rest of this section outlines the main properties, both statistical and numerical, of the estimator $\hat{F}^{(\nu)}_{p}(x)$, and discusses other related issues such as bandwidth selection, variance estimation, and valid (robust bias-corrected) inference. While the smooth CDF estimator $\hat{F}_{p}(x)=\hat{F}^{(0)}_{p}(x)$ is useful, the density and higher-order derivatives estimators are perhaps more relevant for empirical work. In particular, as mentioned before, the density estimator is intuitive and very easy to implement, while also being boundary adaptive. Thus, the estimator $\hat{f}_{p}(x)$ can be computed for all evaluation points on the (possibly restricted) support of $X$ in an automatic and straightforward way. This explains why the main functions in the package \pkg{lpdensity} refer to density estimation.
The estimator $\hat{F}_{p}^{(\nu)}(x)$ can be written in the familiar weighted least squares form: $\widehat{\boldsymbol{\beta}}_p(x) = (\mathbf{X}^\mathsf{T}\mathbf{W}\mathbf{X})^{-1}(\mathbf{X}^\mathsf{T}\mathbf{W}\mathbf{Y})$ with $\mathbf{X}$ the usual polynomial design matrix and $\mathbf{W}$ a diagonal matrix consisting of kernel weights. The only difference relative to standard local polynomial regression is that here the “dependent variable” is estimated: $\mathbf{Y}=[\hat{F}(X_1),\hat{F}(X_2),\cdots,\hat{F}(X_n)]^\mathsf{T}$, where $\hat{F}(x)$ is the possibly weighted empirical CDF. Unlike other local polynomial density estimators proposed in the literature \citep*[e.g.,][]{Cheng-Fan-Marron_1997_AoS,Zhang-Karunamuni_1998_JSPI}, the estimator $\hat{F}^{(\nu)}_{p}(x)$ does not require pre-binning or any other pre-processing of the data beyond constructing the empirical distribution function. As a result, this estimation approach removes the need of choosing the number, position, and length of the bins in a preliminary histogram estimate.
CJM obtained a general stochastic approximation to the bias and variance of the estimator $\hat{F}_{p}^{(\nu)}(x)$, $\nu=0,1,2,\cdots,p$, for all evaluation points $x\in\mathcal{X}$. Here we discuss the leading case of density estimation and derivatives thereof. For any choice of polynomial order $p$, and $1\leq\nu\leq p$, the variance and bias of $\hat{F}_p^{(\nu)}(x)$ are approximately
where $\mathsf{V}_{\nu,p}(x)$, $\mathsf{B}_{1,\nu,p}(x)$, and $\mathsf{B}_{2,\nu,p}(x)$ denote quantities that can be constructed directly using only the data, choice of (preliminary) bandwidth, evaluation point, polynomial order, derivative order, and kernel function. That is, all these quantities are in pre-asymptotic form, which has been shown to offer higher-order distributional refinements in the context of local polynomial regression \citep*{Calonico-Cattaneo-Farrell_2018_JASA,Calonico-Cattaneo-Farrell_2020_CEOptimal}. Furthermore, it can be shown that $\mathsf{V}_{\nu,p}(x)$, $\mathsf{B}_{1,\nu,p}(x)$, and $\mathsf{B}_{2,\nu,p}(x)$ converge (in probability) to well-defined non-random limits.
Since the above approximations are in pre-asymptotic form and are valid for all evaluation points, we can define a generic pointwise MSE-optimal bandwidth choice as \[h_{\mathtt{MSE}, p}(x) = \operatorname*{arg\,min}_{h>0} \mathsf{MSE}[\hat{F}_p^{(\nu)}(x)] = \operatorname*{arg\,min}_{h>0} \left\{\mathsf{Var}[\hat{F}_p^{(\nu)}(x)] + \mathsf{Bias}[\hat{F}_p^{(\nu)}(x)]^2\right\}. \] The optimal bandwidth also depends on $\nu$, but we suppress this in the notation to conserve notation. Under standard regularity conditions, $h_{\mathtt{MSE}, p}(x)$ is MSE optimal in rates for all evaluation points and choices of $p$ and $\nu$; and is MSE-optimal in constants if either (i) $p-\nu$ is odd or (ii) $x$ is a boundary point. See Appendix (ref) for more details.
We also define the IMSE-optimal bandwidth choice as follows: \[h_{\mathtt{IMSE}, p} = \operatorname*{arg\,min}_{h>0} \int\mathsf{MSE}[\hat{F}_p^{(\nu)}(x)]w(x)\mathrm{d}x, \] where $w(x)$ denotes a user-chosen weighting scheme. Dependence on $\nu$ is again suppressed to ease notation. In practice, the integral will be approximated using the grid points specified in \code{lpdensity()} or \code{lpbwdensity()}, allowing for both a uniform weighting ($w(x)=1$) as well as the empirical distribution weighting.
The MSE-optimal and IMSE-optimal bandwidth selectors, $h_{\mathtt{MSE}, p}(x)$ and $h_{\mathtt{IMSE}, p}$, are carefully developed so that they automatically adapt to boundary points, while also retaining their main theoretical features (e.g., rate optimality). In practice, these bandwidths can be computed after replacing the unknown quantities by estimates thereof. We discuss implementation details below.
Both $h_{\mathtt{MSE},p}(x)$ and $h_{\mathtt{IMSE},p}$, as well as their feasible counterparts, denoted by $\hat{h}_{\mathtt{MSE},p}(x)$ and $\hat{h}_{\mathtt{IMSE},p}$, can be used to construct MSE-optimal or IMSE-optimal point estimators for the PDF or its derivatives. The package \pkg{lpdensity} also allows for CDF estimation and computes an (I)MSE-optimal bandwidth, though we do not provide the details here to conserve space: some stochastic approximations change in non-trivial ways because the CDF estimator is $\sqrt{n}$-consistent. See CJM for more details.
As it is well known in the nonparametric literature, standard Wald-type inference is not valid when an (I)MSE-optimal bandwidth is used to construct the nonparametric point estimator. To be specific, the following distributional approximation holds for the standard Wald-type test statistic based on the local polynomial density estimator constructed using an (I)MSE-optimal bandwidth: \[\mathsf{T}_{\nu,p}(x) = \frac{\hat{F}_p^{(\nu)}(x)-F^{(\nu)}(x)} {\sqrt{\mathsf{Var}[\hat{F}_p^{(\nu)}(x)]}} \rightsquigarrow \mathcal{N}(\mathsf{bias},\ 1), \qquad 1\leq\nu\leq p, \] where $\rightsquigarrow$ indicates convergence in distribution, and $\mathcal{N}(\mu, \sigma^2)$ denotes the normal distribution with mean $\mu$ and variance $\sigma^2$. The bias term cannot be dropped in general: if the point estimator $\hat{F}_p^{(\nu)}(x)$ is constructed using an (I)MSE-optimal bandwidth, inference based on the usual “point estimator $\pm$ $z_{1-\alpha/2}$ $\times$ standard error” confidence interval is invalid due to the presence of an asymptotic bias ($z_\alpha$ denotes the $\alpha$ quantile of the standard normal distribution). The mechanical solution to this inferential problem is to undersmooth the point estimator $\hat{F}_p^{(\nu)}(x)$ using an ad hoc bandwidth $h$ smaller than $h_\mathtt{MSE}(x)$ or $h_\mathtt{IMSE}$. Of course, the function \code{lpdensity()} allows for this approach by simply running bandwidth selection, estimation, and inference in separate steps (see Section (ref) for an illustration).
\citet*{Calonico-Cattaneo-Farrell_2018_JASA,Calonico-Cattaneo-Farrell_2020_CEOptimal} showed that undersmoothing is suboptimal for inference under the same assumptions employed to construct an (I)MSE-optimal bandwidth. Instead, it is demonstrably better, in terms of higher-order distributional approximations and asymptotic coverage properties, to employ robust bias correction (RBC). The idea is to bias correct the point estimator and then adjust the variance accordingly. Heuristically, and abusing notation for simplicity, this leads to the Wald-type test statistic \[\mathsf{T}^{\mathtt{RBC}}_{\nu,p}(x) = \frac{\hat{F}^{(\nu), \mathtt{BC}}_p(x)-F^{(\nu)}(x)} {\sqrt{\mathsf{Var}[\hat{F}^{(\nu), \mathtt{BC}}_p(x)]}} \rightsquigarrow \mathcal{N}(0,1),\qquad \hat{F}^{(\nu), \mathtt{BC}}_p(x) = \hat{F}^{(\nu)}_p(x) - \mathsf{Bias}[\hat{F}_p^{(\nu)}(x)], \] which has a valid standard normal distribution even when an MSE, IMSE or a cross-validation-type bandwidth for $\hat{F}_p^{(\nu)}(x)$ is used. Confidence intervals with correct asymptotic coverage can be constructed by inverting the test statistic $\mathsf{T}^{\mathtt{RBC}}_{\nu,p}(x)$. In particular, it can be shown that a RBC confidence interval is equivalent to employing the test statistic $\mathsf{T}_{\nu,p+1}(x)=\mathsf{T}^{\mathtt{RBC}}_{\nu,p}(x)$ for a particular choice of parameters/implementation. Therefore, the function \code{lpdensity()} employs an RBC test statistic by default, assuming an (I)MSE-optimal or cross-validation-based bandwidth for the $p$-th order point estimator is used, denoted generically by $h_{p}$, and therefore forms the test statistic \[\mathsf{T}^{\mathtt{RBC}}_{\nu,p}(x)\equiv \mathsf{T}_{\nu,p+1}(x;h_p) = \frac{\hat{F}^{(\nu)}_{p+1}(x;h_p)-F^{(\nu)}(x)} {\sqrt{\mathsf{Var}[\hat{F}^{(\nu)}_{p+1}(x;h_p)]}}, \] and associated confidence intervals \[\mathsf{CI}_{\nu,p}^{\mathtt{RBC}}(x)\equiv \mathsf{CI}_{\nu,p+1}(x;h_p) = \left[\hat{F}^{(\nu)}_{p+1}(x;h_p)\ \ \pm\ \ z_{1-\alpha/2} \sqrt{\mathsf{Var}[\hat{F}^{(\nu)}_{p+1}(x;h_p)]}\right]. \] The notation makes the bandwidth explicit to distinguish the two polynomial degrees used in constructing the point estimator and the RBC confidence interval/test statistic: (i) a $p$-th order polynomial is used for point estimation (and bandwidth selection), and (ii) a $(p+1)$-th order polynomial is used for inference.
More generally, the package \pkg{lpdensity} implements confidence intervals of the form: \[\mathsf{CI}_{\nu,p}^{\mathtt{RBC},q}(x)\equiv \mathsf{CI}_{\nu,q}(x;h_p) = \left[\hat{F}^{(\nu)}_{q}(x;h_p)\ \ \pm\ \ z_{1-\alpha/2} \sqrt{\mathsf{Var}[\hat{F}^{(\nu)}_{q}(x;h_p)]}\right], \] with $q$ determining the inference approach. The above confidence interval is thus based on inverting the statistic $\mathsf{T}^{\mathtt{RBC},q}_{\nu,p}(x)\equiv \mathsf{T}_{\nu,q}(x;h_p)$, and by default we set $q=p+1$. CJM formally showed that the RBC confidence intervals have asymptotically correct coverage: \[ \lim_{n\to\infty}\mathbb{P}\left[ F^{(\nu)}(x) \in \mathsf{CI}_{\nu,p}^{\mathtt{RBC},q}(x) \right]= 1-\alpha,\qquad \forall x\in\mathcal{X}. \]
In addition to pointwise confidence intervals, the \pkg{lpdensity} package also offers uniform confidence bands for the CDF, PDF, or derivatives thereof. The uniform confidence band for $F^{(\nu)}(x)$ takes a similar form, \[\mathsf{CB}_{\nu,p}^{\mathtt{RBC},q}(\mathcal{G}) \equiv \mathsf{CB}_{\nu,q}(\mathcal{G};h_p) = \left\{\left[\hat{F}^{(\nu)}_{q}(x;h_p)\ \ \pm\ \ z_{\mathcal{G},1-\alpha/2} \sqrt{\mathsf{Var}[\hat{F}^{(\nu)}_{q}(x;h_p)]}\right],\quad x\in\mathcal{G}\right\}, \] with two noticeable differences. First, the confidence band no longer depends on the evaluation point, but rather on a collection of evaluation points, $\mathcal{G}$. Second, the critical value also changes, which is now denoted by $z_{\mathcal{G},1-\alpha/2}$. In practice, the new critical value can be obtained by first simulating a suitable Brownian bridge on the grid $\mathcal{G}$, and then computing the upper $\alpha$ quantile of the supremum of the simulated process. The option \code{CIuniform=TRUE} enables estimation and reporting of the uniform confidence band, which is turned off by default. CJM established a uniformly valid distributional approximation for the stochastic process $\{\mathsf{T}^{\mathtt{RBC},q}_{\nu,p}(x):x\in\mathcal{G}\}$, and proved that a nominal $1-\alpha$ level RBC confidence band is asymptotically valid: \[ \lim_{n\to\infty}\mathbb{P}\left[ F^{(\nu)}(x) \in \mathsf{CI}_{\nu,p}^{\mathtt{RBC},q}(\mathcal{G}),\ \forall x\in\mathcal{G} \right] = 1-\alpha. \] See \citet*{Cattaneo-Jansson-Ma_2021_LocRegDistribution} for technical details, regularity conditions, and additional discussions.
Robust bias correction methods lead to confidence intervals/bands that will not be centered at the density point estimates because of the recentering introduced by the bias correction. That is, different polynomial orders are used for constructing point estimates and confidence intervals/bands. Setting $q$ and $p$ to be equal delivers confidence intervals/bands that are centered at the point estimates, but requires undersmoothing for valid inference (i.e., an (I)MSE-optimal bandwdith cannot be used). Hence the bandwidth would need to be specified manually when $q=p$, and the point estimates will no longer be (I)MSE-optimal. Sometimes the point estimates may even lie outside of the confidence intervals/bands, which can happen if the underlying distribution exhibits high curvature at some evaluation point(s). One possible solution in this case is to increase the polynomial order $p$ or to employ a smaller bandwidth.
The package \pkg{lpdensity} implements several bandwidth selectors through \code{lpbwdensity()}, including MSE-optimal and IMSE-optimal plug-in rules, as well as rule-of-thumb bandwidth selectors based on a normal reference model. We only outline the main aspects of bandwidth selection here, but further details are given in Appendix (ref).
To introduce our bandwidth selectors, recall that the quantities $\mathsf{V}_{\nu,p}(x)$, $\mathsf{B}_{1,\nu,p}(x)$, and $\mathsf{B}_{2,\nu,p}(x)$ are given in pre-asymptotic form, and hence they can be computed from the data directly given a pilot/preliminary bandwidth. As a consequence, to construct the (I)MSE-optimal bandwidth, the only unknown quantities are $F^{(p+1)}(x)$ and $F^{(p+2)}(x)$, which can be consistently estimated using the local polynomial density derivative estimators implemented in \code{lpdensity()} with a pilot/preliminary bandwidth. To be more precise, the MSE-optimal bandwidth is estimated by \[\hat{h}_{\mathtt{MSE}, p}(x) = \operatorname*{arg\,min}_{h>0} \left\{\mathsf{Var}[\hat{F}_p^{(\nu)}(x)] + \widehat{\mathsf{Bias}}[\hat{F}_p^{(\nu)}(x)]^2\right\}, \] with $\widehat{\mathsf{Bias}}[\hat{F}_p^{(\nu)}(x)]$ constructed by replacing $F^{(p+1)}(x)$ and $F^{(p+2)}(x)$ with their estimated counterparts. Similarly, the IMSE-optimal bandwidth selector is given by \[\hat{h}_{\mathtt{IMSE}, p} = \operatorname*{arg\,min}_{h>0} \sum_{g_j\in \mathcal{G}}\left\{\mathsf{Var}[\hat{F}_p^{(\nu)}(g_j)] + \widehat{\mathsf{Bias}}[\hat{F}_p^{(\nu)}(g_j)]^2\right\}, \] where $\mathcal{G}$ is the collection of grid points specified in the function (by default, $\mathcal{G}$ takes on nineteen quantile-spaced values over the support of the data).
While the estimator $\hat{F}_p^{(\nu)}(x)$ is valid for all $\nu\geq0$, our discussion so far focused on the case $\nu\geq1$ because the resulting estimators of the density ($\nu=1$) and its derivatives ($\nu\geq2$) are the main focus of the package. Nevertheless, as a by-product, CJM developed analogous estimation, bandwidth selection and RBC inference results for the smooth CDF estimator $\hat{F}_p(x)=\hat{F}_p^{(0)}(x)$. These results are also implemented in the package \pkg{lpdensity} via the option \code{v=0}. For example, CDF estimation using a local constant approximation is obtained using \code{lpdensity(...,p=0,v=0)}, which employs the corresponding MSE-optimal bandwidth (\code{bwselect="mse-dpi"}) and a local linear approximation for inference (\code{q=p+1}) by default.
We showcase some of the main features of the \pkg{lpdensity} package. The data consists of 2,000 observations simulated from the normal distribution $\mathcal{N}(1,1)$ truncated from below at 0. We create a discontinuity in density at $x=0$ to illustrate the performance of our procedure at boundaries. Panel (a) of Figure (ref) plots a histogram estimate and the true density function.
The function \code{lpdensity()} provides both point estimates as well as RBC inference (confidence intervals and bands) employing the local polynomial density estimator, given a grid of points and a bandwidth choice. If the latter are not provided, then by default the function chooses nineteen quantile-spaced grid points over the support of the data and computes $\hat{h}_{\mathtt{MSE},p}(x)$ at each point.
The following command estimates the density function (\code{v=1}, the default) with fixed bandwidth \code{bw=0.5} at points $0, 0.5, \cdots, 4$, using a local quadratic approximation (\code{p=2}, the default) to the empirical distribution function. RBC confidence intervals over the grid are also computed, in this case using a local cubic approximation (\code{q=3}, the default).
The first part of the output provides basic information on the options specified in the function. For example, the default estimand is the density function, indicated by \code{Density function estimated (v=) 1}. The rest of the output gives estimation results, including (i) \code{Grid}: the grid points; (ii) \code{B.W.}: the bandwidths; (iii) \code{Eff.n}: the effective sample size for each grid point; (iv) \code{Point Est.}: the point estimates using polynomial order \code{p}, and the associated standard errors under \code{Std.} \code{Error}; (v) \code{Robust B.C.[95% C.I.]}: robust bias-corrected 95% confidence intervals. Point estimates, standard errors, and other information can be easily extracted for further statistical analysis. The output is stored in a standard matrix, and can be accessed with the following.
When the argument \code{grid} is suppressed, the evaluation points will be the $0.05, 0.1, \cdots,0.9, 0.95$ quantiles computed from the data. Conventional inference results (i.e., without robust bias correction) can be obtained by setting \code{q=p}. For example (output is suppressed):
It is also possible to suppress the argument \code{bw}, and the function will select the bandwidth automatically by minimizing (an estimated approximation to) the mean squared error, employing \code{lpbwdensity()}. Other bandwidth selection methods are available; we will illustrate data-driven bandwidth selection procedures in an upcoming subsection.
The method \code{summary()} takes six additional arguments. The first one, \code{alpha}, specifies the (one minus) nominal coverage of the confidence interval, with default being 0.05. Another argument is \code{sep}, which controls the horizontal separator. The default value is 5, and hence a dashed line is drawn after every five grid points. This feature can be suppressed by setting it to 0. Sometimes it may be desirable to report only a subset of the estimates, which can be done by using either the \code{grid} or the \code{gridIndex} option. The \code{grid} option allows reporting results for a selected set of grid points originally specified in the \code{lpdensity()} function, while \code{gridIndex} helps achieve the same goal by specifying the indices of the grid points. The last two options are related to confidence bands. By setting \code{CIuniform=TRUE}, a uniform confidence band, instead of pointwise confidence intervals, will be reported. Because the critical values have to be simulated in this case, the number of simulations used is controlled by the option \code{CIsimul} (its default value is 2000). The following example produces the 99% confidence band for four grid points $0$, $0.5$, $1$ and $2$, with dashed lines appearing after every three grid points. (Fixing the random seed allows reproducing the simulated critical values and the confidence intervals.)
Another important argument in \code{lpdensity()} is \code{scale}, which scales the point estimates and standard errors. This is particularly useful if only part of the data is used. For example, assume one would like to estimate the PDF using the two subsamples $\{X_i:X_i<1.5\}$ and $\{X_i:X_i>1.5\}$ separately. Simply splitting the data will not give consistent estimates, as it produces conditional (rather than marginal) density estimates:
The previous commands give point estimates $0.676$ and $1.222$, which are far from the true value $0.418$. To have consistent estimates, we need to scale the estimates by the proportion of the data used for estimation:
The function \code{plot()}, along with many other methods, is supported. This function takes the output from \code{lpdensity()} and produces plots of point estimates and robust bias-corrected confidence intervals/bands over the grid of evaluation points selected. Panel (b) of Figure (ref) shows how plots can be easily generated.
The confidence intervals/bands are not centered at the point estimates in general. As described in Section (ref), by default the point estimates are constructed using MSE-optimal bandwidths, which implies the smoothing bias is non-negligible and hence valid inference should be based on robust bias-corrected confidence intervals.
The function \code{plot()} allows for customization: Figure (ref) illustrates some of the features. (For panel (d), we again fix the random seed to reproduce the simulated critical values and the confidence band.)
The function \code{lpbwdensity()} implements four bandwidth selectors, (i) MSE-optimal plug-in bandwidth selector, denoted by \code{"mse-dpi"} (this is the default option), (ii) IMSE-optimal plug-in bandwidth selector, denoted by \code{"imse-dpi"}, (iii) rule-of-thumb bandwidth selector with a normal reference model, denoted by \code{"mse-rot"}, and (iv) integrated rule-of-thumb bandwidth selector, denoted by \code{"imse-rot"}. We illustrate some of the main features of \code{lpbwdensity()} with the same simulated data used previously.
By default, \code{lpbwdensity()} computes the MSE-optimal bandwidth for estimating the PDF with a local quadratic regression and triangular kernel, on nineteen quantile-spaced grid points: \code{lpbwdensity(..., p=2, v=1, bwselect="mse-dpi", kernel="triangular")}. The output resembles that of \code{lpdensity()}, and provides basic information for the data and options specified, as well as a matrix with three columns: (i) \code{Grid} for grid of evaluation points, (ii) \code{B.W.} for estimated bandwidths, and (iii) \code{Eff.n} for effective sample size at each grid point given the estimated bandwidth. The following is an example with a user-chosen grid of evaluation points.
The estimated bandwidths from this function can be used as input for \code{lpdensity()}, but constructing bandwidths in a separate step is redundant: bandwidth selection can be specified directly through the option \code{bwselect} in \code{lpdensity()}. For example, the following first computes the IMSE-optimal bandwidth and then estimates the density function:
It may be helpful to estimate bandwidths in a separate, first step so that they can be modified prior to estimation and inference (e.g., to implement ad hoc undersmoothing). To show this procedure, we reproduce panel (b) of Figure (ref) with the estimated IMSE-optimal bandwidth as well as ad hoc under-smoothing (where the IMSE-optimal bandwidth is divided by 2). See the following code and Figure (ref).
To prevent the estimated bandwidth from being too small, the default implementation in the \code{lpbwdensity()} function requires the local neighborhood around the evaluation point to contain at least $20+p+1$ (unique) observations. If the resulting neighborhood is not large enough, then the bandwidth is enlarged until the minimum number of observations is met. The default values can be changed through the options \code{nLocalMin}, controlling the minimum number of observations in each local neighborhood, and \code{nUniqueMin}, controlling the minimum number of unique observations in each local neighborhood. This minimal local sample size checking feature can be turned off by setting \code{regularize=FALSE}. Finally, the package \pkg{lpdensity} also includes checks and adjustments for repeated observations of the variable $X$ in the data. This feature can be turned off by setting massPoints=FALSE.
We illustrate the finite-sample performance of our \pkg{lpdensity} package in a simulation study, and compare it with other \proglang{R} packages implementing kernel-based density estimation procedures. The functions/packages we consider are: \code{bkde()} and \code{locpoly()} in the \pkg{KernSmooth} package Wand-Ripley_2015_KernSmooth, \code{kdde()} and \code{kde()} in the \pkg{ks} package Duong_2007_JSS,Duong_ks, \code{npudens()} and \code{npuniden.boundary()} in the \pkg{np} package Hayfield-Racine_2008_JSS,Racine-np, \code{kdrobust()} in the \pkg{nprobust} package \citep*{Calonico-Cattaneo-Farrell_2019_JSS,Calonico-nprobust}, \code{plugin.density()} in the \pkg{plugdensity} package Herrmann-Maechler_2011_plugdensity, as well as the built-in density estimator \code{stats::density()}.
Table (ref) provides a brief summary of their main features. First, three packages offer valid density estimates at (or near) boundaries, including \pkg{KernSmooth}, \pkg{np} and our \pkg{lpdensity}. However, only \pkg{KernSmooth} and \pkg{lpdensity} provide automatic boundary carpentry, while \pkg{np} requires specifying boundary kernels. Second, only two packages, \pkg{KernSmooth} and \pkg{lpdensity}, support higher-order bias reduction. Third, statistical inference is available in \pkg{np}, \pkg{nprobust}, and \pkg{lpdensity}. However, among these three packages, only \pkg{nprobust} and \pkg{lpdensity} account for the possibly leading smoothing bias when constructing test statistics/confidence intervals using (I)MSE-optimal bandwidths, and \pkg{nprobust} is not valid at or near boundary points. In addition, our \pkg{lpdensity} package is the only one that supports constructing uniform confidence bands. Fourth, only three packages, \pkg{KernSmooth}, \pkg{ks} and \pkg{lpdensity}, offer density derivative estimation. In summary, the \pkg{lpdensity} package provides valid density and derivatives estimation for both interior and boundary evaluation points, allows higher-order bias reduction through the use of higher-order local polynomial approximations, and offers several (I)MSE-optimal bandwidth selection methods. For statistical inference, \pkg{lpdensity} takes into account the possibly leading smoothing bias, and hence delivers (asymptotically) valid testings and confidence intervals/bands, both pointwise and uniformly over evaluation points.
We now describe our simulation design. The data consists of a random sample of size $n=1,000$, generated either from the normal distribution $\mathcal{N}(1,1)$ truncated below at $0$ (column “Truncated Normal”), or the exponential distribution with a scale parameter of 1 (column “Exponential”). We consider the estimation of the PDF at three evaluation points: $x=1.5$, $x=0.2$ and $x=0$, corresponding to interior, near boundary and boundary regions, respectively. We employ $2,000$ Monte Carlo repetitions. For the point estimate, we report its bias (column “Bias”), standard deviation (column “SD”) and root mean squared error (column “RMSE”). Whenever available, we also report the empirical coverage probability of a nominal 95% confidence interval (column “EC”) as well as its average length (column “IL”).
Simulation results are reported in Table (ref). At the interior evaluation point, all procedures perform similarly in terms of RMSE and empirical coverage. Point estimates obtained using \code{lpdensity()} have relatively small RMSEs, and the corresponding RBC confidence intervals exhibit satisfactory coverage properties. When the evaluation point is close to or exactly at the boundary, most packages/functions are no longer valid, and hence we only report simulation results for \code{locpoly()}, \code{npuniden.boundary()}, and \code{lpdensity()}. In such cases, \code{lpdensity()} delivers points estimates and confidence intervals with excellent finite-sample performance.
We gave an introduction to the general purpose software package \pkg{lpdensity}, which offers local polynomial regression based estimation and inference procedures for a cumulative distribution function, probability density function, and higher-order derivatives thereof. This package is available in both \proglang{R} and \proglang{Stata} statistical platforms, and further details can be found at \url{https://nppackages.github.io/lpdensity/}.
We thank Sebastian Calonico, David Drukker, Yingjie Feng, the editor, and two anonymous reviewers for thoughtful comments on our software implementation and article. We are also grateful to many users who provided valuable feedback. Cattaneo gratefully acknowledges financial support from the National Science Foundation (SES-1459931 and SES-1947805). Jansson gratefully acknowledges financial support from the National Science Foundation (SES-1459967 and SES-1947662) and the research support of CREATES.