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.
57,470 characters · 17 sections · 34 citation commands
Estimates of derivatives of (log) densities and related objects
We propose a new nonparametric estimator for (the logarithm) of a density function and its derivatives that attains the optimal rate of convergence both in the interior and at the boundary of the support. Our density estimator is available in closed form and is guaranteed to be positive unlike several alternatives, which is appealing in some applications and critical in others, such as in semiparametric maximum likelihood estimation klein1993efficient.\footnote{Klein and Spady square their density estimates to ensure positivity.} The new methodology differs from the previous literature in that it first estimates a function's derivatives, which, if desirable, can then be used to construct an estimate of the function itself. Our general estimation strategy can also be applied to obtain estimates of other quantities of economic interest, including the density in regression discontinuity design models, the (reciprocal) of the propensity score, the inverse bid function in auction models, and any other application in which the density appears inside a logarithm or a denominator.\footnote{The object of interest here is nonparametric in nature, i.e.\ it does not necessarily get averaged out as it does in e.g.\ lewbel2007simple.}
Specifically, we consider an i.i.d.\ sequence of random variables $\Set{\R x_1,\dots, \R x_n}$ with $\R x_i$ distributed according to some unknown distribution $F$ with density $f>0$ on its support $[0,\mathscr{U})$, where $\mathscr{U}$ can be infinite. The standard Rosenblatt--Parzen (RP) kernel density estimator is inconsistent at the boundary and is typically badly biased in finite samples at values of $x$ near the boundary. In contrast, our method employs a local polynomial approximation of $L\parens{x}=\log f\parens{x}$ to obtain asymptotically normal estimates of $L$ and its derivatives, away from, at, or near the boundary.
An advantage of using a polynomial approximation to the log--density instead of the density is that the estimated density can be guaranteed to be positive, which is not true for alternative boundary correction methods that use boundary kernels or a local polynomial approximation of $f$ Cheng1997auto,Karunamuni1998endpoints of $F$ Lejeune1992smooth,Cattaneo2019simple.\footnote{An exception is jones1996simple.} Unlike Loader1996likelihood and Hjort1996local, however, the computation of our estimator does not require solving a nonlinear system of equations that involves numerical integration. In fact, the estimator of derivatives of $L$ may be expressed as the solution to a linear system of (weighted) local averages. Thus, our method can be characterized as a local method of moments, similar in spirit to local likelihood density estimation Loader1996likelihood, Hjort1996local but computationally more similar to local polynomial regression Lejeune1992smooth, Cheng1997auto, Karunamuni1998endpoints,Cattaneo2019simple. We therefore retain the computational ease of a local polynomial regression while eliminating the possibility of negative density estimates. The estimator of $L,f$ itself then obtains in explicit form with estimates of the derivatives of $L$ as inputs.
Apart from its numerical advantages, our estimator for the density has the same first order asymptotic properties when applied with the same bandwidth as the local likelihood estimator. When applied with a larger bandwidth, however, our estimator achieves a smaller asymptotic mean squared error. We cannot generally compare the bias of our method with the biases of methods that use a polynomial approximation to $f$, but our estimator has the same asymptotic variance as traditional methods when they are applied using an optimal kernel and bandwidth sequence. Hence, our local polynomial approximation to $L$ can be expected to outperform alternative estimators for $f$ in finite samples when our bias is smaller, e.g.\ when the log--density is in fact polynomial.
In large enough samples, an asymptotically unbiased version of our estimator achieves a smaller variance and therefore a smaller mean square error than the optimized alternatives. Importantly, our estimator realizes this improved performance without sacrificing nonnegativity and continuity of the estimated density, as would be required to achieve the same asymptotic distribution using alternative methods.\footnote{One could replace negative density estimates produced by alternative methods by zero, but that is both clunky and would not help in cases in which the density must be positive.}
We also note that the log--density or its derivatives may be of direct interest to the researcher, in which case our method may be an attractive alternative to transforming estimates of $f$ and its derivatives to obtain the desired estimates. For instance, the generalized reflection method of Karunamuni2005boundary and Karunamuni2008some requires an estimate of $L'\parens{0}$ which they obtain using a finite difference approximation. Estimates of $L$ can moreover be used as an input into other objects. One case that has already been mentioned is semiparametric maximum likelihood estimation, of which klein1993efficient is a classical example in which the likelihood objective can be written as a function of the log--density. But there are other important examples. For instance, in regression discontinuity design, estimation of the density at the discontinuity point can be of interest Cattaneo2019simple. A second example would be the estimation of propensity scores which are of importance in the estimation of treatment effects. A final example is that of the estimation of auction models for which a version of our estimator can be used to obtain direct estimates of the inverse strategy function; see e.g.\ hickman2015replacing,pinkse2019estimation. These examples and more are discussed in (ref). Finally, we provide code in several languages, including Julia and R at \url{https://github.com/kschurter/logdensity}.
\makeatletter \makeatother
We now discuss our main estimator, postponing the discussion of applications and variants to (ref). Let $L\parens{y} = \log f(y)$ denote the log--density function and assume it to possess at least $S+1\geq 2$ derivatives at $x$, the point at which we wish to estimate $L$. Our estimator will be in the kernel family of estimators and we denote our bandwidth by $h$. To specifically allow for $x$ approximating the boundary, we introduce the notation $z = \min\parens{x/h, 1}$.
In the first step of our estimation procedure, we estimate derivatives of $L$, which are subsequently used to construct an estimate of $L$ itself. Our estimator of derivatives of $L$ is based on the fact that for any differentiable function $g:\mathbb{R} \to \mathbb{R}^S$ with support $\sparens{-z,1}$ and for which $g\parens{-z}=g\parens{1}=0\in\mathbb{R}^{S}$, we have
where $L^{\parens{s}}$ denotes the $s$--th derivative. The above follows from integration by parts under the assumption that $f$ is bounded on the domain of integration and a Taylor expansion of $L'$ around $x$. We define $\beta_s = L^{\parens{s}}\parens{x}$ and gather these coefficients into a vector $\beta = \sparens{\beta_1 \dots \beta_S}^{\mathpalette\raiseT{\intercal}}$. We then estimate integrals on the right and left sides of (ref) by their sample analogs and estimate $\beta$ by solving \< \frac1{nh^2} \sum_{i=1}^n g'\parens[\Big]{\frac{\R x_i-x}{h}} = -\sum_{s=1}^S \boldsymbol{\hat\beta}_s \frac1{nh}\sum_{i=1}^n g\parens[\Big]{\frac{\R x_i-x}{h}}\frac{(\R x_i - x)^{s-1}}{\parens{s-1}!}\,. \> Since (ref) is linear in $\boldsymbol{\hat\beta}$, the solution will generally be unique. Moreover, we can choose $g$ such that our estimator $\boldsymbol{\hat\beta}$ of $\beta$ is in closed form.
There are many functions $g$ that satisfy the desiderata outlined above. For the purpose of providing examples, we choose $g$ to be a vector whose $j$--th element is \[ g_j\parens{t} = \parens{t+z}^{j}\parens{1-t} \mathbb{1}\parens{-z \leq t \leq 1}. \] In (ref) we describe the sense in which this choice of $g$ is in fact optimal.
(ref) illustrates that the integration by parts in (ref) can be viewed as a device for obtaining a feasible set of local moment conditions from an infeasible set of moment conditions involving $L'$. Thus, $g$ fulfills a role similar to a kernel, but the restrictions we impose are different. Indeed, we require $g(-z)=g(1)=0$ so that $g(1)f(x+h)-g(-z)f(x-zh)$ is zero after integration by parts.
Now, once we have an estimator $\boldsymbol{\hat\beta}$ of the derivatives of $L$ at $x$, we can use it to construct estimators of $f\parens{x}$ and $L\parens{x}$. Indeed, substituting our approximation for the log--density in $\int_{x-zh}^{x+h} m_z\parens{x}f\parens{x}\:\mathrm{d} x$ and rearranging suggests an estimator for $\beta_0$ similar to that of Loader1996likelihood: \< \boldsymbol{\hat f}\parens{x} = \frac{\boldsymbol{\hat f}_m\parens{x}}{ \int_{-z}^1 m_z\parens{t} \exp\parens[\Big]{\sum_{s=1}^S \frac{\boldsymbol{\hat\beta}_s t^s h^s}{s!} } \:\mathrm{d} t }, \> where $\boldsymbol{\hat f}_m$ is a RP estimator using a nonnegative kernel $m_z$ with support $\sparens{-z,1}$. It should be apparent that $\boldsymbol{\hat f}\parens{x}$ cannot be negative and is zero only if there are no data on the interval $\sparens{x-zh,x+h}$.
If $z<1$ then $m_z$ can be thought of as a traditional boundary kernel\footnote{$k / \int_{\max\parens{x/h-1,0}} k$; see e.g.\ gasser1979kernel.} such that the bias in the numerator of (ref) is $O\parens{h}$. The role of the denominator is that it is an (asymptotically) biased estimator of the number one. Indeed, the denominator bias compensates for the numerator bias such that for $S=1$, the bias of $\boldsymbol{\hat f}\parens{x}$ is again $O\parens{h^2}$. We define $\boldsymbol{\hat L}\parens{x} =\log \boldsymbol{\hat f}\parens{x}$ and show that its asymptotic bias is $\beta_{S+1} h^{S+1}$ times a constant that is independent of both $n$ and fixed $x$.\footnote{For the case $x=zh$, the asymptotic bias does depend on $z$.}
Computing $\boldsymbol{\hat f}$ is relatively simple because $\boldsymbol{\hat\beta}$ is simply a local least squares statistic and $\boldsymbol{\hat f}$ is a ratio. Although $\boldsymbol{\hat f}$'s denominator contains an integral, for values of $S\leq 2$, which will be the most common scenario, the denominator in (ref) obtains in closed form if $m_z$ is a truncated Epanechnikov; for $S=1$ this is demonstrated in (ref) below.\footnote{For $S=2$ and $\boldsymbol{\hat\beta}_1<0$ it would however involve a normal distribution function and for $\boldsymbol{\hat\beta}_1>0$ a similar integral.} For other kernels and greater values of $S$, an asymptotically equivalent closed form expression can be obtained by expanding the denominator in (ref) in terms of exponential Bell polynomials bell1927partition.\footnote{For $S=1$ we get that that the exponential in the denominator in (ref) can be written as $1+\boldsymbol{\hat\beta}_1 th + \boldsymbol{\hat\beta}_1^2 t^2 h^2 + o_p\parens{h^2}$. For $S=2$ we get $1+\boldsymbol{\hat\beta}_1 th + \parens{2\boldsymbol{\hat\beta}_2 +\boldsymbol{\hat\beta}_1^2}t^2h^2 + \parens{3\boldsymbol{\hat\beta}_1\boldsymbol{\hat\beta}_2+\boldsymbol{\hat\beta}_1^3}t^3h^3 + o_p\parens{h^3}$.} Standard numerical integration methods can also be applied to this integral, in which case an advantage of our method is that this integral only needs to be computed once rather than at each iterate of a maximization routine, as in a local maximum likelihood approach.
The following examples compare the asymptotic behavior of our estimator and traditional approaches to kernel density estimation at and away from the boundary. (ref) obtains a closed form expression for the denominator in (ref), which is used in (ref) to obtain explicit expressions for the special cases $z=1$ and $z=0$ (away from the boundary and at the boundary, respectively).
In (ref) we took $x>h$ and hence $z=1$ to provide intuition. In the following example we consider what happens at the boundary, i.e.\ if $x=z=0$.
(ref) demonstrates that, unlike traditional boundary kernel estimators, the bias of our estimator is $O\parens{h^2}$ at the boundary, also. It may seem odd that the bias in (ref) is less than that in (ref) but note that the variance will be larger at the boundary and one would hence generally choose a greater bandwidth.
We first derive limit results for the vector $\boldsymbol{\hat\beta}$ of estimates of derivatives of $L$. Since our estimator $\boldsymbol{\hat\beta}$ is defined as the inverse of a matrix times a vector, its bias is the inverse of a matrix times a vector, also.
To simplify expressions for the asymptotic bias and variance of our estimator, we introduce the following objects which depend only on the choice of function $g$ (which in turn depends on the proximity to the boundary $z$), as well as a diagonal matrix that depends on $h$. Let \< \Omega_{js}=\Omega\parens{z} = \frac1{\parens{s-1}!}\int_{-z}^{1}-g_{j}\parens{t} t^{s-1}\:\mathrm{d} t = \sum_{t=0}^{s-1} \frac{\parens{-1}^{s-t+1}\parens{s-t} j!}{\parens{j+s+1-t}! t!} \parens{1+z}^{j+s+1-t}, \quad j,s=1,2,\dots, \> and let $\Omega,b,\Lambda$ be defined as \[ \Omega =
, \qquad b = \Omega^{-1}
, \qquad \Lambda =
. \] Let further $V \in \mathbb{R}^{S\times S}$ have $\parens{j,s}$ element equal to \[ \int_{-z}^1 g_j'\parens{t}g_{s}'\parens{t}\:\mathrm{d} t = \frac{2 js \parens{1+z}^{j+s+1}}{\parens{j+s+1}\parens{j+s}\parens{j+s-1}}. \] We are now in a position to state our first theorem.
The “in a neighborhood” condition comes from the fact that we specifically allow $x=zh$.
We now continue with the results for $\boldsymbol{\hat f}\parens{x}$.
Let $c_{msz} = \int_{-z}^1 m_z\parens{t} t^s \:\mathrm{d} t \mathrel/ s!$ and let $c_{mz}$ be a vector with elements $c_{m1z},\dots,c_{mSz}$. Let further \( \Omega_z\parens{t} = m_z\parens{t} - c_{mz}^{\mathpalette\raiseT{\intercal}} \Omega^{-1} g'\parens{t}, \) and define \( \mathscr{V} = f\parens{x}\int_{-z}^1 \Omega_z^2\parens{t} \:\mathrm{d} t \) and \( \mathscr{B} = f\parens{x}\beta_{S+1}\Xi_{f}\int_{-z}^{1}\omega_{z}\parens{t}t^{S+1}\:\mathrm{d} t \mathrel/\parens{S+1}!= f\parens{x}\beta_{S+1} \Xi_f \parens[\big]{ c_{m,S+1,z} - c_{mz}^{\mathpalette\raiseT{\intercal}} b} , \) for $\Xi_f$ a constant defined in the statement of (ref). Because $m_{z}$ integrates to one and $\int_{-z}^{1} g'\parens{t}\:\mathrm{d} t = 0$ and $\int_{-z}^{1}\Omega^{-1}g'(t)t^{s}\:\mathrm{d} t \mathrel/ s! =\Omega^{-1}\Omega_{\cdot s}$ is the $s$--th standard basis vector in $\mathbb{R}^{S}$, the function $\omega_{z}$ is a kernel of order $S+1$ or higher. To be clear, $\omega_{z}$ is not used to compute the density estimate; rather, it is a convenient object that arises in the asymptotic theory.
The asymptotic bias of our estimator is zero in some instances. For example, if $S=2$ and $z=1$ then the asymptotic bias is zero whenever $m$ is a symmetric kernel function; this is natural since this is effectively equivalent to choosing a higher order kernel $\omega_{z}$, albeit that unlike higher order kernel density estimates, our estimates cannot be negative.
The following two examples derive the $\Omega_z$ functions for the case in which both $m_z$ is a uniform and $x$ is at the boundary and the case in which $m_z$ is a truncated Epanechnikov and $x$ is anywhere.
{
}
As the above two examples demonstrate, deriving the asymptotic bias and variance for generic $z$ can be a messy but straightforward exercise.
In this section, we explore the optimal $(g,m_{z})$ in the local linear case $S=1$ and compare our optimized estimator with existing methods. We show that the above choice of $g$ and $m_{z}$ achieve the same asymptotic variance as an optimal RP estimator in the interior ($z=1$), while their respective biases cannot be compared in general. We then consider the optimal choice of $g$ and $m_{z}$ at the boundary ($z=0$), where we show that the truncated Epanechnikov $m_{z}$ paired with $g(t) = (t+z)(t-1)^{2}$ attains the same variance as an optimal boundary kernel Karunamuni1998endpoints, though the biases are again incomparable because our estimator's bias is a function of $f(x)L''(x)$ rather than $f''(x)$ as in the case of RP estimators. We are, however, able to compare the asymptotic performance of our estimation method with a local--likelihood based estimator. We show that our method with the cubic choice $g\parens{t}=\parens{t+z}\parens{1-t}^{2}$ attains the same asymptotic mean squared error (AMSE) in the interior and is more efficient at the boundary than the estimator in Loader1996likelihood with an Epanechnikov kernel.
Letting $\chi^{5} = \lim_{n\to\infty} h^{5}nf\parens{x}L''\parens{x}^{2}$, the AMSE of $\hat f\parens{x}$ only depends on $h$ and $\parens{g,m}$ through a multiplicative constant that can be written in terms of $\chi$ and the second--order kernel $\omega_{z}$: \< \sparens[\bigg]{\chi^{4} \parens[\bigg]{\int_{-z}^{1}\omega_{z}\parens{t}t^{2}/2 \:\mathrm{d} t}^{2}+ \chi^{-1} \int_{-z}^{1}\omega_{z}\parens{t}^{2}\:\mathrm{d} t}\sparens[\bigg]{f\parens{x}^{6/5}L”\parens{x}^{2/5}n^{-4/5}}\,. \> Unlike the typical approach to comparing kernels in kernel density estimation, in which one considers the optimal choice of $\chi$ as a function of $\omega_{z}$, we treat $\chi$ as fixed and seek to minimize the asymptotic MSE over $\omega_{z}$ instead of $(\omega_{z},\chi)$. We do so for two reasons. First, for values of $x$ near but not at the boundary, the function $\omega_{z}$ depends on $h$ through $z$, with the result that the first--order condition for optimality of the bandwidth is generally insufficient for the global minimum of the AMSE as a function of $h$. Second, many combinations of $g$ and $m_{z}$ yield a function $\omega_{z}$ that achieves zero asymptotic bias, which implies that there does not exist a finite optimal $\chi$.\footnote{One could assume an additional derivative of $f$, in which case the optimal bandwidth sequence would be proportional to $n^{-1/7}$ and one might also consider using a quadratic approximation ($S=2$).}
In light of the apparent similarity between (ref) and the corresponding expression for the AMSE of RP estimators, one might expect $\omega_{1}=3(1-t^{2})/4$ to be optimal using our method for the same reason that the Epanechnikov kernel is an optimal second--order kernel for use in RP estimation. Although we will eventually recommend $\omega_{1}=3\parens{1-t^{2}}/4$ for a particular value of $\chi$, our reasoning is different in two important ways. First, we do not require $\omega_{z}\geq 0$ as in Epanechnikov1969, because this restriction is not necessary to guarantee positive density estimates. Nor do we require $\omega_{z}\parens{1} = 0$ and $\omega_{1}\parens{-1}=0$ as in Muller1984 because these restrictions are not needed to ensure the density estimate is continuous in $x$. We place these restrictions on $m_{z}$, instead. Second, these constraints on $m_{z}$ and the maintained assumptions on $g$ are not binding if one minimizes the AMSE over $\parens{\omega_{z},\chi}$; for instance, $m_{1}\parens{t} = 3\parens{1-t}^{2}\parens{1+t}/4$ and $g\parens{t} = \parens{t+z}\parens{1-t}\parens{t^{2} + 2t - 1}$ yields $\omega_{1}\parens{t}=3\parens{3-5 t}^{2}/8$, which is the fourth--order kernel that minimizes the variance conditional on achieving zero bias and minimizes the AMSE as $\chi$ tends to infinity. Hence, when translated into the context of our estimator, the typical constraints on second--order kernels in RP estimation do not yield an interior solution to the optimal choice of $\parens{\omega_{z},\chi}$.
Thus, we seek a pair $\parens{g, m_{z}}$ with $m_{z}\geq 0$ that yields the optimal $\omega_{z}$ given a particular $\chi$, though we do allow $\chi$ to depend on $z$ so that the bandwidth sequence may be larger at the boundary than in the interior. The necessary conditions for the minimizer of the AMSE in (ref) imply that the optimal $\omega_{z}$ is quadratic, which generally implies a quadratic $m_{z}$ and a cubic $g$. The minimizer is not unique, however, because many pairs will produce the same $\omega_{z}$ and therefore the same asymptotic distribution. In fact, if the “first moment” of $m_{z}$ is zero, the AMSE in $\hat f$ does not depend on the choice of $g$ because $c_{m1}=0$ and $\omega_{z}=m_{z}$. A case such as this arises, for example, when $\chi = 15^{1/5}$, $z=1$, and one minimizes the AMSE by letting $m_{z}$ be the Epanechnikov kernel.
We focus on the choice of $\chi=15^{1/5}$ for two reasons. First, the constants in the expressions for the limiting bias and variance of $\hat f$ are the same as those found in the asymptotic bias and variance of the RP estimator for $f$ using the Epanechnikov kernel. This choice of $\chi$ and $m_{z}$ therefore provides a benchmark for comparison with an optimal RP estimator, because, for a fixed bandwidth $h$, the magnitude of our bias to the RP density estimator's bias depends on the ratio of $f(x)L''(x)$ to $f''(x)$ but our asymptotic variances are the same. And, second, the higher--order bias terms may be non-negligible when $\chi$ is too large, which could worsen the asymptotic approximation to the MSE in finite samples. Lacking a useful definition of “too large,” we default to a familiar choice.
Though we treat $\chi=15^{1/5}$ as fixed, this value is the optimal bandwidth scaling factor to use with the Epanechnikov kernel. We note, however, that this does not imply that the Epanechnikov $m_{z}$ and $\chi=15^{1/5}$ attain the minimum over all pairs $(m_{z},\chi)$. One can achieve a smaller AMSE at $z=1$ given a larger bandwidth by using a cubic $m_{z}$ and quartic $g$.\footnote{No symmetric nonnegative $m_{1}$ can improve on the AMSE of the Epanechnikov kernel because $g$ does not affect the limiting distribution. An asymmetric kernel---e.g.\ cubic $m_{1}$ with $m_{1}\parens{-1}=m_{1}\parens{1}=0$---and carefully selected $g$ is needed in order to improve on the AMSE at $z=1$.} Indeed, the above choice of $\omega_{1}(t) = 3\parens{3-5t^{2}}/8$ is an extreme example in which the asymptotic bias is zero. Its AMSE is smaller whenever $\chi > 3\times 15^{1/5}/2$, and its asymptotic variance is smaller as long as $\chi > 15^{6/5}\mathrel/8$, i.e.\ 1.875 times larger than the bandwidth used with the Epanechnikov kernel. Thus, even though our estimator's bias is generally not comparable to the bias of estimators that employ a local constant (RP) or local polynomial Lejeune1992smooth,Cattaneo2019simple approximation to $f$ or $F$, the asymptotically unbiased version of our estimator always has a smaller AMSE than these alternatives if the researcher is willing to use a large enough bandwidth.
Since the above criterion does not inform the optimal choice of $g$ for $z=1$ and $\chi=15^{1/5}$, one can choose $g\parens{t} = 1-t^{2}$ to minimize the AMSE in $\beta_{1}$.
At the boundary ($z=0$), it is perhaps reasonable to use a bandwidth that is twice as large as that used at $z=1$, i.e. $\chi = 2 \times 15^{1/5}$. In this case, the optimal AMSE is attained with the truncated Epanechnikov and $g(t) = t (1-t)^{2}$. Interestingly, this choice implies that $\omega_{0}\parens{t} = 6\parens{1-2t}\parens{1-t}$, which is the optimal boundary kernel derived by Karunamuni1998endpoints. Indeed, the constants in our bias and variance expressions are the same as the kernel--related constants in the limiting distribution of the RP estimator for $f\parens{0}$ given by $\sum_{i=1}^{n}\omega_{0}\parens{x_{i}/h}/\parens{nh}$, indicating that the asymptotic variances of the two estimators are the same for a fixed bandwidth. In contrast to this estimator, however, our proposed estimator is always positive.
As in the interior case ($z=1$), this choice of $m_{z}$ and $g$ is not optimal over all triples $\parens{g,m_{z},\chi}$. One can obtain zero asymptotic bias and a smaller asymptotic variance using the truncated Epanechnikov $m_{z}$, $g(t) = \parens{t+z}\parens{t-1}\parens{t-5/7}$, and any finite $\chi>15^{6/5}/4$.\footnote{These inputs yield the third--order kernel boundary $\omega_{0}\parens{t}= 9 - 36 t + 30 t^{2}$ that minimizes $\int_{0}^{1} \omega_{0}\parens{t}^{2}\:\mathrm{d} t$} But we caution that this relatively large bandwidth---at least 3.75 times larger than the optimal bandwidth in the interior---can limit the usefulness of our asymptotic approximation to the bias in finite samples.
Finally, we note that the asymptotically unbiased $\omega_{0}\parens{t}$ and $\omega_{1}\parens{t}$ are the same as those derived for the $S=2$ case in (ref). Hence, the asymptotically unbiased local linear estimator has the same limiting distribution as an undersmoothed local quadratic estimator, i.e.\ a local quadratic estimator using a bandwidth sequence of order $n^{-1/5}$ instead of $n^{-1/7}$, though they are not numerically equivalent.
The AMSE of the local--likelihood estimator for $f\parens{x}$ in Loader1996likelihood can be written in a form similar to (ref). The relative AMSE of our proposed estimator and the local likelihood estimator using an optimal kernel is then given by the ratio of the multiplicative constants that scale $f\parens{x}^{6/5}L''\parens{x}^{2/5}n^{4/5}$. At $z=1$, if one uses an optimal bandwidth sequence and the Epanechnikov kernel with the local--likelihood based estimator and $\chi=15^{1/5}$ with our optimal estimator, the relative AMSE of our estimators is one. In fact, the limiting distributions are identical. At $z=0$, the optimal kernel to use with the local likelihood estimator is triangular, i.e. $k\parens{t}=\parens{1-|t|}\mathbb1\cparens{|t|\geq 1}$, which yields the same asymptotic bias and variance as our estimator using $\chi = 2\times15^{1/5}$, the truncated Epanechnikov, and $g\parens{t}=\parens{t+z}\parens{1-t}^{2}$.
We should expect our estimator with $g\parens{t}=\parens{t+z}\parens{1-t}^{2}$ and the local--likelihood estimator to perform similarly in finite samples. Thus, our density estimator's computational ease is its more salient advantage over the local--likelihood based approach, given these inputs. Of course, one could make an alternative choice of $g$ and increase the bandwidth to widen the gap in AMSE at the possible expense of a larger finite sample bias. Greater reductions in the AMSE require larger bandwidths and risk worse finite sample performance.
For $S>1$, the optimal $\parens{g,m_{z}}$ can again be reformulated as the optimal choice of a higher order kernel $\omega_{z}$. Unlike in RP estimation, however, the higher order kernel does not necessarily entail the possibility of negative density estimates, since one can achieve a higher order $\omega_{z}$ kernel using a nonnegative $m_{z}$ and a suitable $g$. Moreover, as in the linear case, the restrictions $\omega_{z}\parens{1}=0$ and $\omega_{1}\parens{-1} = 0$ are not necessary in order for the estimated density to be continuous. Without these restrictions we do not obtain an interior solution to the optimal combination of bandwidth and $\omega_{z}$. We would therefore fix the bandwidth when we optimize the AMSE over $\parens{g,m_{z}}$ in the higher order case, as well.
It is well--known hirano2003efficient that under an unconfoundedness assumption the average treatment effect can be expressed as \[ \mathbb{E}\parens[\Big]{\frac{ \R y \R t}{p\parens{\R x}} - \frac{ \R y \parens{1-\R t}}{1-p\parens{\R x}}}, \] where $p$ is the propensity score, $\R y$ the outcome variable, $\R x$ a vector of regressors, and $\R t$ a binary treatment variable. Let $p_1= \mathbb{E} \R t$ be the unconditional treatment probability. Let further $f_1$ denote the regressor density function conditional on treatment and $f_0$ the density conditional on nontreatment. Then $f\parens{x} = f_1\parens{x} p_1 + f_0\parens{x} \parens{1-p_1}$, which produces \[ \frac1{p\parens{x}} = 1 + \frac{f_0\parens{x}}{f_1\parens{x}} \frac{1-p_1}{p_1}, \qquad \frac1{1-p\parens{x}} = 1+ \frac{f_1\parens{x}}{f_0\parens{x}} \frac{p_1}{1-p_1}, \] such that the reciprocals of the propensity scores only depend on the ratios of the densities and the unconditional choice probabilities. In practice, $p\parens{x}$ is often estimated using a logistic functional form and possibly a series expansion in $x$,\footnote{This casual observation is supported by the fact that the built--in propensity score matching estimator in Stata defaults to the logit model.} which implies the logarithm of the odds ratio is a polynomial in $x$. Our log--polynomial approximation to $f_{1}$ and $f_{0}$ similarly imply the odds ratio is log--polynomial. The data $\R x$ will not typically be scalar--valued, so one could for instance use a linear index of regressors instead of the regressors themselves; see (ref) for an example of how one might estimate $p\parens{x}= p^*\parens{x^{\mathpalette\raiseT{\intercal}}\theta_0}$. In any case, our approach is a natural local extension to the logit series estimator for $p\parens{x}$.
In first--price, sealed--bid procurement models with independent private values, it is well--known Guerre2000 that the inverse bid function is of the form $b - \bar F\parens{b} / f\parens{b}$, where $\bar F,f$ are the survivor and density functions of the minimum rival bid. Since the support of the bid distribution is assumed to have a lower bound in this literature (costs cannot be less than zero and hence neither are bids), boundary issues are a serious concern.
So let $Ψ\parens{y} = \bar F\parens{y} / f\parens{y}$ be the object of estimation. One way of estimating $Ψ$ is to estimate $\bar F,f$ separately where $f$ is estimated using the machinery in the main part of this paper and $\bar F$ is estimated by the empirical survivor function. This estimator has all the features of the estimator discussed earlier in the paper. In particular, if the underlying cost distribution is approximately an exponential then so is the bid distribution and our estimator could be expected to work especially well.
The above approach is not specific to auction models. Indeed, consider the hazard function $H\parens{x} = f\parens{x} / \bar F\parens{x}$. $H$ can also be estimated using the machinery developed in our paper.
There are many examples of semiparametric maximum likelihood estimators. Here, we only consider a classical ones, namely the klein1993efficient estimator of the coefficients in a semiparametric binary response model, which maximizes \[ \sum_{i=1}^n \sparens[\big]{ \R y_i \log \boldsymbol{\hat p}\parens{\R x_i^{\mathpalette\raiseT{\intercal}} \theta} + \parens{1-\R y_i} \log \cparens{ 1- \boldsymbol{\hat p}\parens{\R x_i^{\mathpalette\raiseT{\intercal}} \theta}}}, \] where $\boldsymbol{\hat p}$ is an estimator of the choice probability. Klein and Spady apply techniques to ensure that the estimates $\boldsymbol{\hat p}$ are positive and less than one, including trimming and adding a sample--size--dependent constant. Our method could be helpful since \( p\parens{t} = \Pr\condr{\R y_1=1 \nonscript\:\delimsize\vert \allowbreak \nonscript\: \mathopen{} \R x_1^{\mathpalette\raiseT{\intercal}}\theta=t} = f_1\parens{t} p_1 / f\parens{t}, \) where $f_j$ is the density of the linear index for observations with $\R y_i=j$ and $p_1$ is the unconditional choice probability. Since $f\parens{t} = p_1 f_1\parens{t} + \parens{1-p_0} f_0\parens{t}$, the infeasible contribution to the loglikelihood could be written as \[ \R y_i \log \frac{f_1\parens{\R x_i^{\mathpalette\raiseT{\intercal}} \theta}}{f_0\parens{\R x_i^{\mathpalette\raiseT{\intercal}}\theta}} - \log \parens[\bigg]{ p_1\frac{f_1\parens{\R x_i^{\mathpalette\raiseT{\intercal}}\theta}}{f_0\parens{\R x_i^{\mathpalette\raiseT{\intercal}}\theta}} + \parens{1-p_1}} + \text{constant}, \] such that it is only the ratio of $f_1/f_0$ that matters. Obtaining conditions under which our estimator obtains the semiparametric efficiency bound, like the Klein and Spady estimator does, are well beyond the scope of this paper.
One context in which the behavior of estimates near or at the boundary is of special importance is that of regression discontinuity design. For instance, Cattaneo2019simple provide a test of continuity of the density function at the boundary using a boundary density estimator that is similar to the estimator in Lejeune1992smooth in that it is based on a quadratic expansion of the distribution function. Compared to that approach, our method requires that the density be nonzero at the boundary, which is a requirement for the regression discontinuity framework in any case. The bottom line is that our method will work better if the log density is approximately a low order polynomial near the boundary and theirs if the density itself is approximately a low order polynomial. This is borne out by our simulation results.
The boundary correction method of Karunamuni2008some requires a well--behaved estimate of $L'\parens{0}$, which is exactly what our method provides.
The following simulation exercise compares the performance of our estimator for the density and its derivatives near the boundary with alternatives that also employ local polynomial approximations Lejeune1992smooth, Cattaneo2019simple,Loader1996likelihood and the generalized reflection method, which also estimates the derivative of the log--density near the boundary to remove the boundary effects of the RP estimator Karunamuni2005boundary,Karunamuni2008some. For the local polynomial estimators, we use a local linear approximation to the density or log--density, depending on the method.\footnote{This corresponds to a local--quadratic polynomial approximation to the distribution function in Lejeune1992smooth and Cattaneo2019simple.}
We simulate 2000 i.i.d.\ samples of size $n=500$ and estimate $f$ at points within two bandwidths of the boundary in order to compare the estimators away from, near, and at the boundary. The random variables are drawn from each of four parametric distributions whose densities exhibit varying behaviors near their left boundary $x=0$. The first is a beta distribution rescaled to take support on $[0,5]$ (so that right boundary is sufficiently far away from zero) with density $f_{1}\parens{x} = \theta \parens{1-x/5}^{\theta-1}/5$, which is in fact polynomial in $x$ for integer values of $\theta$, which might favor CJM, although that is not reflected in our simulations if $\theta>3$. The second design is a normal distribution with a mean of $\theta/2$ and variance of one, truncated at zero. This density is log--quadratic, which should favor our method and Loader's. The third and fourth designs are $f_{3}\parens{x} = \parens{e^{-x} + \theta x e^{-x}}/\parens{1+\theta}$ and $f_{4}\parens{x} = \parens{e^{-x} + \theta x^{2}e^{-x}}/\parens{1+2\theta}$.
For each simulation design, we estimate $f$ and its derivative using our approach with $g_{j}\parens{t}=\parens{t+z}^{j}\parens{t-1}$, $g_{j}\parens{t}=\parens{t+z}^{j}\parens{1-t}^{2}$, and $g_{j}\parens{t}=\parens{t+z}^{j}\parens{t-1}\parens{t-5/7}$ (PS$_{1}$, PS$_{2}$, PS$_{3}$), Loader's local likelihood estimator (Loader), a local polynomial regression of the empirical CDF (LS--CJM), and a generalized reflection estimator (KZ). Wherever a kernel is required, we use the Epanechnikov kernel $k(u) = 3 \parens{1- u^{2}} /4$ or a truncated version thereof. This choice is not optimal at the boundary for Loader's estimator, but the efficiency loss is quite small.\footnote{The relative efficiency of Loader's estimator using the triangular and Epanechnikov kernels at the boundary is about 1.008, meaning the Epanechnikov kernel requires a sample size 1.008 times larger to achieve the same MSE as the triangular kernel. We therefore expect the Epanechnikov kernel with $n=504$ to have the same MSE as the triangular kernel with $n=500$.}
Where possible, we use the asymptotically optimal bandwidth sequences for $z=1$ and $z=0$. For intermediate values of $z$, we linearly interpolate the bandwidth.\footnote{Specifically, we use a bandwidth $h = h_{0}\parens[\big]{1-\min\cparens[\big]{\frac{x}{h_{1}},1}} + h_{1}\min\cparens[\big]{\frac{x}{h_{1}},1}$, where $h_{0}$ and $h_{1}$ are the asymptotically optimal bandwidths at $z=0$ and $z=1$ and $x$ is the point of evaluation. This bandwidth selection rule implies that the same window is used to estimate the density and its derivatives at all points within $h_{1}$ of the boundary, i.e.\ $x+h = 2h_{1}$ for all $x<h_{1}$.} For the asymptotically unbiased version of our density estimator, PS$_{3}$ with $S=1$, there is no finite optimal bandwidth unless we assume more derivatives of $f$. Instead, we choose the bandwidth at $z=0$ so that the asymptotic variance is the same as the variance of PS$_{2}$ at the boundary. For the generalized reflection method, which requires separate bandwidths and finite--difference approximation to estimate $L'$ in a first step, we do not develop a theory of the asymptotically optimal inputs. Instead, we select a finite--differencing scheme and choose a combination of auxiliary bandwidths so that the pilot estimate of $L'$ at the boundary and the density estimate at $z=1$ have the same asymptotic variances as our method using $g\parens{t}=\parens{t+z}\parens{1-t}$. Specifically, we use a main bandwidth equal to the asymptotically optimal bandwidth for our method at $z=1$, and we estimate $L'\parens{0}$ using $\cparens{\log f_{nh_{L'}}\parens{2h_{L'}}-\log f_{nh_{L'}}\parens{0}}\mathrel/\parens{2h_{L'}}$, where $f_{nh_{L'}}\parens{2 h_{L'}}$ and $f_{nh_{L'}}\parens{0}$ are a kernel and boundary--kernel estimator for the density whose variance is proportional to $1/nh_{L'}$.
Comparing the square root of the mean squared error (RMSE) of the density estimates at zero in (ref) and the RMSE for the derivative in (ref), there is no clear ranking of the estimators. (ref) depicts the bias and RMSE of the local linear estimators for $\theta=4$. The linear approximation of $f$ (LS--CJM) performs well when $f$ is a beta distribution, but has difficulty estimating the truncated normal density and the estimate is often negative. KZ is generally neither the best nor the worst of the estimators we consider here, but we acknowledge that we have not optimized the inputs into the KZ estimator as thoroughly as the other estimator's inputs. We also note that KZ provides a familiar benchmark away from the boundary because it is simply the RP density estimator using an Epanechnikov kernel for $z=1$.
As expected, PS$_{1}$, PS$_{2}$, PS$_{3}$, and Loader perform similarly away from the boundary. The differences between these estimators in the top two figures are due to the relatively large bandwidth, but the curves are nearly indistinguishable in the bottom two figures. At the boundary, Loader appears to consistently achieve a smaller RMSE than PS$_{1}$ as our asymptotic theory predicts. The difference between PS$_{2}$ and Loader is almost indiscernible in this sample size, though we expect PS$_{2}$ to have a slightly smaller asymptotic bias and variance.
While our estimator can closely mimic the performance the local likelihood estimator, our method can also achieve significantly smaller AMSE if we use a larger bandwidth and an alternative $g$. As an extreme example, the simulations results show the bias of PS$_{3}$ is relatively small in all four designs; in fact, it converges to zero faster than $h_{0}^{2}$. But its MSE appears to be roughly the same at the boundary as the other estimators' and is generally larger for $z>0$ in three of the four designs, indicating that the finite--sample costs of the asymptotic benefits do not justify this ambitious choice of $g$. The notable exception is in the case of the truncated normal, where the higher--order bias terms are in fact zero due to the fact that $f_{2}$ is log--quadratic. As a result, PS$_{3}$ has a markedly smaller RMSE at the boundary. One could also eliminate the asymptotic bias and significantly reduce the MSE for other values of $z$. For example, in the interior one could use $m_{z}\parens{t}=\frac{3}{4}\parens{1-t}^{2}\parens{1+t}$, $g\parens{t}=\parens{t+z}{1-t}\parens{t^{2}+2t - 1}$, and a bandwidth that is 1.875 times larger than that used for Loader and our other estimators, as suggested in (ref). In fact, if the density possesses more derivatives than the researcher was willing to assume, the asymptotic gains will typically be realized more quickly in the interior than at the boundary because the third--order bias term is zero, as well.
We interpret these simulation results as a proof of concept that our approach can improve on the asymptotically optimal local polynomial and RP estimators without sacrificing continuity or nonnegativity of the estimated density. Moreover, we demonstrate these gains are possible in empirically relevant sample sizes even when the researcher ambitiously attempts to eliminate the asymptotic bias at the boundary. In practice, however, researchers might prefer less extreme versions of our estimator that do not require such large bandwidths to achieve a lower AMSE than commonly used alternatives. Indeed, the asymptotically unbiased version of our estimator does not minimize the AMSE for any bandwidth sequence on the order of $n^{-1/5}$, and it would only be advisable if the researcher specifically requires an unbiased estimate.
In (ref) we plot the bias and RMSE for the derivative of $f$. In three of the four simulation designs, PS$_{1}$ has the smallest RMSE among the log--linear approximation methods in the interior region, which we expected because $g\parens{t}=\parens{t+z}\parens{t-1}$ minimizes the AMSE in $L'$ using the given bandwidth sequence. At the boundary, however, PS$_{1}$ has a significantly larger AMSE than the alternative choices of $g$, and this is borne out in the simulations to some extent even though the estimator for the derivative converges at the relatively slow rate of $n^{-1/5}$.
We develop an asymptotically normal nonparametric estimator based on a log--polynomial approximation to the unknown density. By approximating the log--density with a polynomial, we can guarantee our estimated density is nonnegative; and by using a polynomial approximation instead of a local constant approximation, we achieve the optimal rates of convergence at the boundary of the support as well as in the interior.
Because our approach allows for a relatively larger degree of customization---the researcher must specify a bandwidth, a kernel, and a vector--valued function $g$ that is zero at the extremes of its support---we explore the optimal set of inputs. Unlike the standard analysis of optimal kernel and bandwidth inputs, our estimator is nonnegative and continuous in the point of evaluation under a relaxed set of constraints. Because these constraints were needed in order to derive the optimal kernel for use with alternative methods, there is no interior solution to the optimal choice of inputs using our approach. If one fixes the bandwidth sequence, however, the choice of kernel and $g$ can be optimized in a straightforward manner, and our method can achieve the same asymptotic variance as the optimized alternatives if one uses the same bandwidth. Moreover, if the researcher uses a slightly larger bandwidth, marginal reductions in the asymptotic mean squared error are possible.
In fact, there is no positive lower bound to the asymptotic mean squared error using our approach if the researcher is willing to adopt a large enough bandwidth sequence. Although such large bandwidths stretch the plausibility of the asymptotic approximation to the mean squared error in finite samples, we provide simulation evidence that demonstrates significant reductions in the mean squared error can be realized in empirically relevant sample sizes.
More generally, our approach is based on a sample analog to partial integration and can be applied to other settings, as well, such as the estimation of hazard functions or propensity scores.