Shape-Enforcing Operators for Point and Interval Estimators

Shape-Enforcing Operators for Generic Point and Interval Estimators of Functions

\editor

Introduction

A common problem in econometrics, statistics, and machine learning is to estimate and make inference on functions that satisfy shape restrictions. These restrictions might arise either from the nature of the function and variables involved or from theoretical reasons. Examples of the first case include distribution functions, which are nondecreasing and range between zero and one, and height growth charts, which are nondecreasing in age. Examples of the second case include demand functions of utility-maximizing individuals, which are nonincreasing in price according to consumer demand theory; production functions of profit-maximizing firms, which are nondecreasing and quasi-concave in input quantities according to production theory and can also be concave in industries that exhibit diminishing returns to scale; bond yield curves, which are monotone and concave in time to maturity; and American and European call option prices, which are concave and monotone in the underlying stock price and increasing in volatility, according to the arbitrage pricing theory.\footnote{Different, but similar shape restrictions apply to put prices, with the American put price being log-concave in the stock price, for example.}

We propose a method to enforce shape restrictions ex post on any initial generic point and interval estimates of functions by applying functional operators. If an operator has reshaping, invariance, order-preserving, and distance-reducing properties, enforcing the shape restriction improves the point estimates and improves the coverage property of the interval estimates. Thus, the shape-enforced point estimates are closer to the target function than the original point estimates under suitable distances, and the shape-enforced interval estimates have greater coverage and shorter length under suitable distances than the original interval estimates. We show that these properties hold for six different operators that enforce the following restrictions: range, convexity, monotonicity, joint convexity and monotonicity, quasi-convexity, and joint quasi-convexity and monotonicity, as well as for combinations of range with all of the above. We impose the range restriction with a natural operator that thresholds the estimates to the desired range. The double Legendre-Fenchel transform enforces convexity by transforming the estimates into their greatest convex minorants. We focus on the monotone rearrangement to enforce monotonicity (though projection on isotone class can also be used in all composition results, as well as convex combinations of isotone projection with rearrangement). We further develop a new operator to enforce quasi-convexity---a shape that has not been well explored in the literature, although it is common in applications. We also show that the compositions of the monotone rearrangement with the double Legendre-Fenchel and the new quasi-convexity operators yield monotone convex and monotone quasi-convex estimates, respectively. In other words, the application of the convex and quasi-convex operators does not affect the monotonicity of the function. We further demonstrate how to modify the operators to deal with concavity, quasi-concavity, their composition with the monotonicity and range operators, and shape restrictions on transformations of the function.

Our method is generic in that it can be applied to any point or interval estimator of the target function. For example, it works in combination with parametric, semi-parametric and nonparametric approaches to model and estimate the target function. It works with modern machine and deep learning methods that are able to exploit either smoothness or structured parsimony (e.g., approximate sparsity) of target functions. Hence our post-processed point and interval estimates automatically inherit the rates of convergence of these estimators and provably improve these properties in finite samples, while also enforcing qualitative shape restrictions brought by scientific reasoning. Moreover, our method applies without modification to any type of function including reduced form statistical objects such as conditional expectation, conditional density, conditional probability and conditional quantile functions, or causal and structural objects such as dose response, production, supply and demand functions identified and estimated using instrumental variable or other methods. The only requirement to obtain consistent point estimators or valid confidence bands is that the source point estimators be consistent or the source confidence bands be valid. There are many existing methods to construct such point estimators and confidence bands under general sampling conditions, including obtained through frequentist, Bayesian or approximate Bayesian methods.\footnote{ Bayesian methods are often used to quantify the uncertainty of complicated methods where the frequentist quantification is intractable, for example, in deep learning problems. Like in the classical approach, one may impose constraints directly during the estimation, though this is often quite cumbersome and is rarely done in practice. The post-processing can be applied to the unconstrained estimates and be justified on pragmatic grounds, ease of computation, or desire to analyze data without restrictions and accept a menu of restrictions ex-post only after validating them.} Under misspecification these requirements may not be satisfied, but the shape-enforcing operators will bring improvements to the point estimators and confidence bands in a sense that we will make precise. To implement our method, we develop algorithms to compute the Legendre-Fenchel transform of multivariate functions and the new quasi-convexity enforcing operator.

We illustrate the theoretical results with two empirical applications to the estimation of a height growth chart for infants in India and a production function for chemical firms in China. In the case of the growth chart, we impose natural monotonicity in the effect of age, together with concavity that is plausible during early childhood. In the case of the production function, we enforce that a firm's output is nondecreasing and quasi-concave in labor and capital inputs according to standard production theory. We also consider imposing concavity in the effect of the inputs. In both applications we use series least squares methods to flexibly estimate the conditional expectation functions of interest, and construct confidence bands using bootstrap. We quantify the size of strict improvements that imposing shape restrictions bring to point and interval estimates in small samples through numerical simulations calibrated to the empirical applications.

\paragraph{Literature Review.} Due to the wide range of applications of shape restrictions, shape-constrained estimation and inference have received a lot of attention in the statistics community. Classical examples include hildreth1954point, ayer1955empirical, brunk1955maximum, vaneeden, grenander1956theory, groeneboom2001estimation, and mammen1991estimating. We refer to barlowstatistical and robertsonorder for classical references on isotonic regression for monotonicity restrictions, and to Koenker:10 for the work on log-concave density estimation. In terms of risk bounds for estimation, please refer to Zhang02, GuntuAnnIso, han2017isotonic, and references therein for recent developments in isotonic regression; and kuosmanen2008representation, seijo2011nonparametric, Guntuboyina:15, and Han:16:convex for convex regression. bellec2018sharp established sharp oracle inequalities for least squares estimators, when only shape restrictions are known to hold. Moreover, hengartner1995finite, dumbgen2003, and anevski2006 considered the construction of confidence bands for univariate functions under monotonicity or convexity restrictions. Please refer to the book groeneboom2014nonparametric and the survey paper guntuboyina2017nonparametric for more comprehensive reviews on estimation and inference under shape constraints.

Most existing works developed constrained methods via maximum likelihood methods for regression or density estimation that impose only shape restrictions and produce constrained estimates without further restrictions. We remark here that these direct approaches deliver advantages over our approach when such target functions are known to satisfy only the qualitative shape constraints. By contrast, our post-processing approach delivers advantages when any generic target function, in addition to satisfying qualitative constraints, satisfies smoothness or other structured parsimony restrictions (e.g., sparsity). Indeed, our method applies to generic problems, and is not tied to statistical parameters such as regression or density estimation. To explain where the advantages arise, we note that the direct isotone univariate regression converges to the true regression function at the $n^{-1/3}$ rate, which is minimax optimal for the parameter space of monotone functions. If the target function is known to lie in the space of smooth functions (H\"older or Sobolev with smoothness $s>1$), the better and optimal rate \textcolor{black}{$n^{-s/(2s+1)}$} can be achieved by an unconstrained estimator (e.g., Stone:80), making the pure isotonic regression suboptimal in this case. To fix the direct isotonic regression in this case, we would need to impose the smoothness constraints in the estimation directly, which ordinarily is not done in practice, let alone theoretically analyzed. (One exception here is Chernozhukov:15:constrained that considered testing shape restrictions in Banach spaces, with the target function being (possibly partially) identified by general conditional moment condition problems, where shape restrictions induce a lattice structure on the space). Smooth cases and other problems, where unconstrained estimators achieve optimal rates, provide the chief motivation for our approach: in such cases, our method automatically inherits the optimal rate and improves the finite sample properties of the estimator through the distance-reducing properties. \textcolor{black}{On the other hand, unlike constrained estimators, unconstrained estimators require delicate choices of tuning parameters to achieve the optimal rate, although adaptive estimation and inference of smooth functions is possible using the method of l91 ls97,gn10b,gn10,cck14.} It is worthwhile noting that cl2018 investigated asymptotic properties of smoothed isotonic estimators, and jw2009 studied the method of rearrangements for obtaining discrete monotone distributions. However, these works only consider very specific classes of shape-constrained estimators, i.e., isotonic and/or discrete estimators.

Another recurrent problem with imposing shape restrictions in estimation is that the derivation of the statistical properties of the constrained estimators is involved and specific to the estimator and shape restriction. As a consequence, there exist very few distributional results, mainly for univariate functions. The results available for the Grenander and isotonic regression estimators show that these estimators exhibit non-standard asymptotics (including relatively slow rates, since smoothness conditions are not exploited); see guntuboyina2017nonparametric for a recent review. \textcolor{black}{Moreover, Horowitz_Lee_2017 and Freyberger2017InferenceUS have recently pointed out the difficulties of developing inference methods from shape-constrained estimators with good uniformity properties with respect to the data generating process. They showed that for shape restrictions defined by inequalities, the distribution of the constrained estimator depends on where the inequalities are binding, which is unknown a priori. Inference based on this distribution therefore becomes sensitive to how close the inequalities are to binding relative to the sample size.} We avoid all of these complications arising from the constrained estimators by enforcing the restrictions ex post and therefore relying on the distribution of the unconstrained estimators (whenever it is available) to construct the confidence bands. Our confidence interval method can also be applied on top of a different constrained estimator to provide potential improvements when the end-point functions of the generated confidence band do not themselves satisfy the restriction. It is worthwhile noting that the idea of ex post confidence bands was mentioned in Section 4.2 in dumbgen2003, which only discussed two cases on the univariate monotone function and univariate convex function. Moreover, the construction of confidence bands for the convex case in dumbgen2003 is quite different from ours (e.g., their lower bound is not necessarily a convex function).

Our paper generally follows the approach introduced in CFG09, which focused on producing improved generic point and interval estimates of monotone functions using the monotone rearrangement. The class of shape enforcing operators covered by our paper is much bigger and much more useful, with analysis being much more challenging, and we view both aspects as a substantial contribution of our paper. Some of the operators that we consider have been analyzed previously in the literature. Dette:08 apply a smoothed rearranged operator to kernel estimators for monotonization purposes and derive pointwise limit theory. CFG10 applied the monotone rearrangement to deal with the quantile crossing problem and BCCF to impose monotonicity in conditional quantile functions estimated using series quantile regression methods and construct monotonized uniform confidence bands. beare2017weak used the double Legendre-Fenchel transform to construct point and interval estimates of univariate concave functions on the non-negative half-line. Other applications of the double Legendre-Fenchel transform include delgado2012distribution, beare2015nonparametric, beare2016empirical. CFG10 and beare2017weak used an alternative approach to make inference on shape-constrained functions. Instead of applying the shape-enforcing operator to a confidence band constructed from an unconstrained estimator, they constructed confidence bands from the estimator after applying the shape-enforcing operator. To do so, they characterized the distribution of the constrained estimator from the distribution of the unconstrained via the delta method, after showing that the shape-enforcing operator is Hadamard or Hadamard directional differentiable. This approach usually yields narrower confidence bands than ours, but it is computationally more involved and requires additional assumptions and non-standard methods. For example, beare2017weak showed that the bootstrap is inconsistent for the distribution of constrained estimators after applying the double Legendre-Fenchel transform when the target function is not strictly concave. Finally, we refer to matzkin1994restrictions, and chetverikov2017econometrics for excellent, insightful up-to-date surveys on the use of shape restrictions in econometrics.

Relative to the literature, we summarize the major contributions of this paper as follows. (1) We introduce an operator to enforce quasi-convexity and deliver improved point and interval estimates of general multivariate quasi-convex functions. Quasi-convexity extends the notion of unimodality to multiple dimensions and generalizes convexity constraints. Despite its importance, the shape restriction of quasi-convexity has not been well studied in the literature and guntuboyina2017nonparametric listed quasi-convexity as an open area in shape-constrained estimation. (2) We extend the use of the Legendre-Fenchel transform to construct improved point and interval estimates of general multivariate convex functions. (3) We show that the composition of the monotone rearrangement with the Legendre-Fenchel transform can be used to construct improved point and interval estimates of monotone convex functions. (4) We show that the composition of the monotone rearrangement with our quasi-convex operator can be used to construct improved point and interval estimates of monotone quasi-convex functions. (The third and fourth contributions proved to be the most challenging and important steps, where the importance stems from shape restrictions often being a composition of monotonicity with convexity or quasi-concavity). (5) We provide a new algorithm to compute the Legendre-Fenchel transform of multivariate functions. (6) We develop an algorithm to compute our quasi-convex operator. The main advantage of our approach is that it works in conjunction with any generic point estimate (e.g., including recent machine and deep learning methods), or any generic interval estimate (that can be a frequentist confidence band or a Bayesian credible region). Because of genericity, it is able to exploit smoothness or other forms of structured parsimony through the use of the appropriate initial estimator. It inherits the rate properties of the initial estimator, while delivering better finite sample properties through distance-reducing inequalities.\\

\paragraph{Outline.} The rest of the paper is organized as follows. Section (ref) introduces the functional shape-enforcing operators and their properties, together with examples of operators that enforce the shape restrictions of interest. Section (ref) discusses the use of shape-enforcing operators to obtain improved point and interval estimates of functions that satisfy shape restrictions. Section (ref) provides algorithms to compute the shape-enforcing estimators. Section (ref) reports the results of two empirical applications and numerical simulations calibrated to the applications. Section (ref) concludes the paper. The proofs of the main results are gathered in the Appendix.

\paragraph{Notation. } For any measurable function $f: \mathcal{X} \to {\Bbb{R}}$ and $p \geq 1$, let $\|f \|_{p} := \left\{\int_{\mathcal{X}} |f(x)|^p dx\right\}^{1/p}$, the $L^p$-norm of $f$, with $\|f \|_{\infty} := \sup_{x \in \mathcal{X}} |f(x)|$, the $L^{\infty}$-norm or sup-norm of $f$. We drop the subscript $p$ for the Euclidean norm, i.e., $\|x\| := \|x\|_2$. For $p\geq 1$, let $\ell^p(\mathcal{X}) := \{f: \mathcal{X} \to {\Bbb{R}} : \|f\|_p < \infty \}$, the class of all measurable functions defined on $\mathcal{X}$ such that the $L^p$-norms of these functions is finite. For $x,x'\in \mathbb{R}^k$, we say $x\geq x'$ if every entry of $x$ is no smaller than the corresponding entry in $x'$. For two functions $f$ and $g$ that map $\mathcal{X} \to {\Bbb{R}}$ we say that $f \leq g$ if $f(x) \leq g(x)$ for all $x \in \mathcal{X}$. We also use $a\vee b:=\max(a,b)$ and $a\wedge b:=\min(a,b)$ for any $a,b \in {\Bbb{R}}$. For two scalar sequences $a_n$ and $b_n$, the notation $a_n \sim b_n$ means that $a_n/b_n \to 1$ as $n \to \infty$. For an operator $\mathbf{O}: \ell^{\infty}_0(\mathcal{X}) \to \ell^{\infty}_0(\mathcal{X})$, we use $\mathbf{{O}^{-}}f:=-\mathbf{O}(-f)$. For two operators $\mathbf{O_1}: \ell^{\infty}_0(\mathcal{X}) \to \ell^{\infty}_0(\mathcal{X})$ and $\mathbf{O_2}: \ell^{\infty}_0(\mathcal{X}) \to \ell^{\infty}_0(\mathcal{X})$, we define $\mathbf{O_{1}O_{2}}$ to be the composition $\mathbf{O_{1}} \circ \mathbf{O_{2}}$.

Functional Shape-Enforcing Operators

Properties of Shape-Enforcing Operators

Assume that the function of interest, $f$, is real-valued with domain $\mathcal{X} \subset {\Bbb{R}}^k$, for some positive integer $k$. Let $\ell^{\infty}(\mathcal{X})$ be the set of bounded measurable functions. Let $\ell^{\infty}_0(\mathcal{X})$ and $\ell^{\infty}_1(\mathcal{X})$ be two subspaces of $\ell^{\infty}(\mathcal{X})$, such that $\ell^{\infty}_1(\mathcal{X}) \subset \ell^{\infty}_0(\mathcal{X})$; and let $\mathbf{O}: \ell^{\infty}_0(\mathcal{X}) \to \ell^{\infty}_0(\mathcal{X})$ be a functional operator. In our case, $\ell^{\infty}_0(\mathcal{X})$ will be the class of unconstrained functions and $\ell^{\infty}_1(\mathcal{X})$ will be the subclass of functions that satisfy some shape restriction. We first introduce three properties that an operator must satisfy to be considered a shape-enforcing estimator.

In addition to these properties, we consider the following “distance contraction" property.

Some particularly interesting cases of constrained classes $\ell^{\infty}_1(\mathcal{X})$ are subsets of functions that satisfy shape restrictions. In this paper, we focus on seven types of shape restrictions: (1) range, (2) convexity, (3) monotonicity, (4) monotone convexity, (5) quasi-convexity, (6) monotone quasi-convexity, and (7) compositions of range with all of the above. Our methods also apply to the restrictions of concavity and quasi-concavity by noting that if $f$ is concave (quasi-concave), then $-f$ is convex (quasi-convex). In the case of monotonicity we focus on the case of monotonically nondecreasing functions. The methods also apply to monotonically nonincreasing functions noting that if $f$ is nondecreasing then $-f$ is nonincreasing.

\textcolor{black}{

}

Range Restrictions

We first consider the subset of range-constrained functions $\ell^{\infty}_{R}(\mathcal{X}):=\{f \in \ell^{\infty}(\mathcal{X}): \underline{f} \leq f(x)\leq \overline f \textrm{ for all } x \in \mathcal{X}\}$ for some constants $\underline{f} \leq \overline{f}$. A natural range-enforcing operator is as follows.

Let $d_p$ be the distance measure induced by the $L^p$-norm, i.e., $d_p(f,g) = \| f - g\|_p$ for any $f,g \in \ell^p(\mathcal{X})$ and $p\geq1$. The following theorem shows that $\mathbf{R}$ is indeed range-enforcing and distance-reducing with respect to $d_p$.

Convexity

Let $\mathcal{X}$ be a convex subset of $\Bbb{R}^k$. Let $\ell^{\infty}_{S}(\mathcal{X}) := \{ f \in \ell^{\infty}(\mathcal{X}) : \liminf_{x' \to x} f(x') \geq f(x) \textrm{ for all } x \in \mathcal{X} \}$, the set of bounded lower semi-continuous functions on $\mathcal{X}$, and $\ell^{\infty}_C(\mathcal{X}) := \{f \in \ell_S^{\infty}(\mathcal{X}): f(\alpha x + (1-\alpha) x') \leq \alpha f(x) + (1-\alpha) f(x') \textrm{ for all } x,x' \in \mathcal{X}, \alpha \in [0,1] \}$, the subset of convex functions on $\mathcal{X}$. We consider the Double Legendre-Fenchel (DLF) transform as a convexity-enforcing operator. To define this operator, we first recall the definition of the Legendre-Fenchel transform (see, e.g., ConvAna:01).

The function $\xi \mapsto f^*(\xi)$ is a closed convex function (see Lemma (ref) in the Appendix) which is also called the convex conjugate of $f$, and the Legendre-Fenchel transform $\mathbf{L}_{\mathcal{X}}$ is also called the conjugate operator. The Legendre-Fenchel transform is a functional operator that maps any function $f$ to a function of its family of tangent planes, which is often referred to as the dual function of $f$.

Lemma (ref) in the Appendix shows that the double Legendre-Fenchel operator maps any lower semi-continuous function $f$ to its greatest convex minorant, i.e., the largest function $g \in \ell^{\infty}_C(\mathcal{X})$ such that $g \leq f$. The left panel of Figure (ref) illustrates this property with a graphical example. We apply the operator $\mathbf{C}$ to the function $f(x) = [10 \exp(3x/2) - \lfloor 10 x \rfloor - 10]/25$ on $\mathcal{X} = [0,1]$, where $\lfloor \cdot \rfloor$ is the floor function. The convexity-enforced function $\mathbf{C} f$ is the greatest convex minorant of $f$, $\mathbf{C} f(x) = 10[f(\overline{x})(x - \underline{x}) + f(\underline{x})( \overline{x} - x)]$ if $\overline{x} \neq \underline{x}$ and $\mathbf{C} f(x) = f(x)$ otherwise, where $\overline{x} = \lceil 10 x \rceil/10$, $\underline{x} = \lfloor 10 x \rfloor/10$, and $\lceil \cdot \rceil$ is the ceiling function.

The following theorem is an immediate consequence of applying known results from convex analysis.

\textcolor{black}{

}

Monotonicity

Let $\ell^{\infty}_M(\mathcal{X}) := \{f \in \ell^{\infty}(\mathcal{X}) : f(x') \leq f(x) \textrm{ for all } x,x' \in \mathcal{X} \textrm{ such that } x' \leq x \}$, the set of bounded nondecreasing, measurable functions on $\mathcal{X}$. We consider the multivariate monotone rearrangement of CFG09 as a monotonicity-enforcing operator.

Proposition 2 of CFG09 showed that the multivariate increasing rearrangement is monotonicity-enforcing and distance-reducing with respect to $d_p$ for any $p \geq 1$. We state this result as a theorem for the purpose of completeness.

\textcolor{black}{

}

Convexity and Monotonicity

Let $\ell^{\infty}_{CM}(\mathcal{X}) := \ell^{\infty}_C(\mathcal{X}) \cap \ell^{\infty}_M(\mathcal{X})$ be the set of bounded convex and nondecreasing functions on $\mathcal{X}$. We consider the composition of the $\mathbf{C}$ and $\mathbf{M}$ operators to enforce both convexity and monotonicity.

In the right panel of Figure (ref), we apply the operators $\mathbf{M}$ and $\mathbf{CM}$ to the function $f(x) = [10 \exp(3x/2) - \lfloor 10 x \rfloor - 10]/25$ on $\mathcal{X} = [0,1]$. The monotonicity-enforced function $\mathbf{M} f$ in dashed line is not convex, whereas the (convexity and monotonicity)-enforced function $\mathbf{CM} f$ is both monotone and convex. Indeed, $\mathbf{CM} f$ is the greatest convex minorant of $\mathbf{M} f$.

Quasi-convexity

\textcolor{black}{Quasi-convexity is a global property of a function, which is weaker than convexity. A convex function must be quasi-convex, but a quasi-convex function is not necessarily convex. Intuitively, a function $f$ defined on a convex domain is quasi-convex if and only if all its level sets are convex. Please also see its definition in Eq. (ref) below. Quasi-convexity is commonly used in economics because it is an ordinal property, preserved by monotone transformations, which represents well economic relationships such as utility and production functions g04,Koenker:10,c17. Moreover, quasi-convex functions have good optimization properties.}

Consider the set of bounded lower semi-continuous quasi-convex functions on $\mathcal{X}$:

We note that $\ell^{\infty}_C(\mathcal{X}) \subset \ell^{\infty}_Q(\mathcal{X})$, and that for any $f \in \ell^{\infty}_Q(\mathcal{X})$, the lower contour sets, $\mathcal{I}_f(y):=\{x \in \mathcal{X} : f(x) \leq y \}$, are convex for all $y \in {\Bbb{R}}$. For any set $\mathcal{Z} \subset {\Bbb{R}}^k$, let $\mathrm{conv}(\mathcal{Z})$ denote the convex hull of $\mathcal{Z}$. We consider the following new operator to impose quasi-convexity:

The operator $\mathbf{Q}$ transforms any bounded lower semi-continuous function into a quasi-convex function. To see this, recall that a function is quasi-convex if its domain and all its lower contour sets are convex. By construction, $x \in \mathrm{conv}[\mathcal{I}_f(y)] $ if and only if $\mathbf{Q} f(x) \leq y$. Therefore, the lower contour set of $\mathbf{Q} f$ at any level $y \in {\Bbb{R}}$ is $\mathcal{I}_{\mathbf{Q} f}(y) = \{x \in \mathcal{X} : \mathbf{Q} f(x)\leq y\} = \mathrm{conv}[\mathcal{I}_f(y)]$, which is a convex set.

The left panel of Figure (ref) shows a graphical example. We apply the operator $\mathbf{Q}$ to the function $f(x) = [10 \exp(3x/2) - \lfloor 10 x \rfloor - 10]/25$ on $\mathcal{X} = [0,1]$. Here we can see that the function $\mathbf{Q} f$ is the greatest quasi-convex minorant of $f$, $\mathbf{Q} f(x) = \min\{f(x), f(\lceil 10 x \rceil /10)\} $.

\textcolor{black}{

}

Quasi-Convexity and Monotonicity

Let $\ell^{\infty}_{QM}(\mathcal{X}) := \ell^{\infty}_Q(\mathcal{X}) \cap \ell^{\infty}_M(\mathcal{X})$ be the set of bounded quasi-convex and partially nondecreasing functions on $\mathcal{X}$. This case is only relevant when $k>1$ because univariate monotone functions are quasi-convex. We consider the composition of the $\mathbf{Q}$ and $\mathbf{M}$ operators to impose both quasi-convexity and monotonicity.

The comments and example in Remark (ref) also apply to the $\mathbf{QM}$-operator. Thus, the assumption that $\mathcal{X}$ is a rectangular set is sufficient to guarantee that the $\mathbf{Q}$-operator preserves monotonicity.

Range and Other Shape Restrictions

The following theorem shows that the operator $\mathbf{R}$ can be composed with $\mathbf{C}$, $\mathbf{M}$ and $\mathbf{Q}$ to produce range-constrained convex, monotone or quasi-convex functions. Let $\ell^{\infty}_{CR}(\mathcal{X}) := \ell^{\infty}_{C}(\mathcal{X}) \cap \ell^{\infty}_{R}(\mathcal{X})$, $\ell^{\infty}_{MR}(\mathcal{X}) := \ell^{\infty}_{M}(\mathcal{X}) \cap \ell^{\infty}_{R}(\mathcal{X})$, $\ell^{\infty}_{QR}(\mathcal{X}) := \ell^{\infty}_{Q}(\mathcal{X}) \cap \ell^{\infty}_{R}(\mathcal{X})$, $\mathbf{CR} := \mathbf{C} \circ \mathbf{R}$, $\mathbf{MR} := \mathbf{M} \circ \mathbf{R}$, and $\mathbf{QR} := \mathbf{Q} \circ \mathbf{R}$.

The operator $\mathbf{MR}$ can be composed with $\mathbf{C}$ and $\mathbf{Q}$ to produce range-constrained monotone convex or quasi-convex functions. The properties of the resulting operators $\mathbf{CMR} := \mathbf{C} \circ \mathbf{MR}$ and $\mathbf{QMR} := \mathbf{Q} \circ \mathbf{MR}$ follow from combining Theorem (ref) with Theorems (ref) and (ref), respectively. Let $\ell^{\infty}_{CMR}(\mathcal{X}) := \ell^{\infty}_{CM}(\mathcal{X}) \cap \ell^{\infty}_{R}(\mathcal{X})$ and $\ell^{\infty}_{QMR}(\mathcal{X}) := \ell^{\infty}_{QM}(\mathcal{X}) \cap \ell^{\infty}_{R}(\mathcal{X})$.

In the right panel of Figure (ref), we apply the operators $\mathbf{MR}$ and $\mathbf{CMR}$ to the function $f(x) = [10 \exp(3x/2) - \lfloor 10 x \rfloor - 10]/25$ on $\mathcal{X} = [0,1]$. We enforce that the range be in the interval $[0.1,0.9]$. The (monotonicity and range)-enforced function $\mathbf{MR} f$ in dashed line satisfies the monotonicity and range restrictions but is not convex. The (convexity, monotonicity and range)-enforced function $\mathbf{CMR} f$ satisfies the three shape restrictions.

Shape Restrictions on Transformations

The shape operators can be combined with other functions to enforce shape restrictions on transformations of the function $f$. An example is log-concavity where we assume that $\log f$ is concave.\footnote{bagnoli2005log discussed applications of log-concavity to economics and statistics, and analyzed the log-concavity properties of common distributions. } Let $h$ be a real-valued bijection with inverse function $h^{-1}$. We consider the operator $\mathbf{O}_{h}$ that applies the operator $\mathbf{O}$ to the transformation $h(f)$ and then recovers the shape-constrained version of $f$ by inversion, that is \[ \mathbf{O}_{h}f=h^{-1}\circ\mathbf{O}(h\circ f). \] For example, if $f$ is log-concave, then $h(x) = \log x$ and $$ \mathbf{O}_{h}f= \mathbf{C}^{-}_{\log}f = \exp[\mathbf{C}^{-}(\log f)]. $$

The following theorem gives conditions under which the transformations $h$ and $h^{-1}$ preserve the properties of the operator $\mathbf{O}$. Define $\ell_{h,j}^{\infty}(\mathcal{X})=\{ f\in\ell_{0}^{\infty}(\mathcal{X}) : h \circ f\in\ell_{j}^{\infty}(\mathcal{X})\} $ for $j \in \{0,1\}$.

Other Ways of Generating Shape Enforcing Operators

When $\ell_0^\infty(\mathcal{X})$ is a Hilbert space \textcolor{black}{and $\ell_1^\infty(\mathcal{X})$ is a closed set in the $L^2$ metric}, it is possible to construct generic shape-enforcing operators via $L^2$-projection in $\ell_1^\infty(\mathcal{X})$:

The $\mathbf{P}$-operator involves an infinite dimensional optimization program that can be computationally challenging except for special cases. For example, when $\ell_1^{\infty}(\mathcal{X}) = \ell_M^{\infty}(\mathcal{X})$, the $\mathbf{P}$-operator corresponds to the isotonization operator that can be computed using the pool adjacent violators algorithm described in barlowstatistical. The following result, discussed on p.45 of chetverikov2017econometrics and stated here as a theorem for the purpose of completeness, shows that $\mathbf{P}$ is range-enforcing and distance-reducing with respect to $d_2$ under some conditions on $\ell_1^\infty(\mathcal{X})$.

\textcolor{black}{ The condition that $\ell_0^\infty(\mathcal{X})$ is a Hilbert space is satisfied when $\mathcal{X}$ is bounded. The sets $\overline{\ell}_C^\infty(\mathcal{X})$, $\overline{\ell}_M^\infty(\mathcal{X})$, $\overline{\ell}_Q^\infty(\mathcal{X})$ and their intersections are convex and closed.\footnote{\textcolor{black}{The set $\overline{\ell}_O^\infty(\mathcal{X})$ denotes the intersection of $\ell_O^\infty(\mathcal{X})$ with the set of uniformly bounded functions, for $O \in \{C,M,Q\}$. The intersection ensures that $\overline{\ell}_O^\infty(\mathcal{X})$ is closed in the $L^2$ metric. Alternatively, we can ensure that $\ell_O^\infty(\mathcal{X})$ is closed in the $L^2$ metric by restricting $\mathcal{X}$ to be a countable set.}} The condition on the maximum and minimum is a more restrictive Hilbert lattice property. It is satisfied by $\overline{\ell}_M^\infty(\mathcal{X})$, but not by $\overline{\ell}_C^\infty(\mathcal{X})$ and $\overline{\ell}_Q^\infty(\mathcal{X})$. Theorem (ref) therefore covers the isotonization operator, but not the convex and quasi-convex projections. }

Improved Point and Interval Estimation

We show how to use shape-enforcing operators to improve point and interval estimators of a shape-constrained function. Let $f_0 : \mathcal{X} \to {\Bbb{R}}$ be the target function, which is known to satisfy a shape restriction, i.e., $f_0 \in \ell_1^{\infty}(\mathcal{X})$. Assume we have a point estimator $f$ of $f_0$, and an interval estimator or uniform confidence band $[f_l,f_u]$ for $f_0$. These estimators are unconstrained and therefore do not necessarily satisfy the shape restrictions, i.e., $f,f_l,f_u \in \ell^{\infty}_0(\mathcal{X})$ but $f,f_l,f_u \not\in \ell^{\infty}_1(\mathcal{X})$ in general.

There are many different ways to obtain these initial estimators, ranging from parametric to modern adaptive nonparametric methods \textcolor{black}{fan96,li07,hastie09. These methods can be tailored to properties of the target function such as smoothness or sparsity.} A common frequentist confidence band for the function $f_0$ is constructed as $$ f_l(x) = f(x) - c_p s(x), \ \ f_u(x) = f(x) + c_p s(x), $$ where $s(x)$ is the standard error of $f(x)$ and $c_p$ is a critical value chosen such that $$ {\mathrm{P}}(f_0 \in [f_l,f_u]) \geq p, $$ for some confidence level $p$, where event $f_0 \in [f_l,f_u] $ means $\{f_0(x) \in [f_l(x),f_u(x)] \text{ for all } x \in \mathcal{X}\} $. Wasserman:2006 provides an excellent overview of methods for constructing the critical value; \textcolor{black}{see also gn10 and cck14 for constructions of adaptive confidence bands in low-dimensional smooth nonparametric models and bkks20 for a recent proposal in high-dimensional generalized additive models.} With a slight abuse of notation, an initial Bayesian credible region $[f_l,f_u]$ can be constructed similarly with the constant $c_p$ determined such that $$ \Pi \{f_0 \in [f_l,f_u] \mid S\} \geq p, $$ where $S$ denotes data (can be a set of statistics derived from data in robust Bayes procedures, for example, means or empirical moment functions), $[f_l,f_u]$ is a measurable function of $S$, and $\Pi (\cdot \mid S) $ denotes posterior distribution of parameter $f_0$ (viewed as a random element in the Bayesian approach), induced by $S$ and a prior distribution over potential values $f_0$ can take. We give empirical and numerical examples in Section 5.

To enforce the shape restriction, we apply a suitable shape-enforcing operator to the original point estimator and end-point functions of the confidence band. The resulting estimator, $\mathbf{O}f$, and confidence band, $[\mathbf{O} f_l, \mathbf{O} f_u]$, improve over $f$ and $[f_l, f_u]$ in the sense that $f$ lies weakly closer to $f_0$ and the width of the band $[\mathbf{O} f_l, \mathbf{O} f_u]$ is weakly smaller than that of $[f_l, f_u]$, while the coverage is weakly greater. These properties of $\mathbf{O}f$ and $[\mathbf{O} f_l, \mathbf{O} f_u]$ are a corollary of Definition (ref):

Part (1) shows that $[\mathbf{O}f_{l},\mathbf{O}f_{u}]$ provides a coverage improvement over $[f_{l},f_{u}]$ in that $[\mathbf{O}f_{l},\mathbf{O}f_{u}]$ contains $f_0$ whenever $[f_{l},f_{u}]$ does. Part (2) shows that the shape-enforced point estimator improves over the original estimator in terms of estimation error measured by the $\rho$-distance between the estimator and the target function. Parts (1) and (3) show that the shape-enforced confidence band not only has greater coverage but also is shorter with respect to the $\rho$-distance than the original band. These improvements apply to any sample size. In particular, they imply that enforcing the shape restriction preserves the statistical properties of the point and interval estimators. Thus, the shape-enforced estimator inherits the rate of consistency of the original estimator, and the shape-enforced confidence band has coverage at least $p$ in large samples if the original band has coverage $p$ in large samples. Corollary (ref) can therefore be coupled with Theorems (ref)--(ref) to yield improved inference on a function that satisfies any of the shape restrictions considered in the previous section. It is also worthwhile noting that further quantifying the exact size of improvement depends on $f_0$ and properties of the obtained estimators $f$ and $[f_l, f_u]$.

Implementation Algorithms

We provide implementation algorithms for the different shape-enforcing operators based on a sample or grid of $n$ points $\mathcal{X}_n = \{x_1, \ldots, x_n\}$ with corresponding values of $f$ given by the array $\mathcal{Y}_n = \{y_1, \ldots, y_n\}$ with $y_i = f(x_i)$. Computation of the $\mathbf{R}$-operator is trivial, as it amounts to thresholding the elements of $\mathcal{Y}_n$ to be between $\underline{f}$ and $\overline{f}$, i.e., $$ \mathbf{R} f(x_i) = (\underline{f} \vee y_i) \wedge \overline{f}. $$ When $k=1$, CFG09 showed that the $\mathbf{M}$-operator sorts the elements of $\mathcal{Y}_n$. Thus, assume that $x_{1} \leq x_{2} \leq \ldots \leq x_{n}$ and let $y_{(1)} \leq y_{(2)} \leq \ldots \leq y_{(n)}$ denote the sorted array of $\mathcal{Y}_n$. Then, $$ \mathbf{M} f(x_i) = y_{(i)}. $$ When $k>1$, each $\mathbf{M}_{j}$-operator in Definition (ref) can be computed by applying the same sorting procedure to the dimension $j$ sequentially for each possible value of the other dimensions. We refer to CFG09 for more details on computation. We next develop new algorithms for the $\mathbf{C}$ and $\mathbf{Q}$ operators.

Computation of $\mathbf{C}$-Operator

When $k=1$, we can obtain the greatest convex minorant using the standard method based on the pool adjacent violators algorithm described in barlowstatistical. We provide an algorithm for the case where $k>1$. By Definitions (ref) and (ref), the DFL transform of $f$ is the solution to $$ \mathbf{C} f(x) = \sup_{\xi \in f^*(\mathcal{X})} \inf_{\tilde x \in \mathcal{X}} \{ \xi'(x - \tilde x) + f(\tilde x) \}. $$ This is a saddle point problem that might be difficult to tackle directly. However, when $\mathcal{X}$ is replaced by the finite grid $\mathcal{X}_n$, the problem has a convenient linear programming representation:

This program can be solved using standard linear programming methods. In particular, the computational complexity of the standard interior point method for solving (ref) is $O((k+1)(n+(k+1))^{1.5})$, where $k+1$ is the number of decision variables and $n$ is the number of constraints.

The following algorithm summarizes the computation of the $\mathbf{C}$-Operator.

Computation of $\mathbf{Q}$-Operator

We propose a method to compute the $\mathbf{Q}$ operator based on solving problem (ref) on a finite grid, namely $$ \mathbf{Q} f(x) = \min \left\{ y \in \mathcal{Y}_n : x \in \mathrm{conv}[\mathcal{I}_{f,n}(y)] \right\}, $$ where $\mathcal{I}_{f,n}(y) = \{ x_i : y_i = f(x_i) \leq y, i = 1,2,\ldots,n\}$. We find the solution to the program using the following bisection search algorithm:

The binary search algorithm for the $\mathbf{Q}$-Operator runs in \textcolor{black}{$O(\log(n))$} iterations. The major computational cost within each iteration is the check of whether $x$ is in the convex hull in step (4). This check does not require construction of the actual convex hull, which is computationally expensive especially in high dimensions. Instead, it is sufficient to check the existence of a feasible solution of a linear program.

We further note that each of the above two algorithms can be run in parallel across the grid points, because the output of the algorithm for one grid point does not depend on the output for any other grid point. This parallelizability allows for efficient computation on nontrivial grids.

Numerical Examples

Univariate Case

We consider an empirical application to growth charts and a calibrated simulation where the target function $f_0$ is univariate.

Height Growth Charts for Indian Children

Since their introduction by Quetelet in the 19th century, reference growth charts have become common tools to assess an individual's health status. These charts describe the evolution of individual anthropometric measures, such as height, weight, and body mass index, across different ages. See cole for a classical work on the subject, and koenker:charts for a detailed analysis and additional references. Here we consider the estimation of height growth charts imposing monotonicity and concavity restrictions. These restrictions are plausible, since an individual's height is nondecreasing in age at a nonincreasing growth rate during early childhood; see, e.g., twt66 and the growth child standards of the World Health Organization at \url{https://www.who.int/childgrowth/en/}.

We use the data from fnh14 and koenker2011 on childhood malnutrition in India. These data include a measure of height in centimeters, $Y$, age in months, $X$, and 22 covariates, $Z$, for 37,623 Indian children. All of the children have ages between 0 and 5 years, i.e., $X\in \mathcal{X} = \{0,1,\ldots,59\}$. The covariates $Z$ include the mother's body mass index, the number of months the child was breastfed, and the mother's age (as well as the square of the previous three covariates); the mother's years of education and the father's years of education; indicator variables for the child's sex, whether the child was a single birth or part of a multiple birth, whether the mother was unemployed, whether the mother's residence is urban or rural, and whether the mother has each of: electricity, a radio, a television, a refrigerator, a bicycle, a motorcycle, and a car; and factor variables for birth order of the child, the mother's religion and quintiles of wealth.

We assume a partially linear model for the conditional expectation of $Y$ given $X$ and $Z$, namely $$ {\mathrm{E}}[Y \mid X=x, Z = z] = f_0(x) + z'\gamma. $$ The target function is the conditional average growth chart $x \mapsto f_0(x)$, which we assume to be nondecreasing and concave. Since $X$ is discrete, we can express $f_0(x) = P(x)'\beta$, where $P(x)$ is a vector of indicators for each value in $\mathcal{X}$, i.e., $P(X) = [1(X=0), 1(X=1), \ldots, 1(X=59)]'$. We estimate $\beta$ and $\gamma$ by least squares of $Y$ on $P(X)$ and $Z$, and construct a confidence band for $f_0$ on $\mathcal{X}$ using weighted bootstrap with standard exponential weights and 200 repetitions praestgaard1993exchangeably,hahn1995bootstrapping. Standard errors are estimated using bootstrap rescaled interquartile ranges chernozhukov2013inference, and the critical value is the bootstrap $0.95$-quantile of the maximal $t$-statistic. Weighted bootstrap is computationally convenient in this application because it is less sensitive than empirical bootstrap to singular designs, which are likely to arise in the bootstrap resampling because $Z$ and $P(X)$ contain many indicators.

Figures (ref) and (ref) report the point estimates and 95% confidence bands of $f_0$ for the entire sample and a random extract with $1{,}000$ observations, respectively. We use the subsample to illustrate the deviations from the shape restrictions that are more apparent when the sample size is small. The original estimates are displayed in the left panels, and the estimates imposing monotonicity and concavity in the right panels. The original estimates in the entire sample are nondecreasing in age except at 45 months, and deviate from concavity in some areas. The $\mathbf{M}$ and $\mathbf{{C}^{-}}$ operators correct these deviations. The estimates in the random extract of the data clearly show deviations from both monotonicity and concavity. The $\mathbf{M}$ and $\mathbf{{C}^{-}}$ operators fix these deviations and produce point estimates that are closer to the estimates in the entire sample.

Calibrated Monte Carlo Simulation

We quantify the finite-sample improvement in the point and interval estimates of enforcing shape restrictions using simulations calibrated to the growth chart application. The child's height, $Y$, is generated by $$ Y_i = \mathbf{{C}^{-}M} [P(X_i)'\hat \beta] +Z_i'\hat\gamma+ \hat\sigma \epsilon_i, \ \ i = 1,\ldots,n, $$ where $P(X_i)$ is the vector of indicators for all the values of $\mathcal{X}$; $\hat \beta$, $\hat\gamma$ and $\hat\sigma$ are the least squares estimates of $\beta$, $\gamma$ and the residual standard deviation in the growth chart data; and $\epsilon_i$ are independent draws from the standard normal distribution. The application of the $\mathbf{{C}^{-}M}$-operator guarantees that the target function $f_0(x) = \mathbf{{C}^{-}M} [P(x)'\hat \beta]$ is monotone and concave. We consider six sample sizes, $n \in \{500, 1{,}000, 2{,}000, 4{,}000, 8{,}000, 37{,}623\}$, where $n = 37{,}623$ is the same sample size as in the empirical application. The values of $X_i$ and $Z_i$ are randomly drawn from the data without replacement. The results are based on 500 simulations. In each simulation we construct point and band estimates of $f_0$ using the same methods as in the empirical application.

Table (ref) reports simulation averages of the $d_{\infty}$-distance between the estimates and target function, coverage of the target function by the confidence band and $d_{\infty}$-length of the confidence band for the original and shape-enforced estimators. We consider enforcing concavity with the $\mathbf{{C}^{-}}$-operator, monotonicity with the $\mathbf{M}$-operator, and both concavity and monotonicity with the $\mathbf{{C}^{-}M}$-operator. The improvements from imposing the shape restrictions are decreasing in the sample size, but there are substantial benefits in estimation error even with the largest sample size. Enforcing monotonicity has generally stronger effects than enforcing concavity, but both help improve the estimates. Thus, the $\mathbf{{C}^{-}M}$-operator produces the best point and interval estimators for every sample size. For the smallest sample size, the reduction in estimation error is almost 37% and the improvement in length of the confidence band is more than 20%. The gains in coverage probability are also substantial, especially for the smaller sample sizes. Overall, the simulation results clearly showcase the benefits of enforcing shape restrictions, even with large sample sizes.

\textcolor{black}{ We compare the $d_{\infty}$-error of several estimators in the simplified design $$ Y_i = \mathbf{{C}^{-}M} [P(X_i)'\hat \beta] + \hat\sigma \epsilon_i, \ \ i = 1,\ldots,n, $$ where $P(X_i)$, $\hat \beta$, $\hat\sigma$, $\epsilon_i$ and $n$ are the same as for Table (ref). We consider unconstrained, shape-constrained, shape-enforced, and combinations of shape-enforced and shape-constrained estimators. The unconstrained estimators include the same estimator as in Table (ref) (Piecewise Constant) and a locally linear estimator with data-driven choice of bandwidth (Locally Linear).\footnote{The piecewise constant estimator can be viewed as a locally constant estimator with bandwidth equal to zero. The locally linear estimator is computed using the package KernSmooth KernSmooth19 with the bandwidth chosen by the plug-in method of rsw95.} We consider two classical shape-constrained estimators: the isotonic regression estimator (Isoreg) that imposes monotonicity and the concave regression estimator (Conreg) that imposes concavity.\footnote{We compute the isotonic regression using the R command isoreg r19, and the concave regression using the package cobs cobs20.} We illustrate how to combine shape-enforced operators with shape-constrained estimators by applying the $\mathbf{{C}^{-}}$-operator to the isotonic regression estimator to enforce monotonicity and concavity. Finally, we compare the $\mathbf{{C}^{-}}$-operator with the $\mathbf{S{C}^{-}}$-operator defined in Remark (ref).}

\textcolor{black}{Table (ref) shows the results based on $5{,}000$ simulations. The comparison between shape-constrained and shape-enforced estimators produces mixed results, which vary with the estimator, shape restriction and sample size. Thus, the $\mathbf{M}$-operator outperforms Isoreg for both unconstrained estimators, whereas Conreg outperforms the $\mathbf{C}$-operator applied to Piecewise Constant. The unconstrained locally linear estimator outperforms Isoreg, Conreg and the shape-enforced estimators applied to Piecewise Constant for most sample sizes, despite the target function not being smooth. This finding highlights the benefit of using estimators that exploit smoothness when the sample size is not large. On the other hand, the shape-enforcing operators are more effective when applied to estimators such as Piece Constant and Isoreg that do not rely on smoothness. Shifting the $\mathbf{{C}^{-}}$-operator to deal with potential bias generally reduces estimation error for all the estimators considered.}

Multivariate Case

We consider an empirical application to production functions and a calibrated simulation where the target function $f_0$ is bivariate.

Production Functions of Chinese Firms

The production function is a fundamental relationship in economics that maps the quantity of inputs, such as capital, labor and intermediate goods, to the quantity of output of a firm. When there are only two inputs, the law of diminishing marginal rate of technical substitution dictates that the production function of a firm is nondecreasing and quasi-concave in the inputs Hicks_Allen_1934. If in addition the industry exhibits diminishing returns to scale, then the production function is concave in the inputs. We use the data from Jacho-Chavez_Lewbel_Linton_2010, and Horowitz_Lee_2017 to estimate the production function of Chinese firms in the chemical industry. These data contain information on real value added (output), real fixed assets (capital) and number of employees (labor) for 1,638 firms in 2001.\footnote{Following Jacho-Chavez_Lewbel_Linton_2010 and Horowitz_Lee_2017, we drop observations with a capital-to-labor ratio below the $0.025$ sample quantile or above the $0.975$ sample quantile.} We estimate a production function using these data and enforce the monotonicity and quasi-concavity restrictions. We provide results from enforcing concavity only for illustrative purposes because the chemical industry might exhibit increasing returns to scale at some levels of the inputs.

Figure (ref) shows 3-dimensional estimates and 95% confidence bands for the average production function, together with upper contour sets for the point estimates. The estimates and bands are displayed in a region defined by the tensor product of two grids for labor and capital. Each grid includes 20 equidistant points from the 10% to the 90% sample percentiles of the corresponding variable. We obtain the unconstrained estimates from least squares with the tensor product of third-degree global polynomials as the two marginal bases for capital and labor. The confidence bands are constructed using weighted bootstrap with standard exponential weights and 500 repetitions. Standard errors are estimated using bootstrap rescaled interquartile ranges and the critical value is the bootstrap $0.95$-quantile of the maximal $t$-statistic. Many of the upper contour sets of the unconstrained point estimates are far from being convex, and thus imply a violation of quasi-concavity. In fact, violations of monotonicity occur over a considerable area---most notably, see the positive slopes of the contour curves at low levels of labor and high levels of capital (the upper-left region of the contour plot).

The second row of Figure (ref) shows the results after the $\mathbf{{Q}^{-}M}$-operator is applied to the point estimates and to each end-point function of the confidence band to ensure monotonicity and quasi-concavity. The contour curves are convex by construction, and thus satisfy the quasi-concavity restrictions. Finally, the third row of Figure (ref) shows the results after the $\mathbf{{C}^{-}M}$-operator is applied to enforce monotonicity and concavity. Although quasi-concavification of a production function estimate is always reasonable, whether restriction to concavity is appropriate depends on prior knowledge of the industry.

Calibrated Monte Carlo Simulation

Similar to the univariate case, we now explore the finite-sample improvements from enforcing shape restrictions via simulations calibrated to the production function application. The output, $Y$, of each firm is generated by $$ Y_i = \hat \gamma + \hat \beta_1 L_i + \hat \beta_2 K_i + \hat \sigma \epsilon_i, \ \ i = 1,\ldots,n, $$ where $\hat \gamma$, $\hat \beta_1$, $\hat \beta_2$ and $\hat \sigma$ are calibrated to the least squares estimates and the residual standard deviation of this linear regression model in the production function data; $\epsilon_i$ are independent draws from the standard normal distribution; and $n$ is the sample size of the simulated data. The vector $(L, K)$ of labor and capital is drawn without replacement from the original data. The target function is $$ f_0(l,k) := {\mathrm{E}}[ Y \mid L = l, K =k] = \hat \gamma + \hat \beta_1 l + \hat \beta_2 k, $$ which is increasing and concave in the capital and labor inputs because $\hat \beta_1 > 0$ and $\hat \beta_2 > 0$. We consider five sample sizes, $n \in \{100, 200, 500, 1{,}000, 1{,}638\}$, where $n = 1{,}638$ is the same sample size as in the empirical application. The results are based on 5,000 simulations. In each simulation we construct point and band estimates of $f_0$ using the same methods as in the empirical application.

Table (ref) reports the same diagnostics as Table (ref). The operators $\mathbf{{Q}^{-}M}$ and $\mathbf{{C}^{-}M}$ perform similarly in this case. Both bring substantial gains in estimation and inference\textcolor{black}{, and the shifted variants bring additional gains. Shifting the $\mathbf{{C}^{-}M}$ operator in particular has a notable effect on estimation error for small sample sizes}. The operators reduce estimation error between 5% and \textcolor{black}{31%} and the width of the confidence band between 14% and 37% in the sup-norm, depending on the sample size. The operators also improve the coverage of the confidence bands, especially for the smaller sample sizes. Indeed, enforcing the constraints compensates for the undercoverage of the unconstrained estimates for most of the sample sizes considered.

Conclusion

In this paper, we investigate a pool of shape-enforcing operators, including range, rearrangement, double Legendre-Fenchel, quasi-convexification, composition of rearrangement and double Legendre-Fenchel, and composition of rearrangement and quasi-convexification operators. We show that enforcing the shape restrictions through these operators improves point and interval estimators, and provide computational algorithms to implement these shape-enforcing operators. It would be useful to develop operators to enforce other shape restrictions, such as \textcolor{black}{supermodularity} or the Slutsky conditions for demand functions. We leave this extension to future research.

Acknowledgments

We are very grateful to Simon Lee for kindly sharing the data for the production function application and to Roger Koenker for kindly making the Indian nutrition data available through his website. We thank the editor Garvesh Raskutti, two anonymous referees, Shuowen Chen and Hiroaki Kaido for comments. We gratefully acknowledge research support from the National Science Foundation and the Spanish State Research Agency MDM-2016-0684 under the Mar\'ia de Maeztu Unit of Excellence Program. Xi Chen is supported by NSF IIS-1845444. Part of this work was completed while Fern\'andez-Val was visiting CEMFI and NYU. He is grateful for their hospitality.