The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
64,706 characters
Frequentist Shrinkage Under Inequality Constraints
\maketitle
\begin{abstract}
This paper shows how to shrink extremum estimators towards inequality constraints motivated by economic theory. We propose an Inequality Constrained Shrinkage Estimator (ICSE) which takes the form of a weighted average between the unconstrained and inequality constrained estimators with the data dependent weight. The weight drives both the direction and degree of shrinkage. We use a local asymptotic framework to derive the asymptotic distribution and risk of the ICSE. We provide conditions under which the asymptotic risk of the ICSE is strictly less than that of the unrestricted extremum estimator. The degree of shrinkage cannot be consistently estimated under the local asymptotic framework. To address this issue, we propose a feasible plug-in estimator and investigate its finite sample behavior. We also apply our framework to gasoline demand estimation under the Slutsky restriction.
\end{abstract}
{\bf Keywords:} James-Stein, extremum estimators, nonlinear models, economic restrictions.
\section{Introduction}
Inequality constraints are common in applied economic research. Typical examples are monotonicity constraints on utility or production functions, restrictions on estimated covariance matrices such as positive definiteness, restrictions on the Slutsky matrix, etc. If the imposed constraints hold, we will get more efficient estimates. If not, the estimates will be biased.
The paper proposes an alternative way to use economic theory in estimation. We introduce a generalized shrinkage estimator that shrinks an estimator that ignores theoretical restrictions towards inequality constraints motivated by theory. The inequality constrained shrinkage estimator (ICSE) takes a simple weighted average form between the unconstrained and inequality constrained estimators, with the data-driven weight inversely proportional to the loss function evaluated at the two estimates. We show that the degree of shrinkage depends on which constraints bind, thus, both the direction and degree of shrinkage are fully data driven.
We show that under certain conditions the ICSE outperforms the unrestricted estimator, regardless of what the true data generating process is and whether the theory is correct or not. We demonstrate that the ICSE has a smaller asymptotic risk than the unrestricted estimator uniformly over the parameter space local to the restricted (shrinkage) parameter space. The theory we present applies to a large set of extremum estimators, such as the Generalized Method of Moments (GMM) estimator, the Maximum Likelihood estimator (MLE), the Minimum Distance (MD) estimator, etc.
We use the local asymptotic framework to analyze the performance of the ICSE. To be precise, we assume that the parameter space is located in a $n^{-1/2}$-neighborhood of the restricted space, reflecting the belief that the imposed theoretical restrictions are only ''approximately correct''. In contrast to the generalized James-Stein estimator, the asymptotic distribution of the ICSE is not normal. Since the ICSE is a weighted average of the unconstrained and inequality constrained estimators, the asymptotic distribution of the shrinkage estimator inherits the non-normality of the inequality constrained estimator.
Under the local asymptotic framework, it is impossible to consistently estimate the optimal degree of shrinkage, as it depends on the local $O(n^{-1/2})$ parameters (see e.g. \citealp{hjort_claeskens2003}). To address this issue, we propose a feasible plug-in estimator based on the asymptotically unbiased estimator of the local parameters. However, this makes the estimated shrinkage parameter asymptotically random, which affects the asymptotic distribution of the averaging weight. As a result, the feasible estimator is not consistent and the dominance result may not hold.
In our Monte Carlo study we investigate the finite sample performance of the feasible ICSE along with the generalized James-Stein estimator of \cite{hansen2016}, the Empirical Bayes (EB) estimator, the unrestricted estimator, and the restricted estimator. Simulations show that the feasible ICSE dominates the unrestricted estimator in terms of mean squared error. Moreover, it also dominates the generalized James-Stein estimator in cases when a subset of the constraints bind. We also show that the ICSE performs better than the EB estimator when the constraints are violated or close to bind, while the EB estimator dominates the ICSE when the constraints are satisfied as strict inequalities.
In our application we consider gasoline demand estimation under the Slutsky restriction. We estimate the demand curves across three income groups corresponding to the first, second, and third quartiles, respectively. We show that the shrinkage effect is more prominent for the low income group, since consumers with low income are less likely to have upward sloping demand curves. In a similar application, \cite{kasy2018} use the Empirical Bayes framework to show that the degree of shrinkage is similar across different groups.
\vspace{1em}
The literature on shrinkage estimation begins with \cite{stein1956} who observed that the unconstrained estimator in a Gaussian location model is inadmissible when the dimension of the parameter vector is greater than two. This lead to a seminal paper by \cite{james_stein1961} where they proposed a shrinkage estimator that dominates the MLE. \cite{baranchik1964} showed that the James-Stein estimator is inadmissible and dominated by its positive part version. However, even the positive part James-Stein estimator is inadmissible. \cite{shao1994} propose a piecewise linear estimator that has even smaller risk. Theory for risk analysis of shrinkage estimators was provided by \cite{stein1981}. \cite{hansen2015} compares the performance of different shrinkage estimators and provides corresponding efficiency bounds.
All of the aforementioned estimators shrink the parameters towards zero. In contrast, \cite{oman1982a,oman1982b} introduce estimators which shrink towards linear subspaces. \cite{delNegro2004} show how to shrink to non-linear subspaces in the Bayesian framework, using a DSGE model-based prior to estimate the VAR impulse response functions. In their recent paper, \cite{kasy2018} provide an Empirical Bayes framework which allows to shrink to various theoretical restrictions in form of both equalities and inequalities. Our paper complements the aforementioned literature by extending the Stein's type shrinkage argument to non-linear inequality constraints.
\cite{james_stein1961} first showed that the shrinkage estimator dominates the unrestricted MLE in exact normal sampling. \cite{hansen2016} provides a generalized James-Stein type estimator for parametric models and shows that it dominates the MLE in a pointwise locally asymptotic sense.\footnote{For a given real vector $c$, the pointwise local asymptotic analysis considers a sequence of localized parameters $\theta_{n} = cn^{-1/2}$, and derives the asymptotic (truncated) risk of the averaging estimator under $\theta_{n}$ for given $c$. Such analysis will produce a pointwise risk function for the shrinkage estimator.} \cite{hansen2017} shows that a shrinkage estimator that shrinks the OLS estimator towards the 2SLS estimator has a smaller asymptotic risk than the ordinary OLS estimator. \cite{ditraglia2016} studies the averaging GMM estimator with the averaging weight based on the focused moment selection criterion. The results in the paper suggest that the averaging estimator does not uniformly dominate the conservative estimator. Unlike the aforementioned papers using the pointwise local asymptotic framework, \cite{cheng_et_al2019} establish the uniform dominance result of the GMM averaging estimator over the conservative estimator.
This paper is also closely related to the frequentist model averaging literature. \cite{hansen2007} introduces a model averaging estimator for linear nested models and shows that it is asymptotically optimal. He proposes to minimize a Mallows criterion to select the model weights, which is asymptotically equivalent to minimizing the squared error. \cite{wan_et_al2010} show that the latter result holds not only for discrete but also for continuous model weights and under a non-nested set-up. \cite{hansen_and_racine2012} show that the optimal weights can be obtained by minimizing the cross validation criterion, which allows for a more efficient use of data. Moreover, their approach allows to easily accommodate for heteroskedasticity. \cite{liu2015} points out that the asymptotic distribution of data dependent weights is non-standard, which complicates the inference. He augments the results from \cite{hjort_claeskens2003} and \cite{claeskens_hjort2008} and proposes a procedure that delivers asymptotically correct coverage probabilities for model averaging estimators. \cite{zhu_et_al2017} use a $J$-fold cross-validation criterion to construct optimal averaging weights for model averaging estimators under inequality constraints.
There is a large literature studying estimation under inequality constraints. \cite{andrews1999boundary} derives the asymptotic distribution of extremum estimators when the parameter of interest is on a boundary of the parameter space. His approach solves an asymptotically equivalent problem by minimizing a stochastic quadratic objective function over a convex cone that approximates the parameter space. The approach follows \cite{chernoff1954}, \cite{feder1968}, \cite{pollard1985}, and \cite{wolak1989}. Andrews extends the results in these papers and allows for cases when the estimator objective function is undefined in the neighborhood of the true parameter.
There has been a growing interest in shrinkage estimators in the modern statistics literature. The main idea there is that shrinkage can be introduced through a penalty imposed on the estimator objective function. The most famous example is LASSO (\citealp{tibshirani1996}), which simultaneously shrinks and selects variables. Another seminal example is a Ridge regression, which shrinks the coefficients to zero, but does not perform selection. More complicated penalties lead to more interesting shrinkage spaces, e.g. a fused LASSO (\citealp{friedman_et_al2007}) can be used to shrink time-varying parameters towards random walk, i.e. it penalizes absolute time deviations of the form $|\theta_{t} - \theta_{t-1}|$. Another example is a nearly isotonic regression (\citealp{tibshirani2011nearly}) which shrinks the sequence of points towards a monotone sequence, i.e. it penalizes only positive part deviations $(\theta_{i} - \theta_{i+1})_{+}$.
The remainder of the paper is organized as follows. Section \ref{sec:model} presents the general framework, describes the choice of shrinkage direction and the local asymptotic framework. Section \ref{sec:estimation} introduces the inequality constrained shrinkage estimator. Section \ref{sec:distribution} derives the asymptotic distribution of the estimator. Section \ref{sec:asymptotic_risk} presents the risk dominance result. Section \ref{sec:weight} provides a feasible estimator for the data-dependent weight. Section \ref{sec:mc} demonstrates the finite sample performance of the ICSE in a series of simulations. In Section \ref{sec:slutsky} we apply the method to estimate gasoline demand under the Slutsky restriction. Section \ref{sec:conclusion} concludes. All the mathematical proofs and additional details are left to the Appendix.
\vspace{1em}
We use the following notation throughout the paper: $\mathcal{I}_{n}$ denotes an $n \times n$ identity matrix. $\mathds{1}\{x \geq a\}$ is the indicator function that equals to one if $x \geq a$ and zero otherwise. We use $(x)_{+} = \max\{0,\,x\}$ to denote the ``positive part'' function. Finally, if $x$ is a vector, we use $x > a$ to denote each vector entry being strictly greater than $a$, the same holds for $x < a$.
\section{Model} \label{sec:model}
Suppose we observe a random array $\bm{X}_{n} = \{X_{in}\}_{i=1}^{n}$ of iid realizations. Let $Q_n(\theta)$ denote an extremum estimator objective function that depends on $\bm{X}_{n}$, for example, a GMM criterion or log likelihood function. The objective function is indexed by a parameter $\theta \in \Theta \subset \mathbb{R}^m$.
The goal is to estimate the parameter of interest $\theta$ in a setting augmented by the belief that the true value of $\theta$ may be close (in a sense to be made clear later) to a restricted parameter space $\Theta_{0} \subset \Theta$ defined by a parametric restriction
\begin{equation} \label{eq:restricted_set}
\Theta_{0} = \{ \theta \in \Theta : r(\theta) \geq 0 \},
\end{equation}
where $r(\theta)$ is a differentiable function that maps $\mathbb{R}^{m} \rightarrow \mathbb{R}^{p}$. Let $R(\theta)$ denote the derivative $\frac{\partial}{\partial \theta'} r(\theta)$.
The pivotal point is that the true parameter value $\theta_{0}$ may not satisfy the restrictions, i.e. $\theta_{0}$ does not necessarily lie within $\Theta_{0}$. The restriction can be rather treated as a reasonable belief or ``prior'' about the likely value of $\theta_{0}$. It means that the empirical implication of the imposed theoretical restrictions are only ``approximately correct''.
\begin{remark}
\textnormal{In this paper I focus on the parameter $\theta$ itself, the presented theory can be extended to functions of $\theta$ using the delta-method approach. However, one has to be cautious, since the level of shrinkage depends on the dimension of the function's output.}
\end{remark}
A common example is sign restrictions on all parameters, i.e. $p = m$. In this case the restricted space is $\Theta_{0} = \{ \theta \in \Theta : \theta \geq 0 \}$, where $r(\theta) = \theta$ and $R$ is simply an $m \times m$ identity matrix. The researcher may want to impose sign restrictions only on a subset of parameters. We can easily allow for that by partitioning the parameter space
\begin{equation*}
\theta =
\begin{pmatrix}
\theta_1 \\ \theta_2
\end{pmatrix}
\quad \quad
\begin{matrix}
m - p \\ p
\end{matrix}
\end{equation*}
then the sign restrictions take the form $r(\theta) = \theta_2$, and $R = [0_{p \times (m-p)} \,\vdots\, \mathcal{I}_{p}]$.
In general, $\Theta_{0}$ may be a non-linear subspace. This can be especially useful for structural estimation when an economic model implies non-linear inequality constraints on structural parameters.
\begin{example} \label{exmp:int_rate}
\textnormal{In macroeconomics inequality restrictions often arise in estimation of DSGE models. \cite{moon2009} study an example of interest rate feedback rules, which we briefly describe here. Consider the following interest rate policy rule}
\begin{equation} \label{eq:int_rate_rule}
R_{t} = \rho_{R}R_{t-1} + (1 - \rho_{R})\psi_{1}\pi_{t} + (1 - \rho_{R})\psi_{2}x_{t} + \varepsilon_{R,\,t},
\end{equation}
\textnormal{where $R_{t}$ is the nominal interest rate in period t, $\pi_{t}$ is the inflation rate, and $x_{t}$ is a measure of real activity, such as output deviations from trend or output growth. The shock $\varepsilon_{R,\,t}$ captures unexpected deviations from the systematic component of the policy rule. To address potential endogeneity of both inflation and output in equilibrium, the researcher needs instrumental variables. Lagged variables of inflation and output are natural candidates. According to a large class of DSGE models, output does not fall in a response to an expansionary monetary shock, which leads to a moment restriction $\mathbb{E}[-x_{t}\varepsilon_{R,\,t}] \geq 0$.}
\textnormal{One can estimate the model using the Generalized Method of Moments.\footnote{For the ease of exposition, we skip the details regarding the representation of a typical DSGE model and its solution.} Let $X_{t} = (R_{t-1},\,\pi_{t},\,x_{t})'$ be the vector of regressors, $Z_{t} = (R_{t-1},\,\pi_{t-1},\,x_{t-1})'$ be the vector of IVs, and $\theta = (\rho_{R},\,(1-\rho_{R})\psi_{1},\,(1-\rho_{R})\psi_{2})'$ be the parameter vector. Based on \eqref{eq:int_rate_rule}, one can form a finite sample moment condition $g_{t}(X_{t},\,Z_{t},\,R_{t};\theta) = T^{-1}\sum_{t=1}^{T} Z_{t}(R_{t} - X_{t}'\theta)$.}
\textnormal{Instead of treating the moment restriction as an additional moment condition, one can impose it directly on the estimation problem. The finite sample analog is $-T^{-1}\sum_{t=1}^{T} x_{t}\varepsilon_{R,\,t} \geq 0$, or more explicitly,}
\begin{equation*}
\rho_{R}\sum_{t=1}^{T}x_{t}R_{t-1} + (1 - \rho_{R})\psi_{1}\sum_{t=1}^{T}x_{t}\pi_{t} + (1 - \rho_{R})\psi_{2}\sum_{t=1}^{T}x_{t}^{2} - \sum_{t=1}^{T}x_{t}R_{t} \geq 0,
\end{equation*}
\textnormal{which imposes a linear inequality constraint on $\theta$.}
\end{example}
\begin{example} \label{exmp:slutzky}
\textnormal{Inequality constraints also arise in many demand models. Consider a consumer who chooses her levels of consumption for different goods $j = 1,\dots,\,J$ by maximizing her utility function with respect to her budget constraint. One can show that the demand functions}
\begin{equation*}
D_{j} = D_{j}(p,\,m|\theta), \quad j = 1,\dots,\,J,
\end{equation*}
\textnormal{where $p$ is a price vector, $m$ is income, and $\theta$ are the structural parameters of interest, are not arbitrary. In particular, they must satisfy the budget constraint}
\begin{equation*}
\sum_{j=1}^{J}p_{j}D_{j}(p,\,m|\theta) = m.
\end{equation*}
\textnormal{Furthermore, since they solve a constrained optimization problem, they must satisfy the Slutsky matrix conditions. Let $S$ denote the Slutsky substitution matrix of size $J \times J$, whose generic entry is}
\begin{equation*}
S_{kj} = \frac{\partial D_{j}(p,\,m|\theta)}{\partial p_{k}} + \frac{\partial D_{j}(p,\,m|\theta)}{\partial m} D_{k}(p,\,m|\theta).
\end{equation*}
\textnormal{Economic theory tells us that such a matrix must be symmetric and negative semidefinite. These conditions imply inequality restrictions on the vector of structural parameters $\theta$.}
\end{example}
To measure the accuracy of an estimator $T_{n} = T_{n}(\bm{X}_{n})$ of $\theta$ we will use a known loss function $\ell(\theta,\,T_{n})$. The corresponding risk is just the expected loss
\begin{equation} \label{eq:risk_fun}
R(\theta,\,T_{n}) = \mathbb{E}_{\theta} \ell(\theta,\,T_{n}).
\end{equation}
The most popular loss function in the literature is weighted quadratic loss,
\begin{equation} \label{eq:quad_loss}
\ell(\theta,\,T_{n}) = (T_{n} - \theta)'W(T_{n} - \theta)
\end{equation}
for some weight matrix $W > 0$. The risk associated with \eqref{eq:quad_loss} is simply weighted mean squared error. In general, the choice of a loss function can be motived by an economic application, see \cite{hansen2016} for more examples.
The choice of a loss function plays a crucial role in the shrinkage estimator's behavior since the weights depend on the loss between the unrestricted and the restricted estimators. We specify the following regularity conditions for the loss function.
\begin{assumption} \label{loss_fun}
The loss function $\ell(\theta,\,T_{n})$ satisfies
\begin{itemize}
\item [(a)] $\ell(\theta,\,T_{n}) \geq 0$
\item [(b)] $\ell(\theta,\,\theta) = 0$
\item [(c)]$W(\theta) = \left. \frac{1}{2}\frac{\partial^2}{\partial T_{n} \partial T_{n}'} \ell(\theta,\,T_{n}) \right\vert_{T_{n} = \theta}$ is continuous in a neighborhood of $\theta_{0}$.
\end{itemize}
\end{assumption}
Assumptions \ref{loss_fun}(a) and (b) are standard properties of any loss function. The dominance result of \cite{james_stein1961} hinges on the quadratic loss function, however, our results hold for a more general family of loss functions. Assumption \ref{loss_fun}(c) requires the loss function $\ell(\theta,\,T_{n})$ to have a second derivative with respect to the second argument. This allows for smooth loss functions, like quadratic loss, and excludes non-smooth loss functions, such as absolute value loss.
The choice of a weight matrix also plays an important role. If one sets $W = \mathcal{I}_{m}$, \eqref{eq:quad_loss} becomes unweighted quadratic loss, which is appropriate for cases where all parameters are roughly identically scaled. However, when it is not the case, a weight matrix that renders a loss function which is robust to rotations of the parameter vector $\theta$ is a more plausible choice. We can fulfill the latter task by setting $W = \Omega^{-1}$, where $\Omega^{-1}$ is the inverse of the asymptotic variance of the unrestricted estimator.
\subsection{Shrinkage direction} \label{sec:shr_dir}
The restriction in \eqref{eq:restricted_set} defining the direction of shrinkage, it is the main building block for the construction of our shrinkage estimator. Inequality constraints impose milder restrictions compared to equality constraints, which makes them harder to deal with. Equality restrictions provide the researcher with a particular shrinkage direction, however, with inequality constraints the shrinkage direction depends on the boundary which the true parameter value is close to. This stems from the properties of the inequality constrained estimator, see Section \ref{sec:distribution} for more details.
The researcher usually believes that restrictions are a reasonable simplification of the unrestricted model specification. And it is well known that if restrictions are correct, the restricted estimator renders more efficient estimates. In contrast, if not, the restricted estimates will be biased. In case of equality restrictions, the researcher can easily test them, however, testing inequality restrictions is an onerous task. Rather than testing inequality constraints, we can use them to construct an Inequality Constrained Shrinkage Estimator, and thereby, improve the efficiency of estimates.
\subsection{Local asymptotic framework}
Our estimation framework is based on the belief that the empirical implications of theoretical restrictions are approximately correct. Put differently, it means that the parameter of interest $\theta_{0}$ does not necessarily lie within the restricted space $\Theta_{0}$, but is localized to it. We model that by assuming that the constraints are local to zero, i.e. $r(\theta_{0}) = cn^{-1/2}$, where $c \in \mathbb{R}^{p}$. In this framework $c$ is a slackness, or localizing, parameter which measures the discrepancy between $\theta_{0}$ and $\Theta_{0}$. When $c > 0$, then the constraints are satisfied and not binding, while if $c < 0$, the constraints are violated.\footnote{We are particularly interested in cases when constraints are locally violated. However, the analysis does not depend on the sign of the localizing parameter.} This modeling assumption ensures that the normalized asymptotic distribution of the ICSE is identical to its finite sample distribution under exact normality (see e.g. \citealp{hansen2016}).
We do not consider distant alternatives of the form $r(\theta_{0}) = \kappa_{n}c$, where $\kappa_{n}$ is $O(n^{-b})$ with $b < 1/2$, since we are interested in the asymptotic distribution of the normalized estimator. For simplicity, assume we have one only constraint. If $c < 0$, then $n^{1/2}r(\theta_{0}) = n^{1/2}\kappa_{n}c \rightarrow -\infty$, meaning that the constraint is violated, and we are better off with the restricted estimator. In contrast, if $c > 0$, then $n^{1/2}r(\theta_{0}) = n^{1/2}\kappa_{n}c \rightarrow \infty$, meaning that the constraint is satisfied as a strict inequality, and we should resort to the unrestricted estimator.
\section{Estimation} \label{sec:estimation}
In order to define the shrinkage estimator, we first need to introduce unrestricted and restricted estimators.
The unrestricted estimator $\hat{\theta}_{n}$ of $\theta$ maximizes the objective function over $\theta \in \Theta$
\begin{equation*}
Q_{n}(\hat{\theta}_{n}) = \sup_{\theta \in \Theta} Q_{n}(\theta).\footnote{We can allow for a numerical error by requiring $Q_{n}(\hat{\theta}_{n})$ to be within $o_{p}(1)$ of the global maximum of $Q_{n}(\theta)$, rather than the exact global maximum. This is a common assumption in the extremum estimators literature, yet fairly technical.}
\end{equation*}
The restricted estimator $\tilde{\theta}_{n}$ is defined analogously
\begin{equation*}
Q_{n}(\tilde{\theta}_{n}) = \sup_{\theta \in \Theta_{0}} Q_{n}(\theta).
\end{equation*}
We assume that the maximum is unique so that $\hat{\theta}_{n}$ and $\tilde{\theta}_{n}$ are well-defined.
The shrinkage estimator is defined as a weighted average of the unrestricted and restricted estimators
\begin{equation} \label{eq:shrink_est}
\hat{\theta}^{*}_{n} = \hat{w}_{n} \hat{\theta}_{n} + (1 - \hat{w}_{n}) \tilde{\theta}_{n},
\end{equation}
where the weight is data driven and takes the form
\begin{equation} \label{eq:weight_def}
\hat{w}_{n} = \left( 1 - \frac{\hat{\tau}_{n}}{n\ell(\hat{\theta}_{n},\,\tilde{\theta}_{n})} \right)_{+},
\end{equation}
where $\hat{\tau}_{n} \geq 0$ is the shrinkage parameter which controls the degree of shrinkage and $n\ell(\hat{\theta}_{n},\,\tilde{\theta}_{n})$ is the scaled loss between the unrestricted and restricted estimators. Under the quadratic loss, the latter becomes $n\left(\hat{\theta}_{n} - \tilde{\theta}_{n}\right)'W\left(\hat{\theta}_{n} - \tilde{\theta}_{n}\right)$.
The shrinkage parameter $\hat{\tau}_{n}$ is set to minimize the asymptotic risk of the ICSE. Thus, we allow $\hat{\tau}_{n}$ to be data-dependent and random, however, require it to converge in probability to a non-negative constant.
\begin{assumption} \label{shrink_par}
$\hat{\tau}_{n} \stackrel{p}{\rightarrow} \tau \geq 0$ as $n \rightarrow \infty$.
\end{assumption}
The degree of shrinkage determines an optimal bias variance tradeoff and depends on the ratio of the shrinkage parameter $\hat{\tau}_{n}$ to the loss $n\ell(\hat{\theta}_{n},\,\tilde{\theta}_{n})$. When the restricted estimator is very close to the unrestricted one, i.e. the loss is small, and $\hat{\tau}_{n} > n\ell(\hat{\theta}_{n},\,\tilde{\theta}_{n})$, we put all the weight on the restricted estimator, $\hat{w}_{n} = 0$ and $\hat{\theta}_{n}^{*} = \tilde{\theta}_{n}$. When $\hat{\tau}_{n} < n\ell(\hat{\theta}_{n},\,\tilde{\theta}_{n})$, then $\hat{\theta}_{n}^{*}$ is a weighted average of the restricted and unrestricted estimators. The larger the loss compared to the shrinkage parameter, the more weight we put on the unrestricted estimator. In other words, it means that if the regularization bias is small, we are better off trading it for a reduction in variance.
\section{Asymptotic distribution} \label{sec:distribution}
It is a well-known fact that the asymptotic distribution of the unrestricted extremum estimator is normal (see e.g. \cite{newey_mcfadden1994}), however, the asymptotic distribution of the inequality constrained estimator takes a more complicated form. Obtaining the restricted estimator requires solving an inequality constrained optimization problem, the solution to which depends on which constraints bind. As a result, the asymptotic distribution will take the form of a sum of truncated normal random variables.
We introduce the following regularity conditions.
\begin{assumption} \label{consistency}
~\begin{itemize}
\item [(a)] For some some function $Q(\theta): \Theta \rightarrow \mathbb{R}$, $\sup_{\theta \in \Theta} |Q_{n}(\theta) - Q(\theta)| \rightarrow_{p} 0$;
\item [(b)] For all $\varepsilon > 0$, $\sup_{\theta \in \Theta / N(\theta_{0},\,\varepsilon)} Q(\theta) < Q(\theta_{0})$, where $N(\theta_{0},\,\varepsilon)$ is an $\varepsilon$-neighborhood of $\theta_{0}$.
\end{itemize}
\end{assumption}
Assumption \ref{consistency}(a) ensures uniform convergence of the sample criterion function to the true criterion function. Assumption \ref{consistency}(b) requires the true criterion function to be uniquely maximized at $\theta_{0}$ in its neighborhood. These conditions guarantee that both the unrestricted and restricted estimators are consistent, i.e. $\hat{\theta}_{n} - \theta_{0}$ and $\tilde{\theta}_{n} - \theta_{0}$ are $o_{p}(1)$. Note that consistency does not depend on whether the estimator is restricted or not, the only thing that changes is the parameter space over which an estimator is defined (see e.g. Theorem 9.1 in \citealp{newey_mcfadden1994}).
\begin{assumption} \label{assumption_QA}
~\begin{itemize}
\item [(a)] $\Theta$ is a compact subset of $\mathbb{R}^{m}$;
\item [(b)] $\theta_{0}$ lies in the interior of $\Theta$;
\item [(c)] $Q_{n}(\theta)$ is twice continuously differentiable in a neighborhood $N(\theta_0,\,\varepsilon)$ of $\theta$;
\item [(d)] $n^{1/2} \frac{\partial}{\partial \theta} Q_{n}(\theta_{0}) \rightarrow_{d} G = \mathcal{N}(0,\,\mathcal{V})$ for some nonrandom positive definite matrix $\mathcal{V}$;
\item [(e)] For $\theta \in N(\theta_{0},\,\varepsilon)$ there exists $\mathcal{J}(\theta)$ that is continuous and non-singular at $\theta_{0}$ and \\ $\sup_{\theta \in N(\theta_{0},\,\varepsilon)} \parallel \frac{\partial^{2}}{\partial \theta \partial \theta'} Q_{n}(\theta) - \mathcal{J}(\theta)\parallel \rightarrow_{p} 0$.
\end{itemize}
\end{assumption}
Assumption \ref{assumption_QA} is a standard set of assumptions to ensure asymptotic normality of extremum estimators (see e.g. \citealp{newey_mcfadden1994}). Note that Assumption \ref{assumption_QA}(b) does not imply that $\theta_{0}$ lies in the interior of the restricted set $\Theta_{0}$, and whether $\theta_{0}$ belongs to the interior of $\Theta_{0}$ or not will affect the asymptotic distribution of both the restricted and shrinkage estimators.
\begin{assumption} \label{constr_assumption}
~\begin{itemize}
\item [(a)] $R(\theta)$ is continuous in some neighborhood of $\theta_{0}$;
\item [(b)] $R(\theta_{0})$ has full row rank.
\end{itemize}
\end{assumption}
Assumption \ref{constr_assumption}(a) allows for applying the continuous mapping theorem, and Assumption \ref{constr_assumption}(b) rules out linearly dependent constraints.
\subsection{Solving an asymptotically equivalent problem} \label{subsec:asy_equiv_problem}
The asymptotic behavior of the unrestricted estimator is easily characterized, however, the distribution of the inequality constrained estimator is more complicated. Recall that in order to obtain the restricted estimator, we have to solve the following problem
\begin{equation} \label{eq:original_problem}
\sup_{\theta \in \Theta} Q_{n}(\theta) \quad \text{s.t.} \quad r(\theta) \geq 0.
\end{equation}
Dealing with non-linear inequality constrained optimization problems typically leads to very cumbersome calculations of the first order conditions. However, it turns out that we do not have to solve the original optimization problem. To derive the asymptotic distribution of the constrained estimator, it is sufficient to solve a simpler, asymptotically equivalent problem (see e.g. Section 21.3.2 in \citealp{gourieroux_monfort1995}).
In our asymptotic analysis we follow \cite{andrews1999boundary} and rely on the quadratic approximation of the objective function around the true parameter value. In particular,
\begin{equation} \label{eq:quad_objective}
Q_{n}(\theta) = Q_{nq}(\theta) + \xi_{n}(\theta).
\end{equation}
where
\begin{equation*}
Q_{nq}(\theta) = Q_{n}(\theta_{0}) + \frac{\partial}{\partial \theta'} Q_{n}(\theta_{0}) (\theta - \theta_{0}) + \frac{1}{2}(\theta - \theta_{0})'\frac{\partial^{2}}{\partial \theta \partial \theta'} Q_{n}(\theta_{0}) (\theta - \theta_{0})
\end{equation*}
and $\xi_{n}(\theta)$ is the approximation error. We need to introduce some additional assumptions ensuring that $\xi_{n}(\theta)$ is of the right order, so that the estimator maximizing $Q_{nq}(\theta)$ has the same asymptotic distribution as of the true maximum.
\begin{assumption} \label{stoch_diff}
For all $\delta_{n} \rightarrow 0$,
\begin{equation*}
\sup_{\theta \in \Theta: ||\theta - \theta_{0}|| \leq \delta_{n}} \frac{|\xi_{n}(\theta)|}{(1 + ||n^{1/2}(\theta - \theta_{0})||^{2})} = o_p(1).
\end{equation*}
\end{assumption}
\cite{pollard1985} refers to Assumption \ref{stoch_diff} as stochastic differentiability, which is a weaker condition than $\xi_{n}(\theta)$ converging to 0 due to the presence of the denominator term $(1 + ||n^{1/2}(\theta - \theta_{0})||^{2})$.
Let
\begin{equation*}
\mathcal{J}_{n} \equiv - \frac{\partial^{2}}{\partial \theta \partial \theta'} Q_{n}(\theta_{0}) \quad and \quad Z_{n} \equiv \mathcal{J}_{n}^{-1}n^{1/2} \frac{\partial}{\partial \theta} Q_{n}(\theta_{0}).
\end{equation*}
The quadratic approximation in \eqref{eq:quad_objective} can be rewritten as
\begin{align*}
Q_{nq}(\theta) & = Q_{n}(\theta_{0}) + n^{-1/2} Z_{n}' \mathcal{J}_{n} (\theta - \theta_{0}) - \frac{1}{2} (\theta - \theta_{0})' \mathcal{J}_{n} (\theta - \theta_{0}) \\
& = Q_{n}(\theta_{0}) + \frac{1}{2n} Z_{n}' \mathcal{J}_{n} Z_{n} - \frac{1}{n}q_{n}(n^{1/2}(\theta - \theta_{0})),
\end{align*}
where
\begin{equation*}
q_{n}(\lambda) \equiv \frac{1}{2}(\lambda - Z_{n})' \mathcal{J}_{n} (\lambda - Z_{n}) \quad \text{and} \quad \lambda \in \mathbb{R}^{m}.
\end{equation*}
Note that under Assumption \ref{stoch_diff}, it is sufficient to minimize $q_{n}(n^{1/2}(\theta - \theta_{0}))$ to obtain a maximum of the quadratic approximation of $Q_{n}(\theta)$. When the parameter space is unrestricted, the estimator $\hat{\theta}_{n}$ equals to $\theta_{0} + n^{-1/2}Z_{n}$. Therefore, $n^{1/2}(\hat{\theta}_{n} - \theta_{0}) = Z_{n}$, and $Z_{n}$ determines the asymptotic distribution of the unrestricted estimator. A lemma below establishes the asymptotic distribution of the re-parameterized quadratic criterion function.
\begin{lemma} \label{lemma:quad_approx}
Under Assumptions \ref{consistency}--\ref{constr_assumption},
\begin{equation} \label{eq:crit_limit}
\begin{aligned}
& \quad Z_{n} \hspace{0.1cm} \stackrel{d}{\rightarrow} \hspace{0.1cm} Z = \mathcal{J}^{-1}G, \\
q_{n}(\lambda) \hspace{0.1cm} \stackrel{d}{\rightarrow} \hspace{0.1cm} &\; q(\lambda) \equiv \frac{1}{2}(\lambda - Z)'\mathcal{J}(\lambda - Z) \quad \forall \lambda \in \mathbb{R}^{m}.
\end{aligned}
\end{equation}
\end{lemma}
Since the restricted estimator $\tilde{\theta}_{n}$ is consistent, its asymptotic distribution depends only on the features of the parameter space around the true parameter value $\theta_{0}$. We use the mean value expansion to approximate the constraints $r(\theta)$ around $\theta_{0}$,
\begin{equation*}
r(\theta) = r(\theta_{0}) + R(\bar{\theta})(\theta - \theta_{0}) = c + R(\bar{\theta})n^{1/2}(\theta - \theta_{0})\geq 0,
\end{equation*}
where $\bar{\theta}$ lies on a segment between $\theta$ and $\theta_{0}$.\footnote{Essentially this approach is the same as approximating the restricted space by a cone of tangents (see e.g. \citealp{chernoff1954}, \citealp{feder1968}, and \citealp{andrews1999boundary}).} Since $\bar{\theta}_{n}$ lies on a segment between $\tilde{\theta}_{n}$ and $\theta_{0}$, under Assumptions \ref{consistency} and \ref{constr_assumption}, $R(\bar{\theta}_{n}) = R(\theta_{0}) + o_{p}(1)$. Let $R \equiv R(\theta_{0})$. As shown in Lemma \ref{lemma:distr_re} below, the asymptotic distribution of $n^{1/2}(\hat{\theta}_{n} - \theta_{0})$ is given by the distribution of
\begin{equation} \label{eq:reparam_problem}
\begin{aligned}
\tilde{\lambda} = \operatornamewithlimits{argmin}_{\lambda \in \Lambda_{c}}\; q(\lambda)
\end{aligned}
\end{equation}
where $\Lambda_{c} \equiv \{\lambda \in \mathbb{R}^{m}: c + R\lambda \geq 0\}$. By approximating the objective function with a quadratic counterpart and linearizing the constraints, we collapsed a potentially highly non-linear problem \eqref{eq:original_problem} to a simple quadratic programming problem.
There are $p$ inequality constraints which form $2^{p}$ different possible combinations of binding and non-binding constraints.\footnote{One can think of these combinations as possible boundaries of the restricted parameter space $\Theta_{0}$.} For each such combination the asymptotic distribution of the restricted estimator is simply a projection of the asymptotic limit of the unrestricted estimator $Z$ on the corresponding boundary. This is exactly the intuition in \cite{andrews1999boundary}, where he shows that under the standard asymptotics the asymptotic distribution of the extremum estimator, when the true parameter value is on a boundary, depends on binding constraints.
Let us introduce some notation simplifying the exposition. Let $L(\iota)$ be a linear subspace of the form $L(\iota) \equiv \{ l \in \mathbb{R}^{m}: c_{\iota} + R_{\iota}l = 0 \}$, where $\iota = 1,\dots,\,2^{p}$ represents one of the possible combinations of binding constraints. Let $\iota = 1$ denote the case when none of the constraints bind. Let $R_{\iota}$ consist of the rows of the Jacobian matrix $R$ corresponding to binding constraints indexed by $\iota$. By analogy, $\tilde{\mu}_{n,\iota}$ denotes a sub-vector of $\tilde{\mu}_{n}$ with entries corresponding to binding constraints indexed by $\iota$. Note that we also have to index the slackness parameter, as only the entries corresponding to binding constraints $c_{\iota}$ will affect the asymptotic distribution.
\begin{lemma} \label{lemma:distr_re}
Suppose that Assumptions \ref{consistency}--\ref{stoch_diff} hold. Then, the asymptotic distribution of the constrained estimator takes the form
\begin{equation} \label{eq:lambda_bind_gen}
n^{-1/2}(\tilde{\theta}_{n} - \theta_{0}) \hspace{0.1cm} \stackrel{d}{\rightarrow} \hspace{0.1cm} \tilde{\lambda} \equiv Z - \sum_{\iota=2}^{2^{p}} P_{L(\iota)}(Z + h_{\iota}) \mathds{1}\{ \tilde{\mu}_{\iota} > 0, \; \tilde{\mu}_{-\iota} \leq 0 \},
\end{equation}
where $\tilde{\mu} = -(R\mathcal{J}^{-1}R')^{-1}(RZ + c)$ is the vector of Kuhn-Tucker multipliers for problem \eqref{eq:reparam_problem},
\begin{equation}
P_{L(\iota)} \equiv \mathcal{J}^{-1}R_{\iota}'\left(R_{\iota}\mathcal{J}^{-1}R_{\iota}'\right)^{-1}R_{\iota}
\end{equation}
is the projection on the linear subspace $L(\iota)$, $h_{\iota} \equiv R_{\iota}^{-1}c_{\iota}$ is the re-parameterized slackness parameter, and $R_{\iota}^{-1}$ is the right inverse of $R_{\iota}$.
\end{lemma}
Note that the distribution in \eqref{eq:lambda_bind_gen} is non-normal and depends on the re-parametarized slackness parameter $h$. The distribution takes the form of a sum of truncated normal random variables. Notice that the indicator functions are random: they depend on the asymptotic distribution of the Kuhn-Tucker multipliers.
The slackness parameter enters the distribution through both the asymptotic bias term $P_{L(\iota)}h_{\iota}$ and the distribution of the Kuhn-Tucker multipliers $\tilde{\mu}$. From \eqref{eq:crit_limit} it follows that if $h_{j} \rightarrow \infty$, then $\tilde{\mu}_{j} \rightarrow -\infty$, implying that the $j^{th}$ constraint is not binding. If, on the contrary, $h_{j} \rightarrow -\infty$, then $\tilde{\mu}_{j} \rightarrow \infty$, resulting into the $j^{th}$ constraint being binding.
The summation starts from $\iota = 2$ since we do not have to project the unrestricted estimator on any subspace when none of the constrains bind. Despite the seemingly complex expression, the basic intuition behind this formula is surprisingly simple. The asymptotic distribution of the inequality constrained estimator is just a projection of the asymptotic limit of the unconstrained estimator onto a boundary defined by the corresponding set of binding constraints.
The following theorem summarizes the analysis above and presents the asymptotic distributions of the unrestricted, restricted, and shrinkage estimators.
\begin{theorem} \label{thm:asy_distr}
Under Assumptions \ref{loss_fun}--\ref{stoch_diff},
\begin{align}
& n^{1/2}(\hat{\theta}_{n} - \theta_{0}) \hspace{0.1cm} \stackrel{d}{\rightarrow} \hspace{0.1cm} Z \sim \mathcal{N}(0,\,\Omega), \label{eq:ur_distr} \\
& n^{1/2}(\tilde{\theta}_{n} - \theta_{0}) \hspace{0.1cm} \stackrel{d}{\rightarrow} \hspace{0.1cm} \tilde{\lambda} \equiv Z - \sum_{\iota=2}^{2^{p}} P_{L(\iota)}(Z + h_{\iota}) \mathds{1}\{ \tilde{\mu}_{\iota} > 0, \; \tilde{\mu}_{-\iota} \leq 0 \}, \label{eq:r_distr} \\
& n\ell(\hat{\theta}_{n},\,\tilde{\theta}_{n}) \hspace{0.1cm} \stackrel{d}{\rightarrow} \hspace{0.1cm} \xi \equiv \sum_{\iota=2}^{2^{p}} (Z + h_{\iota})'P_{L(\iota)}' W P_{L(\iota)}(Z + h_{\iota}) \mathds{1}\{ \tilde{\mu}_{\iota} > 0, \; \tilde{\mu}_{-\iota} \leq 0 \}, \label{eq:loss_distr} \\
& \hat{w}_{n} \hspace{0.1cm} \stackrel{d}{\rightarrow} \hspace{0.1cm} w = \left(1 - \frac{\tau}{\xi}\right)_{+}. \label{eq:weight_distr}
\end{align}
The asymptotic distribution of the inequality constrained shrinkage estimator is
\begin{align} \label{eq:shrink_distr}
& n^{1/2}(\hat{\theta}_{n}^{*} - \theta_{0}) \hspace{0.1cm} \stackrel{d}{\rightarrow} \hspace{0.1cm} wZ + (1 - w)\tilde{\lambda}. \qquad\qquad\qquad\qquad\qquad\qquad\qquad
\end{align}
\end{theorem}
\vspace{1 em}
\section{Asymptotic Risk} \label{sec:asymptotic_risk}
In practice obtaining the restricted estimator still requires solving a potentially complicated non-linear problem. This suggests that having an analytical closed form solution is extremely unlikely. Even if it is possible to derive an analytical solution, this solution will take a complex form, and the ICSE will inherit it. As a result, calculating its finite sample risk may be infeasible. However, we know the asymptotic distribution of the ICSE, which means we can use the asymptotic risk to get a reasonable approximation of the finite sample risk.
Since the ICSE may not have a sufficient number of finite moments, to ensure existence we use an asymptotic trimmed loss. Let $T = \{T_{n}\}_{n=1}^{\infty}$ denote a sequence of estimators. The asymptotic risk of the estimator sequence $T$ is defined as
\begin{equation} \label{eq:asy_risk}
\rho(h,\,T) = \lim_{\zeta \rightarrow \infty} \liminf_{n \rightarrow \infty} \mathbb{E}_{\theta_{0}} \min \left[ n\ell\left(\theta_{0},\,T_{n}\right),\,\zeta \right].
\end{equation}
The loss function is trimmed at $\zeta$, however, the trimming becomes negligible in large samples as $\zeta \rightarrow \infty$ with $n \rightarrow \infty$.
\cite{hansen2016} shows that whenever the loss function is locally quadratic, i.e. satisfies Assumption \ref{loss_fun}, the asymptotic risk, defined in \eqref{eq:asy_risk}, of an arbitrary estimator $T_{n}$, such that $n^{1/2}(T_{n} - \theta_{0}) \hspace{0.1cm} \stackrel{d}{\rightarrow} \hspace{0.1cm} \psi$, where $\psi$ is some random variable, can be calculated as
\begin{equation} \label{eq:quad_risk}
\rho(h,\,T) = \mathbb{E}[\psi'W\psi].
\end{equation}
Equation \eqref{eq:quad_risk} allows us to calculate the asymptotic risk of the unrestricted and shrinkage estimators as expected weighted quadratic loss. Note that $n^{1/2}(\hat{\theta}_{n} - \theta_{0}) \hspace{0.1cm} \stackrel{d}{\rightarrow} \hspace{0.1cm} Z$, hence, the asymptotic risk of the unrestricted estimator is
\begin{equation} \label{eq:asy_risk_ur}
\rho(h,\,\hat{\theta}_{n}) = \mathbb{E}[Z'WZ] = tr(W\mathbb{E}[ZZ']) = tr(W\Omega).
\end{equation}
Define an $m \times m$ matrix $A_{L(\iota)} \equiv W^{1/2\prime}\Omega P_{L(\iota)}'W^{1/2}$, let $\phi_{\max}(A_{L(\iota)})$ denote its largest eigenvalue.
The following theorem establishes the main result of the paper.
\begin{theorem} \label{thm:risk_bound}
Under Assumptions \ref{loss_fun}--\ref{stoch_diff}, if
\begin{equation} \label{eq:tau_bounds}
0 < \tau \leq \sum_{\iota = 2}^{2^{p}} 2\left(tr(A_{L(\iota)}) - 2\phi_{\max}(A_{L(\iota)})\right) \gamma_{\iota},
\end{equation}
where
\begin{equation} \label{eq:weights}
\gamma_{\iota} \equiv \frac{\mathbb{E}[\xi_{L(\iota)}^{-1}]\mathbb{P}(\tilde{\mu}_{\iota} > 0,\,\tilde{\mu}_{-\iota} \leq 0)}{\sum_{\iota = 2}^{2^{p}}\mathbb{E}[\xi_{L(\iota)}^{-1}]\mathbb{P}(\tilde{\mu}_{\iota} > 0,\,\tilde{\mu}_{-\iota} \leq 0)},
\end{equation}
then for any $h$
\begin{equation} \label{eq:risk_opt_bound}
\rho(h,\,\hat{\theta}^{*}_{n}) < \rho(h,\,\hat{\theta}_{n}).
\end{equation}
\end{theorem}
\vspace{1 em}
Equation \eqref{eq:risk_opt_bound} shows that the ICSE has strictly lower asymptotic risk than that of the unrestricted estimator for all values of the slackness parameter $h$, given that the shrinkage parameter $\tau$ satisfies the restriction \eqref{eq:tau_bounds}.
The explicit risk bound for the ICSE is
\begin{equation} \label{eq:explicit_risk_bound}
\rho(h,\,\hat{\theta}^{*}_{n}) < tr(W\Omega) - \tau\sum_{\iota = 2}^{2^{p}} \mathbb{E} \left[ \frac{
2\left(tr(A_{L(\iota)}) - 2\phi_{\max}(A_{L(\iota)})\right) - \tau}{\xi_{L(\iota)}} \right] \mathbb{P}(\tilde{\mu}_{\iota} > 0,\,\tilde{\mu}_{-\iota} \leq 0)
\end{equation}
Since the bound in \eqref{eq:explicit_risk_bound} is quadratic in the shrinkage parameter $\tau$, there exists a unique optimal level of shrinkage $\tau^{*}$ that minimizes this bound,
\begin{equation} \label{eq:opt_shrink}
\tau^{*} = \sum_{\iota = 2}^{2^{p}} \left(tr(A_{L(\iota)}) - 2\phi_{\max}(A_{L(\iota)})\right) \gamma_{\iota}.
\end{equation}
From \eqref{eq:weights} it follows that $\gamma_{\iota} \rightarrow 0$ when either the probability of $\iota^{th}$ event $\mathbb{P}(\tilde{\mu}_{\iota} > 0,\,\tilde{\mu}_{-\iota} \leq 0)$ is close to zero, or the expected inverse loss $\mathbb{E}[\xi^{-1}_{L(\iota)}]$ is approaching zero. The optimal shrinkage parameter puts more weight on events that are more likely to happen and on events where the restricted parameter is close to the unrestricted one. The behavior of $\gamma_{i}$ is ambiguous when $\mathbb{P}(\tilde{\mu}_{\iota} > 0,\,\tilde{\mu}_{-\iota} \leq 0)$ goes to zero and $\mathbb{E}[\xi^{-1}_{L(\iota)}]$ approaches infinity.
When $W = \Omega^{-1}$, \eqref{eq:tau_bounds} simplifies to
\begin{equation} \label{eq:canonical_tau_bounds}
0 < \tau \leq 2\left(\sum_{\iota = 2}^{2^{p}}p_{\iota}\gamma_{\iota} - 2\right),
\end{equation}
which leads to
\begin{equation} \label{eq:canonical_opt}
\tau^{*} = \sum_{\iota = 2}^{2^{p}} p_{\iota}\gamma_{\iota} - 2.
\end{equation}
When $W = \Omega^{-1}$, $tr(A_{L(\iota)}) = tr(W\Omega P_{L(\iota)}') = tr(P_{L(\iota)}') = p_{\iota}$, and $\phi_{\max}(A_{L(\iota)}) = 1$, which gives condition \eqref{eq:canonical_tau_bounds}. This restriction on the shrinkage parameter has the same form as the classical James-Stein condition, $0 < \tau \leq 2(m-2)$, where $m$ is the dimension of the parameter of interest. As long as $m > 2$, the James-Stein estimator will dominate the unrestricted estimator in terms of asymptotic risk. In case of the ICSE, condition \eqref{eq:canonical_tau_bounds} requires $\sum_{\iota = 2}^{2^{p}} p_{\iota}\gamma_{\iota} > 2$. This means that the ``expected'' number of binding constraints must be greater than two for the ICSE to dominate. Constraints that are more likely to bind tell us which boundary of the restricted parameter space we are shrinking to, i.e. they determine the direction of shrinkage.
\section{Data-dependent weights} \label{sec:weight}
As it is pointed out by \cite{hjort_claeskens2003}, model averaging (and shrinkage) optimal weights cannot be consistently estimated in the local asymptotic framework since localizing parameters are $O(n^{1/2})$. And the ICSE is not an exception. Since the weights $\{\gamma_{\iota}\}_{\iota=2}^{2^p}$ depend on the localizing parameter, which is unknown, the optimal shrinkage parameter in \eqref{eq:opt_shrink} is infeasible. Furthermore, the localizing parameter $h$, which is a transformation of the original localizing parameter $c$, cannot be consistently estimated under the local asymptotic framework.
The weights $\{\gamma_{\iota}\}_{\iota=2}^{2^p}$ depend on the localizing parameter through the Kuhn-Tucker multipliers. The distribution of the Kuhn-Tucker multipliers is given by
\begin{equation*}
\tilde{\mu} = -(R\mathcal{J}^{-1}R')^{-1}R(Z + h) = -(R\mathcal{J}^{-1}R')^{-1}(RZ + c) \sim \mathcal{N}(\Psi(c),\,\Xi),
\end{equation*}
where $\Psi(c) = -(R\mathcal{J}^{-1}R')^{-1}c$ and $\Xi = (R\mathcal{J}^{-1}R')^{-1}R\Omega R'(R\mathcal{J}^{-1}R')^{-1}$. We observe that the mean of $\tilde{\mu}$ depends on the localizing parameter, thus, the distribution cannot be consistently estimated, as well as the corresponding probabilities. As a result, the optimal shrinkage parameter is infeasible.
A common approach in the literature is to obtain an asymptotically unbiased estimator of the localizing parameter $c$ (see e.g. \citealp{liu2015}). In our case $\hat{c}_{n} = n^{1/2}r(\hat{\theta}_{n})$ is an asymptotically unbiased estimator of $c$. To see this, approximate $\hat{c}_{n}$ around the true parameter value $\theta_{0}$ using the first-order Taylor expansion,
\begin{equation*}
\begin{aligned}
\hat{c}_{n} = n^{1/2}r(\hat{\theta}_{n}) & = n^{1/2}r(\theta_{0}) + n^{1/2}R(\theta_{0})(\hat{\theta}_{n} - \theta_{0}) + o_p(1) \\
& \hspace{0.1cm} \stackrel{d}{\rightarrow} \hspace{0.1cm} c + RZ \sim \mathcal{N}(c,\,R\Omega R').
\end{aligned}
\end{equation*}
Note that without the normalization a simple plug-in estimator $r(\hat{\theta_{n}})$ is just $O_{p}(1)$.
We propose to use a plug-in estimator of the optimal shrinkage parameter, $\hat{\tau}_{n}^{*} \equiv \tau^{*}(\hat{c}_{n})$. We can replace $R$, $\mathcal{J}$, and $\Omega$ with their consistent estimators $\hat{R}_{n} = R(\hat{\theta}_{n})$, $\hat{\mathcal{J}}_{n}$, and $\hat{\Omega}_{n} = \hat{\mathcal{J}}_{n}^{-1}\hat{\mathcal{V}}_{n}\hat{\mathcal{J}}_{n}^{-1}$. A consistent weighting matrix estimate, $\hat{W}_{n}$, can either be constructed from a specific context (e.g. an identity matrix, $\hat{W}_{n} = \mathcal{I}_{n}$) or as the second derivative of the loss function, i.e. $\hat{W}_{n} = W(\hat{\theta}_{n})$.
We can then estimate $\tau^{*}$ by
\begin{equation} \label{eq:feasible_icse}
\hat{\tau}_{n}^{*} = \sum_{\iota=2}^{2^p} \left(tr(\hat{A}_{n,L(\iota)}) - 2\phi_{max}(\hat{A}_{n,L(\iota)})\right) \hat{\gamma}_{n,\iota},
\end{equation}
where $\hat{A}_{n,L(\iota)} = \hat{W}_{n}^{1/2\prime}\hat{\Omega}_{n}\hat{R}_{n,\iota}'(\hat{R}_{n,\iota}\hat{\mathcal{J}}_{n}^{-1}\hat{R}_{n,\iota}')^{-1}\hat{R}_{n,\iota}\hat{\mathcal{J}}_{n}^{-1}\hat{W}_{n}^{1/2\prime}$ and the weights are constructed as
\begin{equation} \label{eq:feasible_gamma}
\hat{\gamma}_{n,\iota} = \frac{\hat{\mathbb{E}}^{-1}[\xi_{L(\iota)}]\hat{\mathbb{P}}(\tilde{\mu}_{\iota} > 0,\,\tilde{\mu}_{-\iota} \leq 0)}{\sum_{\iota = 2}^{2^{p}}\hat{\mathbb{E}}^{-1}[\xi_{L(\iota)}]\hat{\mathbb{P}}(\tilde{\mu}_{\iota} > 0,\,\tilde{\mu}_{-\iota} \leq 0)}.
\end{equation}
In general, $\xi_{L(\iota)}$ follows a generalized $\chi^{2}$ distribution, which makes estimating its first inverse moment an extremely onerous task.\footnote{For more details on the calculation of inverse moments of the generalized $\chi^2$ distribution see e.g. \cite{jones1986}.} Instead, we proxy $\hat{\mathbb{E}}[\xi_{L(\iota)}^{-1}]$ with $\hat{\mathbb{E}}^{-1}[\xi_{L(\iota)}]$, which tends to work well in practice. We can consistently estimate the expected loss $\mathbb{E}[\xi_{L(\iota)}]$ by
\begin{equation*}
\hat{\mathbb{E}}[\xi_{L(\iota)}] = n(\hat{\theta}_{n} - \tilde{\theta}_{n,\iota})'\hat{W}_{n}(\hat{\theta}_{n} - \tilde{\theta}_{n,\iota}),
\end{equation*}
where $\tilde{\theta}_{n,\iota}$ is the equality constrained estimator given the constraints indexed by $\iota$. Probability estimates $\hat{\mathbb{P}}(\tilde{\mu}_{\iota} > 0,\,\tilde{\mu}_{-\iota} \leq 0)$ are based on the feasible distribution of the Kuhn-Tucker multipliers $\mathcal{N}(\hat{\Psi}_{n},\,\hat{\Xi}_{n})$, where $\hat{\Psi}_{n} = -(\hat{R}_{n}\hat{\mathcal{J}}_{n}^{-1}\hat{R}_{n}')^{-1}\hat{c}_{n}$ and $\hat{\Xi}_{n} = (\hat{R}_{n}\hat{\mathcal{J}}_{n}^{-1}\hat{R}_{n}')^{-1}\hat{R}_{n}\hat{\Omega}_{n}\hat{R}_{n}'(\hat{R}_{n}\hat{\mathcal{J}}_{n}^{-1}\hat{R}_{n}')^{-1}$.
Note, $\{\hat{\gamma}_{n,\iota}\}_{\iota=2}^{2^{p}}$ are not consistent estimates, since they do not converge in probability to their corresponding true values. Instead, they converge in distribution to random limits, which implies that the plug-in estimator of the shrinkage parameter $\hat{\tau}_{n}^{*} \stackrel{p}{\not\rightarrow}\tau^{*}$. Thus, the proposed feasible estimator \eqref{eq:feasible_icse} is not optimal in the sense that it uses the feasible data-driven weight that does not converge in probability to the optimal one. As a result, the dominance over the unrestricted estimator is not guaranteed. Despite that, in the following sections we show that the feasible estimator works well in practice.
\section{Monte Carlo Study} \label{sec:mc}
We demonstrate the finite sample performance of the ICSE in the following numerical simulation. Consider a following linear model. For $i = 1,\dots,\,n$,
\begin{equation*}
y_{i} = x'_{1i}\theta_{1} + x'_{2i}\theta_{2} + \varepsilon_{i}.
\end{equation*}
The regressors $x_{1i}$ and $x_{2i}$ are $k_1 \times 1$ and $k_2 \times 1$, respectively. The vector of regressors, $x_{i}$, is distributed $\mathcal{N}(0,\,\Sigma)$, where $\Sigma_{jj} = 1$ and $\Sigma_{jk} = 0.5$ for $j \neq k$, and the error term, $\varepsilon_{i}$, is $\mathcal{N}(0,\,1)$. The goal is to estimate marginal effects under the belief that $\theta$ may be close to $\Theta_{0} = \{\theta \in \mathbb{R}^{k_1+k_2}:\theta_{1} \geq 0,\,\theta_{2} = 0\}$. For simplicity, in estimation we use a quadratic loss function.
Let $\hat{\theta}_{n}$ denote the unrestricted OLS with $\hat{\Omega}_{n}$ being a consistent estimate of its asymptotic covariance matrix of $n^{1/2}(\hat{\theta}_{n} - \theta_{0})$. Let $\tilde{\theta}_{n}$ be the restricted OLS under $\theta_{1} \geq 0$ and $\theta_{2} = 0$.
We compare the performance of five different estimators of $\theta$. The first is $\hat{\theta}_{n}$, the unrestricted OLS estimator. The second is $\tilde{\theta}_{n}$, the restricted OLS estimator. The third estimator is the generalized James-Stein estimator of \cite{hansen2016}
\begin{equation*}
\hat{\theta}^{JS}_{n} = \hat{w}_{n}\hat{\theta}_{n}, \quad \hat{w}_{n} = \left(1 - \frac{k_1 + k_2 - 2}{n\hat{\theta}_{n}'\hat{\Omega}_{n}^{-1}\hat{\theta}_{n}}\right)_{+},
\end{equation*}
which shrinks both $\theta_1$ and $\theta_2$ to zero.
The fourth estimator is the Empirical Bayes estimator $\hat{\theta}_{n}^{EB}$, which assumes the truncated normal prior $\theta|\nu \sim \mathcal{N}(0,\,1/\nu)\mathds{1}\{\theta \geq 0\}$, where $\nu$ is a hyper parameter that tells us how much weight to put on $\theta$ being equal to zero.\footnote{Further details can be found in Appendix \ref{app:eb}.} The higher the value of $\nu$, the more concentrated is the prior around zero, hence, the more mass is put on zero. The motivation for this estimator comes from the fact that the James-Stein estimator can be represented as an Empirical Bayes estimator (see e.g. \citealp{efron_morris1972a}).
The last estimator is the feasible ICSE, which takes the same form but with the weight
\begin{equation*}
\hat{w}^{*}_{n} = \left(1 - \frac{\sum_{\iota=1}^{2^{k_1}}p_{\iota}\hat{\gamma}_{n,\iota} - 2}{n(\hat{\theta}_{n} - \tilde{\theta}_{n})'\hat{\Omega}_{n}^{-1}(\hat{\theta}_{n} - \tilde{\theta}_{n})}\right)_{+},
\end{equation*}
where $p_{\iota}$ is the total number of binding constraints in $\iota$ case. Note, since there are two equality constraints, if none of the inequality constraints bind, $\iota = 1$, $p_1 = 2$. The weights $\{\hat{\gamma}_{n,\iota}\}_{\iota=1}^{2^{k_1}}$ are estimated by \eqref{eq:feasible_gamma}.
The estimators are compared by the mean square error (MSE), which is calculated based on $N = 2,000$ replications. For the ease of exposition, we normalize the MSE of the unrestricted estimator to be equal to one so that the MSE of other estimators are given relative to the MSE of the unrestricted one.
We set the regression coefficients as $\theta_1 = (1,\,1,\,1,\,b,\dots,\,b)$, and $\theta_{2} = (c,\,c,\dots,\,c)'$. Thus, the remaining control parameters in the model are $k_1$, $b$, $c$, and $n$. The value of $b$ allows us to control the strength of the inequality constraints, i.e. whether they are satisfied or not, and $c$ controls the strength of the equality constraints.
Note that the inequality constraints do not change simultaneously with $b$. When $b$ is negative, the first three constraints are satisfied, while the remaining $k_1 - 3$ constraints are violated. As a result, shrinking towards inequality constraints is fundamentally different from shrinking towards equality constraints.
In Figure \ref{fig:mse_lin2_no_ma}, we display the results for $n = \{200,\,500\}$, $k_1 = \{5,\,7,\,10\}$, and vary $b$ on a $100$-point equispaced grid from $-0.5$ to $0.5$. We set $c=0$ so that the equality constraints are satisfied.
First, the feasible ICSE dominates the unrestricted estimator, while the restricted estimator along with the EB estimator do worse than the unrestricted one when the constraints are violated. Since the posterior is truncated at zero, the EB estimates of $\theta_1$ are always positive, which explains the result.
\begin{figure}[!ht]
\centering
\includegraphics[scale=0.685]{mse_lin2_full.pdf}
\caption[MC results]{\textbf{MC results.} This figure shows normalized MSEs for different combinations of $k_1 = \{5,\,7,\,10\}$ and $n = \{200,\,500\}$.}
\label{fig:mse_lin2_no_ma}
\end{figure}
Second, we observe that the James-Stein estimator exhibits almost no improvement upon the unrestricted estimator. This behavior is expected since the James-Stein estimator shrinks all the constraints towards zero, which is fundamentally different from shrinking towards inequalities. As a result, the James-Stein estimator puts almost no weight on the restricted estimator. If the shrinkage direction is chosen poorly, it will lead to a large bias resulting into poor overall performance. Thus, the shrinkage gains are guaranteed only if the shrinkage direction is chosen properly.
When $b < 0$, the restricted and EB estimators perform worse than the unrestricted one, while the feasible ICSE achieves significant MSE reduction gains. When $b$ approaches zero, the constrained estimator starts to dominate the feasible ICSE. When $b > 0$, the constrained estimator dominates both shrinkage estimators, however, the EB estimator achieves lower MSE when $b$ is slightly greater than zero. As $b$ grows, the EB estimator converges to the unrestricted estimator. When the number of observations increases, the prior gets less weight pushing the drop in the MSE closer to $b = 0$. Notice that the MSE of the restricted estimator does not converge to the one of the unrestricted. Since $c = 0$, the inequality constrained estimator is more accurate than the unrestricted one, which explains the result. Moreover, as the number of inequality constraints grows, the difference between the unconstrained and constrained estimators vanishes, resulting into lower MSE gains of the constrained estimator over the unconstrained one.
Finally, when the number of inequality constraints increases, the shrinkage effect of the feasible ICSE and EB estimators becomes more prominent, which supports the theoretical findings.
\section{Empirical Application: Demand Estimation under the Slutsky Restriction} \label{sec:slutsky}
In our empirical application we consider consumer demand estimation under the Slutsky restriction (see Example \ref{exmp:slutzky} for more details). In this application we build on literature on demand estimation under shape restrictions, especially on the recent results by \cite{blundell2012}, \cite{dette2016}, and \cite{blundell2017}.
Our goal is to estimate price and income elasticities of gasoline demand for different income levels. Slutsky condition is an inequality constraint on the demand function ensuring that the compensated own-price elasticities are negative. Despite the fact that in theory consumer choices should abide the Slutsky restriction, in the data we might find evidence suggesting otherwise. For example, if gasoline prices are too high and households anticipate them to rise further, then households will tend to buy more gasoline now and store it for future use resulting in positive compensated price elasticity, which violates the Slutsky restriction. That is exactly where we expect shrinkage gains. Implementation details can be found in Appendix \ref{app:ea_details}.
We use the same data and sample construction as \cite{blundell2017}, which we briefly describe here.\footnote{Further details on sample construction can be found in Section IV.A of \cite{blundell2017}. A more detailed description of the NHTS dataset is presented in Section 3 of \cite{blundell2012}.} The data are from the 2001 National Household Travel Survey (NHTS). The sample is constructed to reduce heterogeneity by restricting the analysis to households with a white respondent, two or more adults, at least one child under age 16, and at least one driver. Households in the most rural areas and in Hawaii are excluded from the sample, as well as are households with missing relevant variables or without a gasoline based vehicle. The resulting sample contains 3,640 observations, where the key variables of interest are gasoline demand, price of gasoline, and household income.
\begin{figure}[!ht]
\centering
\begin{subfigure}{0.75\textwidth}
\includegraphics[width=\textwidth,height=6cm]{ICSE_demand_estimates_high_income_20_grids.pdf}
\caption{High Income}
\label{subfig:icse_gas_high}
\end{subfigure}
\begin{subfigure}{0.75\textwidth}
\includegraphics[width=\textwidth,height=6cm]{ICSE_demand_estimates_medium_income_20_grids.pdf}
\caption{Medium Income}
\label{subfig:icse_gas_med}
\end{subfigure}
\begin{subfigure}{0.75\textwidth}
\includegraphics[width=\textwidth,height=6cm]{ICSE_demand_estimates_low_income_20_grids.pdf}
\caption{Low Income}
\label{subfig:icse_gas_low}
\end{subfigure}
\vspace{1em}
\caption[Price and income elasticity estimates]{\textbf{Price and income elasticity estimates.} This figure shows the unrestricted, restricted, and ICSE estimates of price and income elasticities.}
\label{fig:icse_gas}
\end{figure}
We demonstrate estimates for low, medium, and high income level groups which correspond to the first, second, and third quartile, respectively. As a base estimator we use the local linear regression (LLR) with 20 grid points in the observe range of values for the log price. We set the bandwidth for log price and log income using the rule of thumb to their respective standard deviations. Further implementation details are left for Appendix \ref{app:ea_details}.
Figure \ref{fig:icse_gas} plots the unrestricted, restricted, and ICSE estimates of price and income elasticities as functions of price, across the income levels. Degree of shrinkage differs across income groups. We estimate the weight on the unrestricted estimator $\hat{w}$ to be 0 for low income group, $0.25$ for medium income group, and $0.75$ for high income group. Thus, consumers from higher income groups are more likely to have upward sloping demand curves, which is consistent with the results in \cite{blundell2012}. However, the Empirical Bayes estimates of \cite{kasy2018}, based on the local linear quantile regression, suggest to shrink more towards the restricted estimates for all income groups. The reason the estimates differ is due to the fact that ICSE shrinks all components of $\hat{\beta}$ by the same factor $\hat{w}$, while the EB estimator provides component-wise shrinkage with different shrinkage factors (for more details see Section 4.1 in \citealp{kasy2018}).
\section{Conclusion} \label{sec:conclusion}
In this paper we have shown how to shrink extremum estimators towards theoretical restrictions in form of inequality constraints. The ICSE asymptotically uniformly dominates the unrestricted estimator. The shrinkage direction depends only on the binding constraints rendering it \textit{ex ante} unknown to the researcher, which is the main difference compared to shrinking towards equality constraints.
An important caveat, however, is that due to the presence of localizing parameters that cannot be consistently estimated we cannot guarantee the risk dominance result in finite samples, which is a common problem in frequentist model averaging and shrinkage literatures. One possible improvement would be to establish uniform dominance of the ICSE, but we leave this for future research.
\printbibliography
\newpage