Shape-Constrained Density Estimation via Optimal Transport

Shape-Constrained Density Estimation Via Optimal Transport

Introduction

Nonparametric density estimation has the advantage over its parametric counterparts of not requiring the underlying population density to belong to a specific family. In the case of distribution functions, Kiefer and Wolfowitz (1956) showed that the empirical distribution function is a maximum likelihood estimator; however, attempting to use this distribution function to directly define a nonparametric density estimate results in a series of point masses located at each of the datapoints. Grenander (1956) provided the first example of a shape constrained density estimator as a way to extricate the maximum likelihood estimator from this “Dirac catastrophe.” Specifically, he showed that maximizing the likelihood function subject to a monotonicity constraint on the estimator results in density estimates without point masses.

A great deal of progress was made in subsequent decades by adding penalty terms to the maximum likelihood objective function to restore the parsimony of the density estimator; for example, see (Silverman, 1986). Parzen (1962) also showed that kernel density estimators resulted in consistent density estimators and derived the rates of convergence. Unlike Grenander's (1956) approach, the performance of maximum penalized likelihood estimators and kernel density estimators is highly dependent on the specification of penalty terms and bandwidths, respectively, which can be difficult to choose.

Partly for this reason, recently there has been a renewed interest in ensuring parsimony of the maximum likelihood density estimator through conditioning on the information provided by the shape of the underlying density. In particular, significant progress has been made on the maximum likelihood density estimator subject to the constraint that the logarithm of the density is a concave function, which defines a $\log-$concave density (Dümbgen and Rufibach, 2009; Cule, Samworth, and Stewart, 2010; Kim and Samworth, 2016).

One early pioneer on the advantages of $\log-$concavity for both statistical testing as well as estimation was Karlin (1968). Suppose the distribution function $F:\mathbb{R}\rightarrow[0,1]$ has a density function denoted by $f:\mathbb{R}\rightarrow\mathbb{R}_{+}.$ Some examples of $\log-$concavity's many implications include that the density $f(x-\theta)$ has a monotonic likelihood ratio if and only if $f(\cdot)$ is $\log-$concave, products and convolutions between $\log-$concave densities are $\log-$concave, and that the hazard function of the $\log-$concave density $f(x),$ defined by $f(x)/(1-F(x)),$ is increasing. Bagnoli and Bergstrom (2005) also provide a survey of economic models in which $\log-$concavity of a density is a sufficient condition for the existence or uniqueness of an equilibrium. Chen and Samworth (2013) as well as Dümbgen, Samworth, and Schuhmacher (2011) provide tests for a population density satisfying $\log-$concavity, and Carroll, Delaigle, and Hall (2011) provide a test for a population density satisfying a more general set of shape constraints.

A wide variety of the random variables in the economics literature, such as annual income or changes in stock prices, are thought to exhibit sub-exponential tail behavior, so a $\log-$concavity constraint would not result in a consistent estimator in these cases. Koenker and Mizera (2010) generalized the $\log-$concave maximum likelihood estimator by maximizing Rényi entropy of order $\rho\in\mathbb{R}$ subject to the $\rho-$concavity constraint,

for all $\alpha\in[0,1].$ This estimator converges to the maximum likelihood estimator subject to a $\log-$concavity constraint in the limit as $\rho\rightarrow0.$\footnote{Maximum likelihood is equivalent to maximizing Shannon entropy, and Rényi entropy of order $\rho$ converges to Shannon entropy as $\rho\rightarrow0$ .} Decreasing $\rho$ corresponds to a relaxation of this shape constraint, so if $f(x)$ satisfies the constraint for some $\rho,$ then it also satisfies the constraint for all $\rho'<\rho.$ Also, this constraint is equivalent to concavity when $\rho$ is equal to one, and the cases of $\log-$concavity and quasi-concavity can be derived in the limit as $\rho\rightarrow0$ and $\rho\rightarrow-\infty$ respectively. Koenker and Mizera (2010) place particular emphasis on the case in which $\rho=-1/2,$ partly because most standard densities are $-1/2-$concave. For example, all Student $t_{v}$ densities with $v\geq1$ satisfy this constraint.

$\rho-$concavity constraints provide a considerable relaxation over $\log-$concavity constraints, while restricting the set of feasible densities sufficiently to ensure parsimony of the density estimator. These constraints are also sufficient conditions for many results in economics, including the uniqueness or existence of equilibria in a variety of models; see for examples, (Ewerhart, 2013; Bagnoli and Bergstrom, 2005). However, in many cases the necessary and sufficient conditions for these results are considerably weaker, so inference and estimation based on these stronger conditions can provide misleading results. For example, inferring whether a population density satisfies these weaker conditions based on tests for their more restrictive counterparts is generally not possible. Things are less straightforward in the case of estimation because shape constraints are generally the source of the density estimator's parsimony. However, since using a shape constraint that is not satisfied by the population density would not result in a consistent estimator, it is prudent to err toward the weakest constraint that theory predicts a population density would satisfy, or when using the estimate in a structural model, the weakest constraint that a model requires a population density to satisfy.

For a concrete example in the economics literature, given a density of private valuations of risk neutral agents, Myerson (1981) defined the virtual valuations function as $x-(1-F(x))/f(x),$ and showed that the first price auction is revenue maximizing if this function is strictly increasing. Sufficient conditions for Myerson's (1981) regularity condition, ordered from strongest to weakest, are $\log-$concavity, a monotonic hazard rate, and $\rho-$concavity for $\rho>-1/2$ (Ewerhart, 2013). While the first two sufficient conditions are commonly cited in mechanism design, they both imply exponential tail behavior. Since it seems plausible that willingness to pay is influenced by ability to pay, allowing valuations to have sub-exponential tail behavior may be a reasonable modeling choice, given that wealth and income are typically modeled with sub-exponential tails. Luckily, Myerson's (1981) regularity condition does not exclude these densities. For example, while the $\log-$normal density has sub-exponential tails, it also satisfies this condition when $\sigma^{2}<2,$ which holds in the structural model provided by Laffont, Ossard, and Vuong (1995).

This paper provides a framework for estimating and performing inference with shape constrained densities using regularized optimal transport (Cuturi, 2013). This objective function has the advantage of having an unconstrained global optimum that is a consistent density estimator, which ameliorates the requirement that the shape constraint is the only source of parsimony. At first we motivate the method using a $\rho-$concavity constraint, but one of the advantages of the method is that the estimator is consistent when this constraint is replaced by a wide variety of alternative shape constraints. We also provide a consistent test for whether or not a population density satisfies a shape constraint based on comparing the objective function at the unconstrained optimum to the constrained optimum.

After introducing density estimation with this more general class of shape constraints, we use the proposed estimator to explore whether or not the California Department of Transportation's decision to use a reverse first price auction to award construction contracts is cost minimizing. To do this, we use the method provided by Guerre, Perrigne, and Vuong (2000) to calculate the firms' marginal costs using data on their bids. We find that a kernel density estimator of these costs does not satisfy Myerson's (1981) regularity condition everywhere; however, the proposed density estimate, subject to the constraint that Myerson's (1981) regularity condition is satisfied, appears to follow the data closely. Our test also fails to reject that the population density satisfies Myerson's (1981) regularity condition.

In addition to the flexibility offered by the proposed framework, there are three other advantages of the proposed method. First, the notion of fidelity to the data that we optimize is independent of the constraint, including the choice of $\rho$ in the case of $\rho-$concavity constraints. Note that the objective function used in Koenker and Mizera's (2010) approach, Rényi entropy of order $\rho,$ is dependent on this constraint parameter. Also, adding a $\rho-$concavity constraint to the maximum likelihood estimator would not provide a convincing way to achieve this goal since this would not provide a convex optimization problem for values of $\rho<0$ and the estimator does not exist when $\rho<-1$ (Doss and Wellner, 2016).

Second, the shape constraints of the estimator only binds in regions in which it would not otherwise be satisfied. This is advantageous when using shape constraints that are not sufficiently restrictive to ensure parsimony by themselves. Moreover, the existence of the unconstrained minimizer ensures that the density estimator exists, regardless of the strength of the shape constraint. As discussed in the preceding paragraph, this is not the case for the maximum likelihood estimator. Although it is not our primary focus, we also provide an option for relying on the shape constraint for parsimony.

Third, the proposed algorithm solves an optimization problem over a set of variables that grows sub-linearly in the sample size, so the time complexity of the proposed algorithm compares favorably to other shape constrained density estimators.

The next section outlines the aspects of optimal transport that are required to formulate our estimator. Galichon (2016) provides a more comprehensive overview of the optimal transport literature, including its many applications in economics. The third section defines the estimator, provides the rate of convergence, and the asymptotic distribution of the estimator. This section also provides a test for the null hypothesis that the population density satisfies the shape constraint. The fourth section proves that the optimization problem defining the estimator is convex and provides an algorithm to calculate the estimator. Note that initializing this algorithm at a reasonable approximation of the density estimator provides a gain in computational efficiency, and an algorithm for finding an approximation is provided in the appendix. The sixth section generalizes the framework presented here to allow for density estimation and inference subject to a much larger class of shape constraints, with a focus on shape constraints that arise in the economics literature. Specifically, this section provides sufficient conditions for each of the results in the paper to hold under an arbitrary set of shape constraints. This generalization is used in the seventh section to provide evidence that the firms bidding on the California Department of Transportation's construction contracts have marginal cost distributions that satisfy Myerson's (1981) regularity condition.

A few notational conventions will be useful in the subsequent sections. For $x\in\mathbb{R}^{m},$ we will denote the vector with an $i^{th}$ element defined by $\exp(x_{i})$ as $\exp(x),$ and a similar convention will be used for $\log(x)$ and $x^{\rho}.$ Also, a diagonal matrix with a diagonal equal to the vector $x$ will be denoted by $D_{x},$ an $m\times1$ vector of ones by $\mathbf{1}_{m},$ the identity matrix by $I,$ element-wise division of the two vectors $x$ and $y$ by $x\oslash y,$ element-wise multiplication by $x\otimes y,$ the Moore-Penrose pseudoinverse of the matrix $A$ as $A^{+},$ the convolution between $f:\mathbb{R}^{d}\rightarrow\mathbb{R}^{1}$ and $g:\mathbb{R}^{d}\rightarrow\mathbb{R}^{1}$ by $g(x)\ast f(x),$ the derivative of $f(x)$ with respect to $x$ by $\triangledown_{x}f(x),$ and a Dirac delta function centered at $z$ by $\delta_{z}(x).$ Also, $sgn(x)$ will be used to denote a function that is $1$ when $x\geq0$ and $-1$ when $x<0.$ Since the proposed method requires discretizing densities, say $\mu:\mathcal{A}\rightarrow\mathbb{R}^{1}$ for $\mathcal{A}\subset\mathbb{R}^{d},$ we will denote the points in the mesh as $\{\mathbf{a}_{i}\}_{i=1}^{m},$ where $\mathbf{a}_{i}\in\mathbb{R}^{d}.$ We will also continue to include parenthesis after functions, as in $\mu(x)$ or $\mu(\cdot),$ and exclude parenthesis when denoting $\mu\in\mathbb{R}^{m}$ with elements $\mu_{i}=\mu(\mathbf{a}_{i}).$

Optimal Transport

Gaspard Monge formulated the theory of optimal transport in the $18^{th}$ century in order to derive the optimal method of moving a pile of sand to a nearby hole of the same volume. Specifically, suppose that both the pile of sand and the hole are defined on $\mathcal{A}\subset\mathbb{R}^{d},$ and we use the measures $\mathcal{M}_{0}:\mathcal{A}\rightarrow\mathbb{R}_{+}$ and $\mathcal{M}_{1}:\mathcal{A}\rightarrow\mathbb{R}_{+}$ to define the volume of the pile and the hole respectively. Monge sought to find a transportation plan, $T:\mathcal{A}\rightarrow\mathcal{A},$ that minimizes transportation costs while ensuring that the hole is completely filled.

Kantorovitch (1958) generalized this problem by describing the transportation plan by the absolutely continuous measure $\psi:\mathcal{A}\times\mathcal{A}\rightarrow\mathbb{R}_{+}.$ For example, given $a_{1}\in\mathcal{A}$ and $a_{2}\in\mathcal{A},$ we can view the Radon-Nikodym derivative of $\psi(\cdot),$ $d\psi(a_{1},a_{2}),$ as the amount of mass moved from $a_{1}$ to $a_{2}$ under the transportation plan, or coupling, $\psi(\cdot).$ Feasibility of $\psi(\cdot)$ simply requires $\psi(a,\mathcal{A})=\mathcal{M}_{1}(a)$ and $\psi(\mathcal{A},a)=\mathcal{M}_{0}(a)$ for all $a\in\mathcal{A}.$

When $\mathcal{M}_{i}(\cdot)$ is absolutely continuous, there exists $\mu_{i}(a)$ such that $\mathcal{M}_{i}(A)=\int_{A}\mu_{i}(a)da$ for each $A\subset\mathcal{A}$ by the Radon-Nikodym theorem. When $\mathcal{M}_{i}(\cdot)$ also satisfies $\mathcal{M}_{i}(\mathcal{A})=1,$ then $\mu_{i}(a)$ is also a probability density function. Optimal transport can be described without assuming $\mathcal{M}_{0}(\cdot)$ and $\mathcal{M}_{1}(\cdot)$ satisfy these conditions; however, since our goal is density estimation, we will generally restrict our attention to these cases for the rest of the paper. In addition we will define (or constrain) all density functions to be continuous, with the exception of $\delta_{z}(\cdot).$ We will use this notation to define $\Psi(\mu_{0}(\cdot),\mu_{1}(\cdot))$ as the set of feasible couplings.

The most common cost function in optimal transport is simply squared Euclidean distance. In this case the cost of moving one unit of earth from $a_{1}\in\mathcal{A}$ to $a_{2}\in\mathcal{A}$ is proportional to $\left\Vert a_{1}-a_{2}\right\Vert ^{2}.$ The resulting minimization problem is then given by

which we will refer to as the squared Wasserstein distance (Mallows, 1972). $W_{0}(\mu_{0}(\cdot),$ $\mu_{1}(\cdot))$ has many desirable properties, one being that $\sqrt{W_{0}(\mu_{0}(\cdot),\mu_{1}(\cdot))}$ satisfies all of the usual properties of a distance metric. $W_{0}(\mu_{0}(\cdot),\mu_{1}(\cdot))$ also metrizes weak convergence and convergence in the first two moments. In other words, given a sequence of densities, $\{\mu_{0}(\cdot),\mu_{1}(\cdot),...,\mu_{n}(\cdot)\},$ we have $\lim_{i\rightarrow\infty}W_{0}(\mu_{0}(\cdot),\mu_{i}(\cdot))=0$ if and only if $\mathcal{M}_{i}$ converges weakly to $\mathcal{M}_{0}$ and the first two moments of $\mu_{i}$ converge to the first two moments of $\mu_{0}.$

The Wasserstein distance between two distributions can be viewed as a measure of distance over the domain of the densities rather than in the direction of their range. For example, the Fréchet mean of two Dirac delta functions, centered at $a$ and $b,$ in the spaces of densities equipped with an $L^{2}$ norm is given by $\text{arg min}_{\nu(x)}\;\left\Vert \delta_{a}(x)-\nu(x)\right\Vert ^{2}+\left\Vert \delta_{b}(x)-\nu(x)\right\Vert ^{2}=\delta_{a}(x)/2+\delta_{b}(x)/2,$ while a similar notion of average in the spaces of densities equipped with the Wasserstein distance is $\text{arg min}_{\nu(x)}\;W_{0}(\delta_{a}(x),\nu(x))+W_{0}(\delta_{b}(x),\nu(x))=\delta_{a/2+b/2}(x).$ To make this intuition more explicit, when $\mathcal{A}\subset\mathbb{R}^{1},$ one can show $W_{0}(\mu_{0}(\cdot),\mu_{1}(\cdot))$ can also be expressed as $\int_{0}^{1}(Q_{0}(\tau)-Q_{1}(\tau))^{2}d\tau,$ where $Q_{0}(\tau)$ and $Q_{1}(\tau)$ are the quantile functions corresponding to $\mu_{0}(\cdot)$ and $\mu_{1}(\cdot)$ respectively. Thus, $(Q_{0}(\tau)-Q_{1}(\tau))^{2}$ represents a squared distance between two points in $\mathcal{A}$ (Villani, 2003).

In practice augmenting the Wasserstein distance with a regularization term ameliorates some numerical difficulties, which will be described below in more detail. The regularized squared Wasserstein distance is a generalization of $W_{0}(\mu_{0}(\cdot),\mu_{1}(\cdot)),$ and is defined by

where $\gamma\geq0$ and $H(\psi(\cdot)):=-\int_{\mathcal{A}\times\mathcal{A}}\log\psi(a_{1},a_{2})d\psi(a_{1},a_{2})$ is the Shannon entropy of $\psi(\cdot)$ (Cuturi, 2013; Cuturi and Doucet, 2014). In the shape constrained density estimation setting, the addition of this entropy term is advantageous for several reasons. First, the objective function is strictly convex when $\gamma>0,$ so the optimal coupling will always be unique. Second, in practice $\mathcal{A}$ must be discretized before finding the unregularized Wasserstein distance, and the computational cost of solving for the optimal coupling scales at least cubically in the number of points in the mesh. Third, after discretizing, the minimizer of $W_{\gamma}(\mu_{0},\mu_{1})$ with respect to $\mu_{0}$ is often a more accurate representation of the minimizer of $W_{0}(\mu_{0}(\cdot),\mu_{1}(\cdot)),$ when $\gamma$ is set to a reasonably small value.\footnote{ Minimizing $W_{0}(\mu_{0},\mu_{1})$ with respect to $\mu_{0}$ generally results in a minimizing density with many large discrete changes. For more detail, see Figures 3.1, 3.2, and the accompanying explanation in (Cuturi and Peyré, 2016). } Four, using $W_{\gamma}(\mu_{0}(\cdot),\mu_{1}(\cdot))$ allows us to avoid assumptions in the next section regarding the existence of the second moments of $\mu_{0}$ and $\mu_{1}.$ Lastly, we can find the minimizer of (2) with a very computationally efficient algorithm after discretizing, which we will describe next.

To introduce the discretized counterparts of $d\psi(\cdot),\mu_{1}(\cdot),$ and $\mu_{0}(\cdot),$ recall our uniform mesh over $\mathcal{A}$ contains the vertices $\{\mathbf{a}_{i}\}_{i=1}^{m},$ and let $\mu_{0},\mu_{1}$ define $\mu_{0}(\mathbf{a}_{i}),\mu_{1}(\mathbf{a}_{i})$ respectively. Also, let $M_{m\times m}$ so that $M_{ij}:=\left\Vert \mathbf{a}_{i}-\mathbf{a}_{j}\right\Vert ^{2},$ and $\psi_{m\times m}$ so that $\psi_{i,j}:=d\psi(\mathbf{a}_{i},\mathbf{a}_{j}).$ After discretizing, (2) can be written as

The corresponding Lagrangian is given by

and the first order conditions imply

In other words, there exists $v,w\in\mathbb{R}_{+}^{m}$ such that the optimal coupling has elements $\psi_{ij}=K_{ij}w_{i}v_{j},$ where $K_{ij}:=\exp\left(-\left\Vert \mathbf{a}_{i}-\mathbf{a}_{j}\right\Vert ^{2}/\gamma\right),$ a symmetric $m\times m$ matrix. This can also be written as,

so adding the entropy term to the objective function reduces the dimensionality of the optimization problem from $m^{2}$ to $2m.$ Sinkhorn (1967) shows that $\psi$ is unique. The iterative proportional fitting procedure (IPFP) is an efficient method of computing these vectors; see Krupp (1979). This method iteratively redefines $w$ so that $D_{w}Kv=\mu_{0},$ and subsequently redefines $v$ so that $D_{v}Kw=\mu_{1},$ as summarized in Algorithm 1. Note that after combining these equalities we have $D_{w}K(\mu_{1}\oslash(Kw))=\mu_{0},$ which will be used in subsequent sections.

In the rest of the paper we will make substantial use of the dual of (3)-(5). Cuturi and Doucet (2014) show that the dual is given by the unconstrained optimization problem,

Note that $K=\exp\left(-M/\gamma\right),$ so the first order conditions of (8) can be written as,

Note that after replacing $\exp\left(x/\gamma\right)$ with $w$ and $\exp\left(y/\gamma\right)$ with $v,$ these formulas are equivalent to the updates of $w$ and $v$ given in Algorithm 1. Also, given $x$ and $y$ that satisfy (9)-(10), consider the vectors $\tilde{x}:=x-c$ and $\tilde{y}:=y+c,$ where $c\in\mathbb{R}.$ Since $\mu_{0}$ and $\mu_{1}$ have the same sum, the objective function of (8), evaluated at $\tilde{x}$ and $\tilde{y}$ must equal the objective function evaluated at $x$ and $y.$ Since $\exp\left(\tilde{y}_{j}/\gamma\right)=\exp\left(y_{j}/\gamma\right)\exp\left(c/\gamma\right)$ and $\exp\left(\tilde{x}_{i}/\gamma\right)=\exp\left(x_{i}/\gamma\right)\exp\left(-c/\gamma\right),$ $\tilde{x}$ and $\tilde{y}$ must also satisfy (9)-(10). In other words, while $v$ and $w$ are unique up to a multiplicative constant on $w$ and one over this constant on $v,$ $y$ and $x$ are unique up to the additive constant $c.$

A few comments regarding the effect of $\gamma$ on the optimal coupling will also be useful in subsequent sections. Higher values of $\gamma$ correspond to placing a higher penalty on the negative entropy of the coupling, so the optimal coupling becomes more dispersed as this parameter is increased. Also, in the limit $\gamma\rightarrow0,$ $W_{\gamma}(\mu_{0},\mu_{1})$ converges to $W_{0}(\mu_{0},\mu_{1})$ at the rate $O(\exp(-1/\gamma))$ and $\psi$ converges to the optimal unregularized coupling at this same rate; see Benamou et al. (2015) and Cuturi (2013). In the next section will use $W_{\gamma}(\cdot)$ as an objective function to define the proposed estimator and show how $\gamma$ impacts this estimator in more detail.

Shape-Constrained Density Estimation

The input of the density estimator proposed in this paper is a kernel density estimator, $\mu,$ based on $N$ i.i.d. datapoints, $\{z_{i}\}_{i=1}^{N},$ drawn from a uniformly continuous population density, $\mu^{\star}:\mathbb{R}^{d}\rightarrow\mathbb{R}_{+}^{1},$ with a bandwidth of $\sigma\geq0.$ Our results also hold when $\gamma$ and $\sigma$ are redefined to be functions of the form $c_{1}(\{z_{i}\}_{i})\gamma$ and $c_{2}(\{z_{i}\}_{i})\sigma,$ where $c_{j}(\{z_{i}\}_{i})\stackrel{p}{\rightarrow}c_{j}$ for $j\in\{1,2\}$ at any polynomial rate. In the interest of the brevity of notation, we will only write the parameters in this way when their randomness would have a non-trivial impact on the result. $\gamma,\sigma$ and $m$ will also be dependent on $N,$ but we will suppress this input throughout the paper and discuss our recommendations for defining them after Theorem 2.

With this notation in mind we can define the shape-constrained density estimator as,

where $\mathcal{K}$ is the cone of concave functions. Although $\mathcal{K}$ is a convex set, the set of $\rho-$concave densities is generally not convex. To ameliorate this problem, we will use a similar formulation as Koenker and Mizera (2010) and solve for $g:=f^{\rho}.$ Note that no generality is lost in doing so, as there is a one to one correspondence between $g$ and $f.$

For the clarity of the derivations in the next section, we will also define $g$ to be a vector of length $m-1$ and refer to the element of the vector $f,$ of length $m,$ that corresponds to this omitted value as $f_{k}.$ We will set $f_{k}$ so that the density sums to $m.$ In other words, the elements of the density $f$ that do not correspond to this $k^{th}$ element, denoted by $f_{-k},$ will be set equal to $g^{1/\rho}$ and $f_{k}$ will be set equal to $m-\sum_{i}g_{i}^{1/\rho}.$

Unlike optimizing over $g$ rather than $f,$ this can be seen as a slight modification of our initial optimization problem, since we will also exclude the shape constraints that depend on the $k^{th}$ element of $f.$ However, as discussed in the next section in more detail, one can choose $k$ to correspond to an element on the boundary of the mesh so that the estimator satisfies the shape constraint everywhere on the interior of its domain.

In a slight abuse of notation, we will also denote the objective function as $W_{\gamma}(g^{1/\rho},\mu).$ In summary, we will define $W_{\gamma}(g^{1/\rho},$ $\mu)$ by,

and the final form of our optimization problem is,

where $\alpha_{i}\in\mathbb{R}^{1},$ $\beta_{i}\in\mathbb{R}^{d},$ and $d$ is the dimension of the support of $\mu$ and $f.$ These inequality constraints are used by Afriat (1972) to estimate production functions with concavity constraints. They tend to provide a gain in numerical accuracy relative to local concavity constraints.

The following Lemma provides the limiting distribution of the estimator after removing the shape constraints. This is achieved by showing that the resulting density can be viewed as a kernel density estimator with a bandwidth of $\sqrt{\sigma^{2}+\gamma/2}$ away from the edges of the mesh over $\mathcal{A}.$ To avoid these boundary value effects, we recommend enlarging the domain of $\mu$ and $\hat{f}$ to include regions within approximately $3\sqrt{\sigma^{2}+\gamma/2}$ of the datapoints. Alternatively, one could replace the matrix $K$ with the direct application of a Gaussian filter. This approach has the added benefit of reducing the computational complexity of Algorithm 1 to $O(m\log m),$ so this is our recommended approach when $d>2;$ see Solomon et al. (2015) for more details on this method.

$\vphantom{c}$

Lemma 1: Suppose $\mu$ is a kernel density estimate, generated with a Gaussian kernel and a bandwidth $\sigma$ and $\mu^{\star}(\cdot)$ is uniformly continuous. Also, suppose $\sigma,$ $\gamma$ and $m$ are chosen so that $\sqrt{\sigma^{2}+\gamma/2}\stackrel{p}{\rightarrow}0,$ $N\sqrt{\sigma^{2}+\gamma/2}\stackrel{p}{\rightarrow}\infty,$ and $\min_{i\neq j}\left\Vert {\bf a}_{i}-\mathbf{a}_{j}\right\Vert /\sqrt{\gamma}\rightarrow0$ as $N\rightarrow\infty.$ Then there exists $c\in\mathbb{R}^{1}$ so that the limiting density of $f_{unc}:=\text{arg min}_{f}\;W_{\gamma}(f,\mu)$ is given by,\footnote{When $\mu^{\star}(\cdot)$ is differentiable at $\mathbf{a}_{i},$ $c=\triangle\mu_{i}^{\star}(x)/2\mid_{x=\mathbf{a}_{i}},$ where $\triangle\mu_{i}^{\star}(x)$ denotes the Laplacian, $\sum_{j=1}^{d}\triangledown_{x_{j},x_{j}}\mu^{\star}(x).$ However, $c$ can be defined without assuming $\mu^{\star}(\cdot)$ is differentiable; Karunamuni and Mehra (1991) provide more details on this approach. }

Proof: To find $f_{unc},$ consider the optimization problem,

The corresponding Lagrangian is

and the first order conditions imply that there exists $v\in\mathbb{R}^{m}$ such that

Note that convexity of negative entropy implies that the optimal coupling will correspond to a minimizer. After combining this equality with the constraints, we have

which implies

Let $\kappa:=K\mathbf{1}.$ Now we can find $f_{unc}$ by finding the other marginal of $\psi,$

Let $\phi_{\eta}:\mathbb{R}^{d}\rightarrow\mathbb{R}^{1}$ be a Gaussian density with variance $\eta I_{d}$ and mean $\mathbf{0}_{d}.$ Suppose $\nu:\mathbb{R}^{d}\rightarrow\mathbb{R}^{1}$ is a continuous function. Then,

by the definition of $K,$ so $K$ can be viewed as a linear operator that discretizes the convolution $\nu(\cdot)\ast\phi(\cdot)\sqrt{\pi}/\gamma.$ Since the distance between adjacent points decreases at a rate that is faster than $\sqrt{\gamma/2},$ and the convex hull of the data is a strict subset of the convex hull of $\{\mathbf{a}_{i}\}_{i},$ the error of this approximation converges to zero on the convex hull of the data. Thus, as $N$ diverges we have $\kappa_{i}\rightarrow\sqrt{\pi}/(m\gamma)$ for all $i$ in the convex hull of the data.

Since $\mu$ is a kernel density estimator, it can be expressed as $\mu=\left(\sum_{i}\delta_{z_{i}}(\cdot)\ast\phi_{\sigma^{2}}/N\right),$ $(\mathbf{a})$ which implies that $f_{unc}$ converges to $\left(\sum_{i}\delta_{z_{i}}(\cdot)\ast\phi_{\sigma^{2}}\ast\phi_{\gamma/2}/N\right)(\mathbf{a}).$ The convolution of two Gaussian densities is $\phi_{\sigma^{2}}\ast\phi_{\gamma/2}=\phi_{\sigma^{2}+\gamma/2},$ so $f_{unc}=\left(\sum_{i}\delta(z_{i})\right)\ast\phi_{\sigma^{2}+\gamma/2}(y)/N,$ which defines a kernel density estimator with a bandwidth of $\sqrt{\sigma^{2}+\gamma/2},$ and Parzen (1962) provides the limiting density of the kernel density estimator when $\sqrt{\sigma^{2}+\gamma/2}\rightarrow0$ and $N\sqrt{\sigma^{2}+\gamma/2}\rightarrow\infty.$

The following theorem provides the limiting density of the estimator with the primary additional assumption that $(\mu^{\star}(x))^{\rho}$ is strictly concave. Afterward, we will move onto results that relax this assumption.

$\vphantom{c}$

Theorem 2: Suppose the assumptions from Lemma 1 hold. If $\mu^{\star}$ is in the interior of the feasible set, $N\sqrt{c_{2}(\{z_{i}\}_{i})\sigma^{2}+c_{1}(\{z_{i}\}_{i})\gamma/2}/\log(N)\rightarrow\infty,$ and for $j\in\{1,2\},$ $c_{j}(\{z_{i}\}_{i})$ converges in probability to a constant, then there exists $c\in\mathbb{R}^{1}$ so that,

Proof: Since $W_{\gamma}(\cdot)$ is differentiable, we can apply the envelope theorem to the dual problem (8) to show that $\triangledown_{f}W_{\gamma}(f,\mu)=x.$ This implies that the $\tilde{f}\in\mathbb{R}^{m}$ is a critical point of $f\mapsto W_{\gamma}(f,\mu)$ if and only if it has a corresponding dual variable in (8) of $x=\mathbf{0},$ so $f_{unc}$ is the unique minimizer of $W_{\gamma}(\cdot)$ by the proof of Lemma 1 and the definition $w:=\exp(x/\gamma).$\footnote{Uniqueness of $f_{unc}$ can also be shown using strict convexity of $W_{\gamma}(f,\mu)$ in $f,$ which will be shown in Theorem 6. } This implies $\hat{f}=f_{unc}$ when $f_{unc}$ is feasible. Einmahl and Mason (2005) show that $f_{unc}\stackrel{a.s.}{\rightarrow}\mu^{\star}$ under the additional assumptions of the theorem when using data dependent bandwidths, as in the statement of the theorem. The result follows from the fact that $\mu^{\star}$ is in the interior of the feasible set.

Remark 1: Algorithm 1 can be slow to converge when $\gamma$ is chosen to be too small, but our assumption that $\gamma\rightarrow0$ as $N\rightarrow\infty$ is not problematic when the assumptions of the theorem hold. To see this note that $f_{unc}$ can be written as $K(\mu\oslash(K\mathbf{1})).$ Evaluating Algorithm 1 at input densities $\mu$ and $f_{unc}$ would result in the algorithm first initializing $w$ as $\mathbf{1}_{m}$ and then define $v$ to be $\mu\oslash(K\mathbf{1}).$ The $w-$update would redefine $w$ to be $f_{unc}\oslash K(\mu\oslash(K\mathbf{1})),$ but since this is also equal to $\mathbf{1},$ the algorithm has already converged. For input densities $f$ and $\mu,$ it is also generally the case that Algorithm 1 tends to converge at a faster rate for densities that are closer to $f_{unc}.$

Remark 2: Our more general convergence result (Theorem 4) provides the same rate of convergence as Theorem 2 without requiring that $\mu^{\star}$ is an interior point of the feasible set, so we can compare this rate with the corresponding rate of convergence of shape constrained maximum likelihood estimators. Seregin and Wellner (2010) show that the pointwise minimax absolute error loss of the $\log-$concave constrained maximum likelihood estimators at $x\in\mathcal{A}$ is $N^{-2/(d+4)}$ when $\mu^{\star}(x)$ is twice differentiable, $x$ is in the interior of the domain of $\mu^{\star}(\cdot),$ and the Hessian has full rank. Theorem 2 implies that our estimator also obtains these bounds when $\sqrt{\sigma^{2}+\gamma/2}$ is chosen so that it converges to zero at the optimal rate, in the sense of minimizing mean squared error, of $O_{p}(N^{-1/(d+4)}).$

In practice, $\sqrt{\sigma^{2}+\gamma/2}$ has a more noticeable impact on the resulting density estimator than the individual values of $\sigma$ and $\gamma,$ so we recommend fixing $\gamma/\sigma^{2}$ to be a constant and focusing on the choice of $\sqrt{\sigma^{2}+\gamma/2}.$ Since increasing $\gamma/\sigma^{2}$ tends to result in an increase in the rate of convergence and numerical stability of Algorithm 1, our recommendation is $\gamma/\sigma^{2}=8,$ which was used in all of the applications and examples given below.

From a numerical standpoint, it is possible to set $\sqrt{\sigma^{2}+\gamma/2}$ so that the smoothing provided by $\gamma$ and $\sigma$ is negligible and parsimony is almost entirely ensured by the shape constraints. Since the stability of Algorithm 1 is the only limiting factor on how small we can make $\sqrt{\sigma^{2}+\gamma/2},$ we can use this fact to estimate a lower bound on the value $\sqrt{\sigma^{2}+\gamma/2}.$ Specifically, we can estimate the lowest possible value of $\sqrt{\sigma^{2}+\gamma/2}$ that still results in convergence of Algorithm 1 to within a given tolerance in a fixed number of iterations, say 100, which can be achieved using a root finding algorithm. This requires an approximation of the density estimate, and Appendix B provides a method for finding an approximation in a computationally efficient manner. When running the main optimization algorithm described in the next section, we recommend increasing $\gamma$ after this root is found by approximately one fourth and increasing the number of iterations used in Algorithm 1 by a factor of approximately ten to ensure the accuracy of the Hessian.

In our tests this approach resulted in density estimates that were surprisingly similar to estimates using the method described by Koenker and Mizera (2010). However, it is likely that the condition $N\sqrt{\sigma^{2}+\gamma/2}\rightarrow\infty$ would not hold in this case, which is a requirement of Theorem 2. If $N\sqrt{\sigma^{2}+\gamma/2}\rightarrow0$ and we redefine the location of the vertices in the mesh, $\mathbf{a},$ to be equal to the datapoints, then asymptotically $\mu_{i}$ and $f_{unc,i}$ would only be influenced by a single datapoint. The following theorem shows that we can at least ensure $\hat{f}\overset{p}{\rightarrow}\mu^{\star}$ in the case $\gamma=\sigma=0.$

$\vphantom{c}$

Theorem 3: Suppose the locations of the vertices in the mesh $\mathbf{a}$ are defined to be the datapoints, $\mu^{\star}(z)>0$ on its bounded domain $\mathcal{A},$ and $\mu_{i}=1/N$ for $i\in\{1,2,...,N\}.$ If there exists $\epsilon>0$ such that $\int_{\mathcal{A}}\left|x\right|^{2+\epsilon}d\mu^{\star}<\infty,$ $\mu^{\star}$ is $\rho-$concave, and $\gamma=\sigma=0,$ then $\hat{f}\overset{p}{\rightarrow}\mu^{\star}.$

Proof: The Wasserstein distance between density functions, denoted by $W_{0}(\nu_{1}(z),$ $\nu_{2}(z)),$ is commonly approximated by discretizing, using $\nu_{j,i}=\nu_{j}(z_{i})$ for $j\in\{1,2\}$ and $z_{i}\in\{\mathbf{a}_{i}\}_{i},$ and solving a linear programming problem, as discussed in the previous section. As long as $\min\{r\mid\mathcal{A}\subset\cup_{i=1}^{m}B_{r}(\mathbf{a}_{i})\}$ goes to zero asymptotically, $W_{0}(\nu_{1},\nu_{2})$ converges to its non-discretized counterpart, $W_{0}(\nu_{1}(z),\nu_{2}(z));$ see for example, (Cuturi and Peyré, 2016). Since $\mu^{\star}(z)>0$ for all $z\in\mathcal{A}$ and $\{\mathbf{a}_{i}\}_{i}=\{z_{i}\}_{i},$ this condition holds asymptotically, so showing $W_{0}(\hat{f}(z),\mu^{\star}(z))\stackrel{p}{\rightarrow}0,$ where $\hat{f}(z)$ is the linear interpolation of $\{(\mathbf{a}_{i},\hat{f}_{i})\}_{i=1}^{N},$ will imply the result.

Let $\mu(z)=\sum_{i}\delta_{z_{i}}(z)/N$ and $f(z)$ be an arbitrary proper density with $\int_{\mathcal{A}}\left|x\right|^{2+\epsilon}df.$ After using the triangle inequality twice, we can bound $W_{0}(f(\cdot),\mu(\cdot))$ by $W_{0}(f(\cdot),\mu^{\star}(\cdot))$ $\pm W_{0}(\mu(\cdot),\mu^{\star}(\cdot))$ for all $N>1.$ Note that Wasserstein distances metrize weak convergence, in the sense that $W_{0}(\nu_{0}(z),\nu_{N}(z))\stackrel{p}{\rightarrow}0$ if and only if the distribution corresponding to $\nu_{N}(z)$ converges weakly to that of $\nu_{0}(z)$ and the first two moments of $\nu_{N}(z)$ converge to the first two moments of $\nu_{0}(z);$ see for example, Theorem 7.12 of (Villani, 2003). Thus, $W_{0}(\mu(\cdot),\mu^{\star}(\cdot))\stackrel{p}{\rightarrow}0$ by the result of Kiefer and Wolfowitz (1956), so our bounds on $W_{0}(f(\cdot),\mu(\cdot))$ imply that $W_{0}(f(\cdot),\mu(\cdot))\stackrel{p}{\rightarrow}W_{0}(f(\cdot),\mu^{\star}(\cdot))$ for all such $f(\cdot).$

Note that $W_{0}(f(\cdot),\mu^{\star}(\cdot))$ obtains its minimum value with respect to $f(\cdot)$ at densities with distribution functions that are arbitrarily close to the distribution function of $\mu^{\star}(\cdot).$ Since $\mu^{\star}$ is $\rho-$concave by definition, and $\hat{f}(\cdot)$ is also by our constraint, this distribution has well defined and continuous density function by Aleksandrov's theorem. Thus, $\hat{f}(\cdot)\overset{p}{\rightarrow}\mu^{\star}(\cdot)$ by the standard argument for consistency of an M-estimator; see for example, Theorem 5.7 of (van der Vaart, 2000).

For the five reasons discussed in the previous section, our primary focus is on cases in which $\mu,$ and $f_{unc},$ are also consistent estimators, so we will move onto our recommended approach for setting $\sqrt{\sigma^{2}+\gamma/2}$ and $m.$ Like kernel density estimators, the mean squared error of $f_{unc}$ is minimized when $\sqrt{\sigma^{2}+\gamma/2}$ is $O(N^{-1/(d+4)}).$ Also, any of the standard techniques for setting the bandwidth of kernel density estimators are reasonable methods of setting $\sqrt{\sigma^{2}+\gamma/2},$ including cross validation or a rule-of-thumb bandwidth estimator; for examples of rule-of-thumb bandwidth estimators see (Silverman, 1986; Scott, 1992). Since the shape constraint also helps to ensure the parsimony of the density estimate, these rules should generally be regarded as upper bounds on $\sqrt{\sigma^{2}+\gamma/2}.$ In practice, using Scott's (1992) rule-of-thumb multiplied by $2/3$ works well. After dividing each dimension of the dataset, $\{z_{i,j}\}_{i}$ for $j\in\{1,2,...,d\},$ by $\min(\text{IQR}(\{z_{i,j}\}_{i})/1.349,\hat{\sigma}(\{z_{i,j}\}_{i})),$ combining this with $\gamma/\sigma^{2}=8,$ results in $\sigma\approx N^{-1/(d+4)}/3$ and $\gamma\approx N^{-2/(d+4)}4/5.$

The choice of $m$ is likely less critical than that of $\sqrt{\sigma^{2}+\gamma/2}.$ The effect of $m$ on the estimator can be viewed as analogous to the effect of the size of a Gaussian filter on its output. In practice, these filters are often constructed by discretizing a Gaussian density over a mesh with a resolution equal to the standard deviation. Using this as a lower bound in our setting yields, $m\sqrt{\gamma/2}/d\geq1.$ The gain in accuracy from increasing $m$ beyond the point $m\sqrt{\gamma/2}/d=1$ generally diminishes fairly quickly when $m$ is set to larger values, so we recommend choosing $m$ so that $m=\left\lceil N^{\beta}d/\sqrt{\gamma/2}\right\rceil ,$ where $\left\lceil \cdot\right\rceil $ is the ceiling function, for any $\beta>0.$ In the interest of specificity, we recommend using $m=\left\lceil N^{1/5}d/\sqrt{\gamma/2}\right\rceil .$

If the set of constraints that bind asymptotically is known, it is straightforward to find the asymptotic distribution of $\hat{f}$ when $\mu^{\star}$ is not $\rho-$concave; however, this is rarely the case in practice. For this reason, finding confidence intervals for estimators subject to shape constraints is an active area of research in nonparametric econometrics. The difficulties in these cases are due to the estimators, when viewed as functions of the data, not being sufficiently smooth (either differentiable or Hadamard differentiable) at points in which an inequality constraint goes from slack to active, and these notions of smoothness are requirements of many of the commonly used methods for deriving limiting distributions (Andrews, 2000).

The following theorem establishes the limiting distribution of the estimator conditional on the set of constraints that bind asymptotically. Note that the set of active constraints converges when $\sqrt{\sigma^{2}+\gamma/2}$ converges to zero at its optimal rate. Thus, the limiting distribution part of the proof might be useful with a rather large sample size. However, since the set of constraints that binds in a finite sample may not be equal to the set of constraints that bind asymptotically, in many cases the primary value of this result is that it provides the limiting value of the estimator without assuming that $(\mu^{\star})^{\rho}$ is strictly concave.

There are several possible methods to avoid conditioning on the set of active constraints, which is an area of ongoing research. One interesting possibility would be to adopt a similar method as Horowitz and Lee (2017) to the present setting. Their technique involves explicitly defining the set of active constraints as those constraints that would bind otherwise or that would “nearly bind.” Asymptotically correct coverage follows from consistency of this conservative estimate of the active set.

$\vphantom{c}$

Theorem 4: Suppose the assumptions from Lemma 1 hold. In addition, suppose $\sigma N\rightarrow\infty.$ Then, conditional on the set of constraints that are active asymptotically, say $\Omega,$ there exists $G,$ as defined in (23-25), and $c\in\mathbb{R}^{m}$ such that

where $\tilde{\mu}$ is the Wasserstein projection of $\mu^{\star}$ onto the set of feasible densities. Specifically, if $\tilde{g}$ is the minimizer of

then $\tilde{\mu}_{-k}:=\tilde{g}^{1/\rho}$ and $\tilde{\mu}_{k}=m-\sum_{i}\tilde{g}_{i}^{1/\rho}$ in the case of a $\rho-$concavity constraint.

Proof: We only consider the case in which $d=1$ and $k=m$ for the sake of clarity, but generalizing the proof to higher dimensions is straightforward. Since $d=1,$ we will assume all sets containing the numerical labels of vertices in the mesh are ordered sequentially. For example, $\Omega_{i}$ will be viewed as the $i^{th}$ largest vertex label in which the constraint binds. Let $\dot{\Omega}:=\{i\mid i\notin\Omega\land i\pm1\in\Omega\},$ the boundary of the complement of $\Omega.$

Since $\sqrt{\sigma^{2}+\gamma/2}\rightarrow0,$ we have $\sigma\rightarrow0.$ Since this is the case, and since $\sigma N\rightarrow\infty$ and $\mu^{\star}$ is continuous, $\mu$ converges uniformly to $\mu^{\star}$ (Parzen, 1962). We can view $\hat{f}$ as a continuous function of $\mu,$ so the continuous mapping theorem implies $\hat{f}$ converges to $\tilde{\mu}.$

Let the matrix $A_{\left|\Omega\right|\times(\left|\Omega\right|+\left|\dot{\Omega}\right|)}$ be defined so that $Ag_{\Omega\cup\dot{\Omega}}$ is the second difference of $g_{\Omega\cup\dot{\Omega}},$ so the binding constraints can be denoted by $Ag_{\Omega\cup\dot{\Omega}}\leq\mathbf{0}.$\footnote{In the interest of the simplicity of the exposition, we are defining $A$ using local concavity constraints rather than Afriat's (1972) formulation of concavity constraints. The nonzero elements of each row of $A$ are given by $(A_{i,j-1},A_{i,j},A_{i,j+1})=sgn(\rho)(1,-2,1).$ Note that $j>1,$ and, since the sets $\Omega$ and $\Omega\cup\dot{\Omega}$ are ordered by the vertex labels, this matrix is in row echelon form by construction. Thus, $A$ has full rank, so $AA^{T}$ is nonsingular. } Since we assumed that $\Omega$ is known, $\hat{g},$ the vector of Lagrange multipliers, $\lambda,$ and $x$ are defined by the following system of equations,

where the final equality results from combining the first order conditions of the optimization problem defining $W_{\gamma}(g^{1/\rho},\mu),$ which are given by (9) and (10).

Theorem 5 will establish that $W_{\gamma}(g^{1/\rho},\mu)$ is strictly convex in $g.$ Since this is the case, the solution of this system is unique. The implicit function theorem, applied to (19)-(22), implies that we can view $\hat{g},$ and thus also $\hat{f},$ as a differentiable function of $\mu.$ After simplifying this system, we will find this limiting density using the delta method.

We will begin by deriving $\triangledown_{\mu}x$ using (22). Recall $w:=\exp(x/\gamma),$ and let $h_{1}(x,\mu)$ be defined as,

(22) can be written as $h_{1}(x,\mu)=\mathbf{0}.$ Recall the following equalities from the previous section: $\psi=D_{w}KD_{v},$ $f=D_{w}Kv,$ and $\mu=D_{v}Kw.$ These imply

and

The implicit function theorem implies,

Since $x=\gamma\log(w),$ we have,

(21) implies that each element in $\hat{g}_{\Omega}$ can be expressed as a mean of its neighbors, so we can express all of the elements of $\hat{g}_{\Omega}$ as a weighted mean of $\hat{g}_{\dot{\Omega}}.$ In other words, there exists $C$ such that $\hat{g}_{\Omega\cup\dot{\Omega}}=C\hat{g}_{\dot{\Omega}}.$ (19) implies that, given $x$ and $g_{\dot{\Omega}},$ $\lambda$ is given by,

We will use the additional $\left|\Theta\right|$ equations in (19) to define $\hat{f}_{\dot{\Omega}}.$ After replacing each instance of $\hat{g}_{\dot{\Omega}}$ with $\hat{f}_{\dot{\Omega}}^{\rho}$ and using this definition of $\lambda,$ we can write (19) as $h_{2}(x,\hat{f}_{\dot{\Omega}})=\mathbf{0},$ where $h_{2}(x,\hat{f}_{\dot{\Omega}})$ is defined as

and $\tilde{A}:=A_{\cdot,\dot{\Omega}}^{T}\left(AA^{T}\right)^{-1}A.$ This implies

and

where $B_{\left|\dot{\Omega}\right|\times\left|\Omega\cup\dot{\Omega}\right|}$ is defined so that $B_{i,j}$ is equal to $\hat{f}_{\dot{\Omega}}^{1-\rho}$ when $i,j$ satisfies $\dot{\Omega}_{i}=\{\Omega\cup\dot{\Omega}\}_{j}$ and zero otherwise. The implicit function theorem implies

$\triangledown_{x_{\Omega\cup\dot{\Omega}}-x_{k}}\hat{f}_{\dot{\Omega}}=-\left(D_{(x_{\dot{\Omega}}-x_{k})}\hat{f}_{\dot{\Omega}}^{-\rho}+\tilde{A}D_{(x_{\Omega\cup\dot{\Omega}}-x_{k})\otimes(C\hat{f}_{\dot{\Omega}}^{\rho})^{1/\rho-2}}CD_{\hat{f}_{\dot{\Omega}}^{\rho-1}}\right)^{-1}$

so

$-\left(D_{(x_{\dot{\Omega}}-x_{k})}\hat{f}_{\dot{\Omega}}^{-\rho}+\tilde{A}D_{(x_{\Omega\cup\dot{\Omega}}-x_{k})\otimes(C\hat{f}_{\dot{\Omega}}^{\rho})^{1/\rho-2}}CD_{\hat{f}_{\dot{\Omega}}^{\rho-1}}\right)^{-1}$

where $P_{\left|\Omega\cup\dot{\Omega}\right|\times m}:=\triangledown_{x}(x_{\Omega\cup\dot{\Omega}}-x_{k})$ is defined so that $P_{i,j}$ is equal to $1$ when $i,j$ satisfies $\{\Omega\cup\dot{\Omega}\}_{i}=j,$ $-1$ when $j=m,$ and zero otherwise.

Note that (20) implies $x_{i}=x_{k}$ for all $i\in\{j\mid j\notin\Omega\cap j\notin\dot{\Omega}\}.$ The proof of Theorem 2 shows that this condition defines $f_{unc,i},$ so we have $f_{unc,i}=\hat{f}_{i}$ for all such $i.$ This implies,

Lastly, writing the equations defining $\hat{f}_{\Omega}$ in terms of $\hat{f}_{\dot{\Omega}}$ yields $\hat{f}_{\Omega}=(C\hat{f}_{\dot{\Omega}}^{\rho})^{1/\rho},$ so we have,

There are also a few options for testing if a population density satisfies a shape constraint. Hypothetically, one could consistently test the null hypothesis that a population density is $\rho-$concave using any consistent shape constrained density estimator. This can be done by estimating the population distribution subject to the shape constraint and then using one of the classic tests for whether or not the empirical distribution of the data is equal to this estimate; see for example (Smirnov, 1948; Anderson and Darling, 1952). In these cases, choosing a test with a statistic that is closely related to the fidelity criterion used for estimation allows for a more straightforward interpretation of the result. For example, if the test statistic is equal to the fidelity criterion that the estimator optimizes, we would reject the null if and only if we would also reject the null for every density that satisfies the shape constraint.

The following theorem provides the distribution of $W_{\gamma}(f_{unc},\mu)$ and a consistent test for the null hypothesis that $\mu^{\star}(x)$ satisfies the shape constraint based on this distribution. The method also has the straightforward interpretation from the preceding paragraph, so, if we denote the set of $\rho-$concave densities by $\mathcal{K}_{\rho},$ the null is rejected $\min_{f\in\mathcal{K}_{\rho}}W_{\gamma}(f,\mu)-W_{\gamma}(f_{unc},\mu)$ is statistically significant. Since $W_{\gamma}(f,\mu)$ is differentiable in $f,$ this can be achieved without conditioning on the set of active constraints.

$\vphantom{c}$

Theorem 5: Suppose the assumptions from Lemma 1 hold and that $\sigma N\rightarrow\infty.$ Let $T$ be defined as $N(\sigma^{2}+\gamma/2)^{d/2}\left(W_{\gamma}(\mu^{\star},\mu)-W_{\gamma}(f_{unc},\mu)\right),$ $\psi$ as the optimal coupling between $\mu$ to $f_{unc},$ $\psi^{\star}$ as the optimal coupling between $\mu^{\star}$ and itself, $B^{\star}$ as $\gamma(D_{\mu^{\star}}-\psi^{\star}D_{\mathbf{1}\oslash\mu^{\star}}\psi^{\star T})^{+}/2,$ and $B$ as $\gamma(D_{f_{unc}}-\psi D_{\mathbf{1}\oslash\mu}\psi^{T})^{+}/2.$ Then, $T\overset{d}{\rightarrow}Z^{T}B^{\star}Z,$ where the iid elements of $Z\in\mathbb{R}^{m}$ are distributed $Z_{i}\sqrt{N(\sigma^{2}+\gamma/2)^{d/2}}\sim N\left(\triangle\mu_{i}^{\star}/2,\mu_{i}^{\star}/\left(2\sqrt{\pi}\right)^{d}\right).$

Also, the hypothesis that $\mu^{\star}$ satisfies the shape constraint can be consistently tested at a significance level of $\alpha$ by rejecting the null when $N(\sigma^{2}+\gamma/2)^{d/2}\left(W_{\gamma}(\hat{f},\mu)-\right.$ $\left.W_{\gamma}(f_{unc},\mu)\right)\geq c_{\alpha},$ where $c_{\alpha}$ satisfies $P(X^{T}BX\geq c_{\alpha})=\alpha$ and the iid elements of $X\in\mathbb{R}^{m}$ are distributed $X_{i}\sqrt{N(\sigma^{2}+\gamma/2)^{d/2}}\sim N\left(0,\mu_{i}/\left(2\sqrt{\pi}\right)^{d}\right).$

Proof:

The gradient and Hessian of $W_{\gamma}(f,\mu)$ with respect to $f$ are $\triangledown_{f}W_{\gamma}(f,\mu)=x_{\mu,f}$ and $\triangledown_{f,f}W_{\gamma}(f,\mu)=\triangledown_{f}x_{\mu,f}=\gamma(D_{f}-\psi_{\mu,f}D_{\mathbf{1}\oslash\mu}\psi_{\mu,f}^{T})^{+},$ which are derived in the proof of the next theorem. Thus, the second order Taylor series expansion about $\mu^{\star}=f_{unc}$ is given by,

$W_{\gamma}(\mu^{\star},\mu)=W_{\gamma}(f_{unc},\mu)+x_{\mu,f_{unc}}^{T}(\mu^{\star}-f_{unc})$

Since $f_{unc}:=\text{arg min}_{f}W_{\gamma}(f,\mu)$ and $\triangledown_{f}W_{\gamma}(f,\mu)=x,$ we have $x=\mathbf{0}.$ Given the general discretized densities $\mu_{0},\mu_{1}\in\mathbb{R}^{m},$ the next theorem will also shows that the matrix $D_{\mu_{0}}-\psi D_{\mathbf{1}\oslash\mu_{1}}\psi^{T}$ has one eigenvalue that is zero and $m-1$ eigenvalues that are strictly positive. Note that the pseudo inverse is a continuous function when its domain is restricted to the set of matrices with the same rank (Stewart, 1969). Since $\mu$ and $f_{unc}$ converge in probability to $\mu^{\star},$ the Slutsky theorem implies $B\overset{p}{\rightarrow}B^{\star}.$ Lemma 1 implies that $N(\sigma^{2}+\gamma/2)^{d/2}\left\Vert f_{unc}-\mu^{\star}\right\Vert ^{3}$ converges to zero in probability at a rate of $O_{p}\left(N^{-1/2}(\sigma^{2}+\gamma/2)^{-d/4}\right).$ Combining these results with the limiting distribution of $f_{unc}$ given in Lemma 1 implies that the Taylor series expansion given above can be written as,

Since $\mu$ is a consistent estimator for $\mu^{\star},$ we also have $X\overset{d}{\rightarrow}Z,$ so we also have $T\overset{d}{\rightarrow}X^{T}B^{\star}X.$

To show the test is consistent, suppose $\mu^{\star}(x)$ is not $\rho-$concave. Recall that $W_{\gamma}(\hat{f},\mu)$ converges to the metric $W_{0}(\hat{f}(x),\mu(x))$ asymptotically (Benamou et al, 2015) and that Wasserstein distance metrizes weak convergence of distributions (and convergence in the first two moments). Continuity of the functions $\hat{f}(x),\mu(x)$ and $\mu^{\star}(x)$ implies that the distributions corresponding to these densities converge weakly if and only if the densities themselves converge to one another in probability. Since $\hat{f}$ is in the feasible set and $\mu^{\star}$ is not, $\hat{f}(x)$ cannot converge to $\mu^{\star}(x)$ in probability, so $W_{0}(\mu^{\star},\hat{f})$ does not converge to zero. Also, asymptotically we have $W_{0}(\hat{f},\mu)-W_{0}(\mu^{\star},\mu)>W_{0}(\mu^{\star},\hat{f})$ by the triangle inequality, so $N(\sigma^{2}+\gamma/2)^{d/2}\left(W_{\gamma}(\hat{f},\mu)-W_{\gamma}(f_{unc},\mu)\right)$ diverges under the alternative hypothesis.

The Optimization Algorithm

In this section we will derive a trust region algorithm to find the global minimum of (13) and (14). To do this, we will require the gradient and the Hessian of $W_{\gamma}(\mu,g^{1/\rho}).$ The following Theorem provides these values and shows that the optimization problem is convex for cases in which $\rho\neq0.$ Appendix A contains these derivations for the $\log-$concave case, which corresponds to the limit as $\rho\rightarrow0.$ For notational convenience, the Hessian given below corresponds to the case in which the index $k$ is set equal to $m;$ although this is not a requirement of the theorem.

$\vphantom{c}$

Theorem 6: The gradient of the $W_{\gamma}(\mu,g^{1/\rho})$ is

}

and the Hessian is

}

where $A:=\left

\right],$ $B:=\gamma(D_{g^{1/\rho}}-\psi D_{\mathbf{1}\oslash\mu}\psi^{T})^{+},$ and $C:=\frac{1-\rho}{\rho^{2}}$ $

D_{x_{-k}-x_{k}}.$ In addition, $W_{\gamma}(f,\mu)$ is strictly convex in $f$ when $k$ is chosen to be $arg min_{i}\;x_{i}$ and $\gamma>0,$ and the optimization problem given in (13)-(14) is convex.}

Proof: Since $W_{\gamma}(\mu,g^{1/\rho})$ is differentiable, the envelope theorem implies that the gradient of the objective function in (12) is equal to the gradient of the function given in (13), so $\triangledown_{g}W_{\gamma}(\mu,g^{1/\rho})=D_{g^{1/\rho-1}/\rho}(x_{-k}-x_{k}).$

The derivative of (26) yields the sum of two matrices. Specifically,

Next we will begin by deriving the first term, which will require several intermediate derivatives.

First, since, $f=(g^{1/\rho},m-\mathbf{1}\cdot g^{1/\rho}),$ we have

Second, we will view $w$ as a function of $f$ in order to find $\triangledown_{f}w(f).$ (9) and (10) can be combined to yield the equality $f-D_{w}K\left(\mu\oslash(Kw)\right)=0,$ and implicit differentiation of this equality implies,

Third, the definition $w:=\exp\left(x/\gamma\right)$ implies $x=\gamma\log(w),$ so

Lastly, let $\tilde{x}(x):=x_{-k}-x_{k},$ so $\triangledown_{x}\tilde{x}(x)=\left[

\right]_{m-1\times1}.$ After combining all four derivatives, we have

Since $\triangledown_{g}g^{1/\rho-1}/\rho=(1-\rho)/\rho^{2}D_{g^{1/\rho-2}},$ the second matrix in (28) is given by $C.$

Convexity requires that this Hessian is positive semidefinite. If $k$ is chosen to be $\text{arg min}_{i}\;x_{i},$ then $x_{-k}-x_{k}\geq0.$ Since this is the case, $C$ is a diagonal matrix with nonnegative diagonal elements, so $C$ is positive semidefinite.

Next we will establish that $ABA^{T}$ is positive definite in several steps. First, note that $ABA^{T}$ is symmetric if $B$ is symmetric. Also, since $D_{f}$ and $\psi D_{\mathbf{1}\oslash\mu}\psi^{T}$ are symmetric $B^{+}$ is also symmetric. Since the Moore-Penrose pseudo inverse preserves symmetry, $B$ is also symmetric. We will proceed by establishing a few intermediary results on the components of $ABA^{T}.$

Since $D_{\mathbf{1}\oslash\mu}$ is positive semidefinite, $\psi D_{\mathbf{1}\oslash\mu}\psi^{T}$ is as well. Since $\psi$ is a coupling of the densities $\mu$ and $f,$ we have

In other words, $\mathbf{1}$ is an eigenvector of $D_{\mathbf{1}\oslash f}\psi D_{\mathbf{1}\oslash\mu}\psi^{T}$ with a corresponding eigenvalue of 1. The Perron-Frobenius theorem states that an $m\times m$ matrix with all positive elements and columns that sum to one has a unique eigenvalue that is equal to one and $m-1$ eigenvalues that are strictly less than one. Note that each eigenvalue, say $\lambda,$ and corresponding eigenvector, say $p,$ of $D_{\mathbf{1}\oslash f}\psi D_{\mathbf{1}\oslash\mu}\psi^{T}D_{\mathbf{1}\oslash f}$ satisfies,

so, if $p$ is an eigenvector of $D_{\mathbf{1}\oslash f}\psi D_{\mathbf{1}\oslash\mu}\psi^{T}D_{\mathbf{1}\oslash f},$ then $p$ is an eigenvector of $I-D_{\mathbf{1}\oslash f}\psi D_{\mathbf{1}\oslash\mu}\psi^{T}D_{\mathbf{1}\oslash f},$ with an eigenvalue corresponding to $\tilde{\lambda}:=1-\lambda.$ This implies that $I-D_{\mathbf{1}\oslash f}\psi D_{\mathbf{1}\oslash\mu}\psi^{T}$ is a positive semidefinite matrix with rank $m-1.$ Since multiplication by a positive definite matrix and applying the pseudoinverse preserve both the rank and the signs of the eigenvalues, $B=\gamma\left(D_{f}\left(I-D_{\mathbf{1}\oslash f}\psi D_{\mathbf{1}\oslash\mu}\psi^{T}\right)\right)^{+}$ is also a positive semidefinite matrix with rank $m-1.$

Observation 7.1.8 in Horn and Johnson (1990) implies that the nullspace of $ABA^{T}$ is the same as the nullspace of $BA^{T}.$ Since the eigenvector of $B$ that corresponds to the eigenvalue of zero is $\mathbf{1},$ $ABA^{T}$ is positive definite if there does not exist $v\in\mathbb{R}^{m-1}$ such that $A^{T}v=\mathbf{1}_{m}.$ This system of equations is equivalent to $g_{i}^{1/\rho-1}v_{i}=\rho$ for all $i\in\{1,...,m-1\}$ and $\sum_{i}g_{i}^{1/\rho-1}v_{i}=-\rho,$ which does not have a solution, so $ABA^{T}$ is positive definite.

This, along with the fact that the constraints in (14) are equivalent to constraining $sgn(\rho)g$ to be in the set of concave functions, which is a convex cone, implies that the optimization problem is convex.

Remark: Choosing $k$ to be $\text{arg min}_{i}\;x_{i}$ is a sufficient, but not necessary, condition to guarantee convexity. Choosing $k$ to correspond to the element on the boundary of the mesh over $\mathcal{A},$ with the lowest corresponding value of $w_{i},$ ensures that the density estimate satisfies the shape constraints everywhere on the interior of its domain and rarely results in the objective function being nonconvex along the convergence path.

As described below, we initialize our algorithm at a very good approximation of $\hat{f}.$ Note that this method does not require the specification of $k.$ When the mesh is enlarged beyond the convex hull of the data, so that $\mu_{i}$ and $\hat{f}_{i}$ are lowest when $\mathbf{a}_{i}$ corresponds to a point on the boundary of the domain, $w$ and $v$ generally obtain their minimum value on the boundary points. Thus, initializing the algorithm at a density that is near $\hat{f}$ almost always results in $w$ falling to a sufficiently low value at points corresponding to the boundary of $\mathcal{A}$ to ensure strict positive definiteness of $H$ along the entire convergence path.\footnote{It may be possible to modify the algorithm we use to compute the approximation, as described in Appendix B, to provide a convergence guarantee. This would allow for the formulation of the estimator using $g\in\mathbb{R}^{m}$ rather than $g\in\mathbb{R}^{m-1}.$ This is an ongoing area of research. }

Having already derived the gradient and Hessian of $W_{\gamma}(\mu,g^{1/\rho}),$ it is straightforward to create a trust region algorithm. The algorithm takes an initial density estimate, $f^{(0)},$ as input and in iteration $i$ the algorithm solves

If the value of the objective function evaluated at $g_{-k}^{(i-1)}+\Delta$ results in an improvement over its value at $g_{-k}^{(i-1)},$ then $g_{-k}^{(i)}$ is defined to be $g_{-k}^{(i-1)}+\Delta.$ If the improvement was significant, then the radius of the trust region, $c,$ is increased, and otherwise it is decreased and the value of $g_{-k}^{(i)}$ is defined to be $g_{-k}^{(i-1)}.$ This is described in Algorithm 2.

Figures 1 and 2 illustrate two examples of the output of Algorithm 2. Figure 1 provides density estimates of the rotational velocity of stars that are constrained to be $-2-$concave and $-1/2-$concave, respectively. Figure 2 provides a plot of a two dimensional density; to illustrate the tail behavior of the density more clearly, the logarithm of the density is shown. This estimator uses a dataset containing the height and left middle finger length of 3,000 British criminals that was analyzed by Macdonell (1902) and Student (1908).

The data used to generate both figures are also used by Koenker and Mizera (2010) to compare $\log-$concave density estimates with $-1/2-$concave density estimates. In the case of the dataset for the rotational velocity of stars, they show that the former provides a monotonic density in the region in which the speed of rotation is strictly positive, while the latter density has a peak near the mode of the kernel density estimator shown in Figure 1. This peak is also present in both of the shape-constrained densities shown in Figure 1.

For the dataset used in Figure 2, Koenker and Mizera (2010) show that the logarithm of the maximum likelihood density estimator subject to a $\log-$concavity constraint is below $-24$ near the observation at the very top of Figure 2, so observations this far from the rest of the data would be fairly unlikely to occur if the density was in fact $\log-$concave. The logarithm of the $-1/2-$concave density given below is approximately $-7.2$ near this observation.

To compute these estimates as well as the estimates given in the application section, each iteration of Algorithm 2 used MOSEK, a highly optimized quadratic program solver, to find $\Delta$ efficiently; however, this step is the limiting variable in terms of the time complexity of the algorithm. The computational efficiency of the algorithm can be improved by initializing $f^{(0)}$ at a good approximation of the final density estimate. Also, using the best available approximation of $\hat{f}$ to define $f^{(0)}$ ensures that $w$ is as close as possible to $\mathbf{1},$ for the reasons given in Remark 1 of Theorem 2. This results in an increase in the numerical accuracy of the gradient and Hessian in Algorithm 2, and thus increased stability of the algorithm. The appendix provides a method for finding a particularly accurate and computationally efficient approximation. In many cases initializing Algorithm 2 at this input density results in convergence within two to four iterations, bringing the total computation time of our implementation down to approximately 10 seconds in our MATLAB implementation and 35 seconds in our R implementation when $m=200.$\footnote{We have not finished optimizing these implementations. We believe there is an opportunity for an improvement in this computation time by at least an order of magnitude, particularly in the version written in R. }

The next section provides results on a generalization of the optimization problem given in (13)-(14). One of these results provides sufficient conditions for convexity of this more general optimization problem. It is worth noting that a similar sequential quadratic programming algorithm can also be used to find the global optimum in this more general setting when the derivatives corresponding to $H$ and $r$ exist. Alternatively, one could use gradient descent to solve this more general optimization problem, which would only require the user to provide the first derivative of the transformation of $g.$

Shape-Constraints More Generally

Although $\log-$concavity, and $\rho-$concavity more generally, are the most commonly studied shape constraints, the results provided in the previous sections are also applicable to a large class of new shape constraints. Specifically, this section considers solutions to the general optimization problem given by,

where $\alpha_{i}:\mathbb{R}^{m-1}\rightarrow\mathbb{R}_{+}$ and $\beta_{i}:\mathbb{R}^{m-1}\rightarrow\mathbb{R}.$ Analogous to the estimator described in the preceding sections, we will denote the minimizer of (29)-(30) as $\tilde{g}$ and the generalized density estimator as $\tilde{f}:=f(g),$ where $f(g):=[\alpha(g)^{T},m-\mathbf{1}^{T}\alpha(g)]^{T}.$ The following theorem provides sufficient conditions for the results provided in Theorems 2, 4, and 5 to hold in this more general setting. The statements of these theorems were purposefully ambiguous in regards to the constraint, so we simply provide additional requirements on the functions $\alpha(\cdot)$ and $\beta(\cdot)$ for these generalizations. We also provide sufficient conditions for Theorem 6, but this requires restating the result in its entirety using the new notation. Lastly, point (4) of the theorem provides a new result for strict convexity of $W_{\gamma}(\mu,\alpha(g))$ in the neighborhood of its global minimum, without requiring the assumption that $\alpha(\cdot)$ is convex or concave.

$\vphantom{c}$

Theorem 7: Let $i(x)$ be defined as $\text{arg min}_{i}\left\Vert x-\mathbf{a}_{i}\right\Vert ,$ $\Theta$ as the set of proper and uniformly continuous density functions, $\Omega_{N}$ as $\{g\mid\cap_{i}\beta_{i}(g)\leq0\land\alpha(g)\in\mathbb{R}^{m}\land\sum_{i}\alpha_{i}(g)\leq m\},$ the set $\Lambda$ so that $f(x)\in\Lambda$ if and only if there exists a sequence $\{g^{(N)}\}_{N}^{\infty}$ such that $g^{(N)}\in\Omega_{N}$ for all $N$ sufficiently large and $f(x)=\lim_{N\rightarrow\infty}[\alpha(g^{(N)})^{T},m-\mathbf{1}^{T}\alpha(g^{(N)})]_{i(x)}.$

If $\mu^{\star}(x)\in\Theta,$ both of the sets $\Theta\cap\Lambda$ and $\Omega_{N}$ are nonempty, the function $\alpha(\cdot)$ is invertible, then the following additional conditions are sufficient for the applicability of Theorems 2, 4, 5, and 6 to $\tilde{f}.$

(1) If the codomain of the function $g\mapsto[\alpha(g)^{T},m-\mathbf{1}^{T}\alpha(g)]$ contains $f_{unc}$ then Theorems 2 and 5 hold for $\tilde{f}.$

(2) If each of the functions $\beta_{i}(g),$ as well as $\alpha(g),$ are differentiable in the neighborhood of $\tilde{g}$ asymptotically, $\tilde{g}$ converges to a point on the interior of the domains of these functions, and $\triangledown_{g}\beta_{i}(g)\mid_{g=\tilde{g}}\neq0$ for each $i$ and $\triangledown_{g}\alpha(g)\mid_{g=\tilde{g}}\neq0,$ then Theorem 4 holds for $\tilde{f}.$

(3) Suppose $\{g\mid\beta_{i}(g)\leq0\}$ is convex for all $i,$ and $\alpha_{j}(g)$ is convex (respectively, concave) for all $j.$ If $k:=\text{arg min}_{i}\;x_{i}$ ($k:=\text{arg max}_{i}\;x_{i}$), then (29)-(30) is a convex optimization problem. Also, the gradient and Hessian of $W_{\gamma}(\mu,\alpha(g))$ exist almost everywhere, and at these points they are given by

where $A:=\left

\right],$ $B:=\gamma(D_{\alpha(g)}-\psi D_{\mathbf{1}\oslash\mu}\psi^{T})^{+},$ and $C:=\sum_{j}(x_{j}-x_{k})\triangledown_{g}^{2}\alpha_{j}(g).$ }

(4) Let $d:\mathbb{R}^{m-1}\times\mathbb{R}^{m-1}\rightarrow\mathbb{R}_{+}^{1}$ denote an arbitrary distance measure. If $\alpha(g)\in C^{2}$ and $\triangledown_{g}\alpha(g)$ has full rank, then there exists $\delta>0$ such that $W_{\gamma}(\mu,\alpha(g))$ is convex in $g$ when $d([\alpha(g)^{T},m-\mathbf{1}^{T}\alpha(g)],f_{unc,-k})\leq\delta.$

Proof: (1): The proof of theorems 2 and 5 do not use the $\rho-$concavity constraint, so they still hold as long as there exists $g_{unc}$ such that $f_{unc}=f(g_{unc}).$ Since $f_{unc}$ is the global minimum of $f\mapsto W_{\gamma}(f,\mu),$ $\tilde{f}$ will be equal to $f_{unc}$ whenever $g_{unc}$ is feasible.

(2): The proof given for Theorem 4 uses the first order delta method, which requires the conditions given in the theorem. These conditions are also sufficient for the application of the continuous mapping theorem, which was used to show convergence in probability of the estimator.

(3): Since each $\alpha_{j}(g)$ is convex, Aleksandrov's theorem implies that $\triangledown_{g}\alpha_{j}(g)$ and $\triangledown_{g}^{2}\alpha_{j}(g)$ exist almost everywhere. We will begin by deriving the gradient and Hessian at these points. Let $\omega(x,y,g)$ be defined as

so that we can write $W_{\gamma}(\mu,\alpha(g))$ as

Since $\omega(x,y;g)$ is differentiable in all of its arguments, the envelope theorem implies

By the same logic used in the proof of Theorem 4, we can write the Hessian as,

$ABA^{T}$ is positive definite by the same argument used in Theorem 6, which is also expanded on further in the proof of statement (2). When each $\alpha_{j}(g)$ is convex and $\triangledown_{g}^{2}\alpha_{j}(g)$ exists, then $\triangledown_{g}^{2}\alpha_{j}(g)$ is a positive semidefinite matrix. In this case, $k:=\text{arg min}_{i}\;x_{i},$ so each element of $(x_{-k}-x_{k}\mathbf{1})$ is nonnegative. This implies each term in the sum defining $C$ is a positive semidefinite matrix, so $C$ is positive semidefinite. Also, when each $\alpha_{j}(g)$ is concave and $\triangledown_{g}^{2}\alpha_{j}(g)$ exists, then $\triangledown_{g}^{2}\alpha_{j}(g)$ is negative semidefinite. Choosing $k:=\text{arg max}_{i}\;x_{i}$ implies $(x_{j}-x_{k})\triangledown_{g}^{2}\alpha_{j}(g)$ is a positive semidefinite matrix, so $C$ is also positive semidefinite in this case.

Since $\{g\mid\beta_{i}(g)\leq0\}$ is convex for all $i,$ the intersection of these sets is also convex. This, combined with the strict convexity of $W_{\gamma}(\mu,\alpha(g))$ in $g,$ implies convexity of the optimization problem given in (29)-(30).

(4): The argument given in Remark 1 following Theorem 2 implies that each of the elements of $x$ are the same when $f(g)=f_{unc}.$ When this occurs, we have $x_{-k}-x_{k}=0,$ so $C=\mathbf{0}_{m-1\times m-1}$ when $f(g)=f_{unc}.$

The argument in the second to last paragraph of Theorem 6 implies that $ABA^{T}$ is positive definite when $\triangledown_{g}\alpha(g)$ has full rank and there does not exist $v\in\mathbb{R}^{m-1}$ such that $A^{T}v=\mathbf{1}_{m}.$ This system of equations requires $\left(\triangledown_{g}\alpha(g)\right)v=\mathbf{1}_{m-1}$ and $-\mathbf{1}_{m-1}^{T}\triangledown_{g}\alpha(g)v=1.$ However, if $v$ satisfies $\left(\triangledown_{g}\alpha(g)\right)v=\mathbf{1}_{m-1},$ then $-\mathbf{1}_{m-1}^{T}\triangledown_{g}\alpha(g)v=1-m,$ so there is not a solution to this system of equations. This, combined with the fact that $\triangledown_{g}\alpha(g)$ has full rank, implies $ABA^{T}$ is positive definite.

Since $\alpha(g)\in C^{2},$ the eigenvalues of the $\triangledown_{g}^{2}W_{\gamma}(\mu,\alpha(g))$ are continuous in $g.$ Since we have shown that these eigenvalues are strictly positive when $f(g)=f_{unc},$ continuity implies that there exists $\delta>0$ such that all eigenvalues are nonnegative when $d(f(g),f_{unc})\leq\delta.$

\subparagraph*{Remark 1:}

Some of the assumptions given above were made to simplify the exposition rather than necessity. For example, we can replace assumptions regarding differentiability and rank for all $g$ with similar assumptions in the neighborhood of $\tilde{g}.$ The assumption regarding the existence of the inverse of $\alpha(\cdot)$ is worth mentioning in particular. Cases in which either $\alpha(\cdot)$ or $\alpha^{-1}(\cdot)$ cannot be expressed in a closed form appear to be fairly common, but closed form solutions are not a requirement of the theorem since they can be replaced with their numerical counterparts. The application provided in the next section is one example of this case.

Mechanism design appears to be a particularly fruitful source for applications of this generalized shape constrained density estimator. For example, consider a private values auction model in which bidders have valuations that are drawn from the density $f(x).$ Myerson (1981) defines the virtual valuations function as $J_{f}(x):=x-(1-F(x))/f(x),$ and shows that if $J_{f}(x)$ is monotone increasing, then an auction that awards the item to the highest bidder is optimal in the sense of maximizing expected revenue. It is common in mechanism design to make the stronger assumption that the hazard function, defined by $f(x)/(1-F(x)),$ is increasing or the even stronger assumption that $f(x)$ is $\log-$concave. McAfee and McMillan (1987) show that monotonicity of $x-(1-F(x))/f(x)$ is equivalent to convexity of the function $g(x)=1/(1-F(x)),$ which is used in the the next section to formulate a density estimate subject to the constraint that $f(x)$ satisfies Myerson's (1981) regularity condition.\footnote{This example also demonstrates that $\alpha(\cdot)$ and $\beta(\cdot)$ do not need to be unique, since we could constrain the discretized counterparts of $1/(1-F(x))$ to be convex, $\triangledown_{x}1/(1-F(x))$ to be monotonic, or $\triangledown_{x,x}1/(1-F(x))$ to be positive. }

In addition, Myerson and Satterthwaite (1983) show that bilateral bargaining between a buyer and a seller will only result in trade when the virtual valuation of the buyer and the virtual cost of the seller, defined by $x+F(x)/f(x),$ are both increasing functions. Note that we can define $H(x):=1-F(x)$ and $h(x):=H'(x)=-f(x)$ to show that this last condition is equivalent to monotonicity of $x-(1-H(x))/h(x).$ A reformulation of McAfee and McMillan's (1987) condition for monotonicity of $J_{f}(x)$ shows that this is equivalent to convexity of $g(x):=1/(1-H(x)).$ This allows for the formulation of this shape constraint in an analogous manner as the method used to formulate the shape constraint in the next section.

It would also be interesting to explore constraining a density to have an increasing hazard function. Wellner and Laha (2017) show that this is equivalent to constraining $g(x)=-\log(1-F(x))$ to be convex. In all three of the examples given above, guaranteeing that $C$ is positive definite requires the density estimate to satisfy a set of inequalities that do not appear to have an obvious interpretation. Regardless, statement (2) in Theorem 5 implies that $W_{\gamma}(\mu,\alpha(g))$ is still convex as long as $f(g)$ is sufficiently close to $f_{unc}.$ Initializing the density near this unconstrained minimizer and then checking for convexity in each iteration often results in local convexity of $W_{\gamma}(\mu,\alpha(g))$ along the path of convergence.

Since the eigenvalues of the positive definite matrix $ABA^{T}$ are increasing in $\gamma,$ $\sqrt{\gamma/2+\sigma^{2}}$ can also be increased to ensure that $W_{\gamma}(\mu,\alpha(g))$ is convex over a larger set. This has the added benefit of increasing the dispersion of $f_{unc},$ which results in $f_{unc}$ moving closer to the feasible set in the case of most shape constraints, including all the examples discussed so far. In some cases ensuring convexity by increasing $\sqrt{\gamma/2+\sigma^{2}}$ may result in the density being too disperse. If this occurs, it would be best to compare the resulting density estimate with an estimate subject to a stronger constraint, that allows for the formulation of a convex optimization problem, and check which density estimate fits the data more closely. For example, Ewerhart (2013) shows that a sufficient condition for a density to satisfy Myerson's (1981) regularity condition is $\rho-$concavity for $\rho>-1/2,$ and a $\log-$concavity constraint can be used to ensure that the hazard function is monotonic.

Many other examples of constraints that are commonly imposed on densities in economics are given by Ewerhart (2013). Even though the examples cited in this paper all constrain $g(x)$ to be concave or convex, this is by no means a requirement of Theorem 7. For example, one could define a density estimate of the form given by (29)-(30) to estimate densities subject to any of the examples of shape constraints that are given by Ewerhart (2013).

Myerson's (1981) Regularity Condition

The California Department of Transportation (Caltrans) uses first price auctions to allocate construction contracts. In this section we use data on the bids submitted to Caltrans in 1999 and 2000 to explore whether or not this choice of auction format minimizes the costs to the state of California. As discussed in the previous section, if $f(x)$ is the density of private valuations for the bidders, with a distribution function denoted by $F(x),$ and if the bidders are risk neutral, Myerson (1981) shows that auctions that award the item to the person with the highest bid are optimal when the virtual valuations of the bidders is monotonically increasing.

To examine whether this condition is plausible, we need to estimate the valuations (or, in this case, marginal costs) of the construction firms. Guerre, Perrigne, and Vuong (2000) used the fact that the best response function of bidders in a first price sealed bid auction is an increasing function of the bidders's valuations to show that the valuation of bidder $i$ can be estimated by

where $l$ is the numbers of bidders participating in the auction, $b_{i}$ is $i$'s bid, $\hat{f}_{b}(\cdot)$ is a consistent estimate of the density of bids, and $\hat{F}_{b}(\cdot)$ is its corresponding distribution function. To control for the size of each project, we normalize each bid by Caltrans's engineers's estimates of the cost of each project before estimating $\hat{f}_{b}(\cdot)$ and $\hat{F}_{b}(\cdot)$ for each auction size.

Bajari, Houghton and Tadelis (2006) use the same dataset to estimate the costs of each firm. We follow a similar estimation strategy but make some modifications because our focus is on the costs to the state of California. Specifically, we did not subtract transportation costs from the cost estimates or treat bids from small firms differently than larger firms. Each bid consists of a unit cost bid on each item that the contract requires, and the total bid is equal to the dot product of the number of items required and the unit bid of each item. If small modifications are made to the contract after it is awarded, the final payment to the firm is equal to the dot product of the modified quantities and the unit costs in the original bid. Bajari, Houghton and Tadelis (2006) found evidence that firms are able to make accurate forecasts of these final quantities, so we follow their recommendation and replace the first term in (26) with the final amount that is paid to the firm (normalized by the Caltran's engineers' estimate of the project cost). We also exclude all auctions in which these modifications resulted in a change in the payment received by the firm by more than 3%. After excluding these auctions we were left with 1,393 bids. Lastly, Hickman and Hubbard (2015) showed that the accuracy of the valuations estimates can be improved by applying a boundary correction to $\hat{f}_{b}(\cdot),$ which we also employed in our estimation procedure.

After we estimated the valuations for each firm, we estimated $f_{unc}$ by setting $\sqrt{\gamma/2+\sigma^{2}}$ using Scott's (1992) rule of thumb; however, the resulting virtual valuations function was not monotonic. This could be an innocuous idiosyncrasy of the data or it could be evidence that Caltran's choice of auction format is suboptimal.

To investigate which possibility is more plausible, we find the proposed density estimate subject to Myerson's (1981) regularity condition. To define $\alpha(\cdot)$ we solved for $F(x)$ in the equation introduced in the previous section, $g(x)=1/(1-F(x)).$ This derivative is $f(x)=mg'(x)/g(x)^{2},$ and after discretizing, we defined $\alpha_{j}(g)$ as $(g_{j}-g_{j+1})/g_{j}^{2}.$ The convexity of the objective function was maintained along the entire path of convergence, without requiring that we increase $\sqrt{\gamma/2+\sigma^{2}}$ above the recommendation given in the third section. The input density and the output of the algorithm are shown in Figure 3.

We also performed the test described in Theorem 5. We failed to reject the null hypothesis that the objective function, evaluated at $\hat{f},$ was equal to the objective function evaluated at its unconstrained counterpart, $f_{unc},$ with a $p-$value of 0.93. This is similar to the result of the Kolmogorov-Smirnov (1948) test and the Anderson-Darling (1952) test for the null hypothesis that the sample was drawn from the distribution function of $\hat{f};$ these tests also failed to reject the null with $p-$values of 0.98 and 0.94, respectively. In this case the constraints are inactive at all but 24 points in the right tail in a mesh of 300 points.\footnote{Myerson's (1981) regularity condition can also be expressed as $f(x)^{2}+f'(x)(1-F(x))\geq0,$ so it is always satisfied when the density is increasing. In this case, $f_{unc}$ decreases rapidly to the right of the mode, as shown in Figure 3, so it is not in the set of feasible densities. }

Since the density already appears parsimonious, there is little reason to decrease $\sqrt{\gamma/2+\sigma^{2}};$ however, in the interest of comparing these three tests further, we also estimated the density using Scott's (1992) rule of thumb multiplied by 1/2 rather than 2/3. In this case the $p-$value of our test decreased to 0.32, while the $p-$values of other two tests both increased to 0.99. This divergence in $p-$values underlines the difference between these two approaches. Specifically, as we decrease the smoothing, the distribution function converges to the empirical distribution function over the vast majority of the domain, so tests based on comparisons between a distribution and its empirical counterpart are less likely to reject. In contrast, our statistic measures the discrepancy between the global unconstrained minimum of $f\mapsto W_{\gamma}(\mu,f)$ and the set of feasible densities. The test is most reliable when $f_{unc}$ is a reasonable estimate of $\mu^{\star},$ so we do not recommend setting $\sqrt{\gamma/2+\sigma^{2}}$ to a value that under-smooths $f_{unc}$ in this way.

Conclusion

This paper proposes a density estimator that is defined as the density that minimizes a regularized Wasserstein distance from the input kernel density estimator subject to $\rho-$concavity constraints. This framework provides the advantages of convexity and consistency, and it allows for a generalization that is capable of estimating densities subject to a large class of alternative shape constraints. In addition, it allows for a test of the impact of the shape constraints on the fidelity criterion.

The framework presented here can also be extended to allow $\gamma$ and $\sigma$ to take different values at each column of the matrix $K,$ which would be appealing in two situations. When one would like $f^{\star}$ to be as close as possible to $\mu,$ $\gamma$ and $\sigma$ can be decreased below what would have otherwise been possible in regions where the shape-constrained density estimator is closer to $\mu,$ without interfering with convergence of Algorithm 1. Secondly, this would allow for the development of methods that set $\sqrt{\gamma/2+\sigma^{2}}$ using an adaptive approach that is similar to the one described by Sheather and Jones (1991) for kernel density estimators.

Another promising area for future research that has not already been mentioned would be to extend this framework to allow for the estimation of a regression and the density of residuals simultaneously. Dümbgen, Samworth, and Schuhmacher (2011) showed that this does not result in a convex optimization problem in the maximum likelihood setting, so verifying convexity of the objective function in this case is an active area of research. Note that extending the framework presented here to estimating the mode of a data generating process conditional on covariates, or a modal regression, is straightforward. For example, this could be done by imposing a $\rho-$concavity constraint on the conditional density of the dependent variable. Using a relatively low value of $\rho,$ say $-1$ or $-2,$ could be viewed as similar to a quasi-concavity constraint.\footnote{Note that quasi-concavity is equivalent to uni-modality if and only if $d=1.$} Convexity of the optimization problem in this case follows from Theorem 7.