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.
113,741 characters
ReLU-Based and DNN-Based Generalized Maximum Score Estimators
\global\long
\global\long
\global\long
\global\long
\global\long
\global\long
\global\long
\global\long
\global\long
\global\long
\global\long
\global\long
\global\long
\global\long
\global\long
\global\long
\global\long
\global\long
\global\long
\global\long\global\long\global\long\global\long
\global\long\global\long
\global\long\global\long
\global\long
\global\long
\global\long\global\long
\global\long\global\long\global\long\global\long\global\long
\global\long
\global\long\def\mc#1{\mathscr{#1}}
\global\long\global\long
\global\long
\global\long
\global\long
\global\long
\global\long
\global\long
\global\long
\global\long
\global\long
\global\long
\global\long
\global\long
\global\long\def\abs#1{\left|#1\right|}
\global\long\def\norm#1{\left\Vert #1\right\Vert }
\global\long\def\rest#1{\left.#1\right|}
\global\long\def\bracket#1#2{\left\langle #1\middle\vert#2\right\rangle }
\global\long\def\sandvich#1#2#3{\left\langle #1\middle\vert#2\middle\vert#3\right\rangle }
\global\long\def\turd#1{\frac{#1}{3}}
\global\long
\global\long\def\sand#1{\left\lceil #1\right\vert }
\global\long\def\wich#1{\left\vert #1\right\rfloor }
\global\long\def\sandwich#1#2#3{\left\lceil #1\middle\vert#2\middle\vert#3\right\rfloor }
\global\long\def\abs#1{\left|#1\right|}
\global\long\def\norm#1{\left\Vert #1\right\Vert }
\global\long\def\rest#1{\left.#1\right|}
\global\long\def\inprod#1{\left\langle #1\right\rangle }
\global\long\def\ol#1{\overline{#1}}
\global\long\def\ul#1{\underline{#1}}
\global\long\def\td#1{\tilde{#1}}
\global\long\def\bs#1{\boldsymbol{#1}}
\global\long
\global\long
\global\long
\global\long
\global\long
\setlength{\abovedisplayskip}{6pt} \setlength{\belowdisplayskip}{6pt}
\title{ReLU-Based and DNN-Based \\Generalized Maximum Score Estimators\thanks{We thank Karun Adusumilli, Joel Horowitz, Simon Lee, Charles Manski, Elie Tamer, Yuanyuan Wan, and seminar participants at Columbia and Penn for helpful comments and suggestions.}}
\author{Xiaohong Chen\thanks{Chen: Department of Economics and Cowles Foundation for Research in Economics,
Yale University, 28 Hillhouse Ave, New Haven, CT 06511, USA, [email removed]}$\ $, Wayne Yuan Gao\thanks{Gao: Department of Economics, University of Pennsylvania, 133 S 36th
St., Philadelphia, PA 19104, USA, [email removed].}$\ $ and Likang Wen\thanks{Wen: Department of Applied Mathematics and Statistics, Johns Hopkins University, 3400 N Charles St,
Baltimore, MD 21218, USA, [email removed].}\textbf{}\\
\textbf{~}}
\maketitle
\begin{abstract}
\noindent We propose a new formulation of the maximum score estimator
that uses compositions of rectified linear unit (ReLU) functions, instead of indicator functions as in \citet*{manski1975maximum,manski1985semiparametric},
to encode the sign alignment restrictions. Since the ReLU function
is Lipschitz, our new ReLU-based maximum score criterion function is substantially easier to optimize using standard
gradient-based optimization pacakges. We also show that our ReLU-based maximum score (RMS) estimator can be generalized to an umbrella framework defined by multi-index single-crossing (MISC) conditions, while the original maximum score estimator cannot be applied. We establish the $n^{-s/(2s+1)}$ convergence rate and asymptotic normality for the RMS estimator under order-$s$ Holder smoothness. In addition, we propose an alternative estimator using a further reformulation of RMS as a special layer in a deep neural network (DNN) architecture, which allows the estimation procedure to be implemented via state-of-the-art software and hardware for DNN. \\
\textbf{Keywords:} semiparametric estimation, maximum score, discrete choice, rectified linear unit, deep neural network, multi-index
\end{abstract}
\section{\label{sec:Intro}Introduction}
In a sequence of papers, \citet*{manski1975maximum,manski1985semiparametric}
proposed and analyzed the properties of the \emph{maximum-score estimator}
in the context of semiparametric discrete choice models. To be specific,
consider the following canonical binary choice model
\begin{equation}
y_{i}=\mathbf{\mathbbm1}\left\{ X_{i}^{'}\theta_{0}\geq\epsilon_{i}\right\} \label{eq:BinChoice}
\end{equation}
under the conditional median restriction $\text{med}\left(\rest{\epsilon_{i}}X_{i}\right)=0$.
The key idea underlying the maximum score estimator is to exploit
the following identifying restriction,
\begin{equation}
h_{0}\left(X_{i}\right):=\mathbb{E}\left[\rest{y_{i}-\frac{1}{2}}X_{i}\right]\gtrless0\quad\Leftrightarrow\quad X_{i}^{'}\theta_{0}\gtrless0,\label{eq:Mono_Equiv}
\end{equation}
which is a \emph{sign alignment} restriction between the function
$h_{0}$ and the the parametric index $X_{i}^{'}\theta_{0}$. \citet*{manski1975maximum,manski1985semiparametric,manski1987semiparametric}
encodes this sign alignment restriction into the following population
criterion function,
\begin{align}
Q_{MS}\left(\theta\right) & :=\mathbb{E}\left[h_{0}\left(X_{i}\right)\mathbf{\mathbbm1}\left\{ X_{i}^{'}\theta>0\right\} \right]=\mathbb{E}\left[\left(y_{i}-\frac{1}{2}\right)\mathbf{\mathbbm1}\left\{ X_{i}^{'}\theta>0\right\} \right],\label{eq:Q_MS}
\end{align}
which is constructed by multiplying the function $h_{0}$ with an
indicator function of the index $X_{i}^{'}\theta_{0}$ along with the
Law of Iterated Expectation. Then, $\theta_{0}$ is a maximizer of $Q_{MS}\left(\theta\right)$.
To see more clearly why $\theta_{0}$ maximizes $Q_{MS}\left(\theta\right)$,
consider the following decomposition $h_{0}\left(X_{i}\right)\equiv\left[h_{0}\left(X_{i}\right)\right]_{+}-\left[-h_{0}\left(X_{i}\right)\right]_{+},$
where $\left[t\right]_{+}:=\max\left(t,0\right)$ denotes the rectified
linear unit (ReLU) function. Then, the population criterion $Q_{MS}$
can be correspondingly decomposed as $Q_{MS}\left(\theta\right)=Q_{MS+}\left(\theta\right)+Q_{MS-}\left(\theta\right)$
with
\begin{align*}
Q_{MS+}\left(\theta\right) & =\mathbb{E}\left[\left[h_{0}\left(X_{i}\right)\right]_{+}\mathbf{\mathbbm1}\left\{ X_{i}^{'}\theta>0\right\} \right]\leq\mathbb{E}\left[\left[h_{0}\left(X_{i}\right)\right]_{+}\right]=Q_{MS+}\left(\theta_{0}\right),\\
Q_{MS-}\left(\theta\right) & =-\mathbb{E}\left[\left[-h_{0}\left(X_{i}\right)\right]_{+}\mathbf{\mathbbm1}\left\{ X_{i}^{'}\theta>0\right\} \right]\leq0=Q_{MS-}\left(\theta_{0}\right).
\end{align*}
In words, the multiplication of $h_{0}$ with the indicator on $X_{i}^{'}\theta_{0}$
precisely extracts the positive part of $h_{0}$ at the true $\theta_{0}$,
and hence $Q_{MS}\left(\theta\right)\leq\mathbb{E}\left[\left[h_{0}\left(X_{i}\right)\right]_{+}\right]=Q_{MS}\left(\theta_{0}\right).$
In this paper, we propose a different population criterion function
that encodes exactly the same sign alignment restriction \eqref{eq:Mono_Equiv}
above. However, instead of using multiplication with indicator functions
on $X_{i}^{'}\theta_{0}$ as in \eqref{eq:Q_MS}, our new formulation
employs compositions of ReLU functions. Specifically, define
\begin{align}
g_{+,\theta,h}\left(x\right):=\left[h\left(x\right)-\left[-x^{'}\theta\right]_{+}\right]_{+},\quad & g_{-,\theta,h}\left(x\right):=\left[-h\left(x\right)-\left[x^{'}\theta\right]_{+}\right]_{+},\label{eq:def_g}
\end{align}
with $Q_{+}\left(\theta\right):=\mathbb{E}\left[g_{+,\theta,h_{0}}\left(X_{i}\right)\right]$,
$Q_{-}\left(\theta\right):=\mathbb{E}\left[g_{-,\theta,h_{0}}\left(X_{i}\right)\right]$,
and
\begin{equation}
Q\left(\theta\right):=Q_{+}\left(\theta\right)+Q_{-}\left(\theta\right),\label{eq:def_Q}
\end{equation}
Clearly, both $g_{+}$ and $g_{-}$, and thus $Q_{+}$ and $Q_{-}$,
are by construction nonnegative.
To see why $\theta_{0}$ is also a maximizer of $Q_{MS}\left(\theta\right)$,
first consider the case when $h_{0}\left(X_{i}\right)>0$. By \eqref{eq:Mono_Equiv},
\begin{align*}
h_{0}\left(X_{i}\right)>0\ & \Leftrightarrow\ X_{i}^{'}\theta_{0}>0\ \Leftrightarrow\ \left[-X_{i}^{'}\theta_{0}\right]_{+}=0\ \Rightarrow\ h_{0}\left(X_{i}\right)-\left[-X_{i}^{'}\theta_{0}\right]_{+}=\left[h_{0}\left(X_{i}\right)\right]_{+},
\end{align*}
and thus
\[
0\leq g_{+,\theta,h_{0}}\left(X_{i}\right)=\left[h_{0}\left(X_{i}\right)-\left[-X_{i}^{'}\theta\right]_{+}\right]_{+}\leq\left[h_{0}\left(X_{i}\right)\right]_{+}=g_{+,\theta_{0},h_{0}}\left(X_{i}\right).
\]
Furthermore, when $h_{0}\left(X_{i}\right)>0$, the negative part
degenerates to 0, i.e.,
\[
g_{-,\theta,h_{0}}\left(x\right)=\left[-h_{0}\left(X_{i}\right)-\left[X^{'}\theta\right]_{+}\right]_{+}\equiv0,
\]
regardless of the parameter value $\theta$. Similarly, the opposite holds
for the case of $h_{0}\left(X_{i}\right)<0$. Together, we have
\begin{align*}
Q_{+}\left(\theta\right) & =\mathbb{E}\left[\left[h_{0}\left(X_{i}\right)-\left[-X_{i}^{'}\theta\right]_{+}\right]_{+}\right]\leq\mathbb{E}\left[\left[h_{0}\left(X_{i}\right)\right]_{+}\right]=Q_{+}\left(\theta_{0}\right),\\
Q_{-}\left(\theta\right) & =\mathbb{E}\left[\left[-h_{0}\left(X_{i}\right)-\left[X_{i}^{'}\theta\right]_{+}\right]_{+}\right]\leq\mathbb{E}\left[\left[-h_{0}\left(X_{i}\right)\right]_{+}\right]=Q_{-}\left(\theta_{0}\right),
\end{align*}
which implies that $Q\left(\theta\right)\leq\mathbb{E}\left[\left|h_{0}\left(X_{i}\right)\right|\right]=Q\left(\theta_{0}\right).$
Hence, our ReLU-based criterion $Q$, even though different from the
original maximum score criterion $Q_{MS}$ above, also incorporates
the identifying restriction about $\theta_{0}$ and can thus serve as
a valid population criterion.
More generally, in a $J$-index setting we let
\[
X_i := (X_{i1},\dots,X_{iJ}) \in \mathcal{X} \subset \mathbb{R}^{d\times J},
\]
and write $x=(x_1,\dots,x_J)$ for a generic realization. For a generic function
$h:\mathcal{X}\to\mathbb{R}$ and direction $\theta\in\Theta\subset\mathbb{S}^{d-1}$,
we define
\begin{align}
g_{+,\theta,h}(x_1,\dots,x_J)
&:= \biggl[h(x_1,\dots,x_J) - \Bigl(\min_{1\le j\le J}(-x_j'\theta)_+\Bigr)\biggr]_+,
\nonumber \\
g_{-,\theta,h}(x_1,\dots,x_J)
&:= \biggl[-h(x_1,\dots,x_J) - \Bigl(\min_{1\le j\le J}(x_j'\theta)_+\Bigr)\biggr]_+.
\label{eq:def_g_J}
\end{align}
The corresponding $J$-index RMS population criterion is $Q_J(\theta) := Q_J^+(\theta)+Q_J^-(\theta)$ with
\begin{equation}
Q_J^+(\theta) := E\bigl[g_{+,\theta,h_0}(X_i)\bigr],\qquad
Q_J^-(\theta) := E\bigl[g_{-,\theta,h_0}(X_i)\bigr].
\label{eq:def_Q_J}
\end{equation}
In the single-index case $J=1$, $X_i$ reduces to a single vector $X_i\in\mathbb{R}^d$,
$g_{+,\theta,h}$ and $g_{-,\theta,h}$, and $Q_J(\theta)$ reduce to those defined in \eqref{eq:def_g} and \eqref{eq:def_Q}.
The main focus of this paper is to show how this new ReLU-based population
criterion $Q$, as defined by \eqref{eq:def_g}-\eqref{eq:def_Q_J},
can be used for the identification, estimation and inference of $\theta_{0}$, and demonstrate that this new approach relates to, differs from, and
improves upon the existing approach based on $Q_{MS}$.
~
We first focus on the binary choice
setting in Section \ref{sec:BinChoice}, which is not only a topic of important
interest on its own, but also serves as a canonical setup where our
new ReLU-based estimator can be related to the original maximum score
(MS) estimator and its previous variants in a clear manner.
Under the binary choice setting, we propose the \emph{ReLU-based maximum
score} (RMS) estimator as a semiparametric two-stage M-estimator based
on the population criterion $Q$. Specifically, in the first stage,
we obtain an estimator $\hat{h}$ of $h_{0}$ via nonparametric regression
of $Y_{i}-\frac{1}{2}$ on $X_{i}$. Then, we define the sample criterion
function $\hat{Q}$ as the sample analog of $Q$ with $\hat{h}$ plugged
in for $h_{0}$, and obtain the RMS estimator $\hat{\theta}$ as the maximizer
of the sample criterion function $\hat{Q}$ in the second stage. We
establish the convergence rate and asymptotic normality for the RMS estimator under lower-level conditions on
the primitives of the binary choice model, with $\hat{h}$ given by
kernel or linear series estimators.
In particular, we show that, under appropriate conditions, the RMS
estimator is asymptotically normal with rate of convergence as fast as $n^{-\frac{s}{2s+1}}$ (with $s$ being the imposed order of smoothness).
This rate is slower than the $\sqrt{n}$ rate but faster than the
$n^{1/3}$-rate of the original MS estimator \citep{kim1990cube},
and it coincides with the rate of the \emph{smoothed maximum score}
(SMS) estimator in \citet*{horowitz1992smoothed}. The RMS and SMS
estimators are conceptually similar in the sense that both exploit
additional smoothness conditions (on $h_{0}$, in particular) relative
to the original MS estimator, which leads to the accelerated convergence
rates.\footnote{Recall also from \citet*{horowitz1992smoothed} that the rate $n^{-\frac{s}{2s+1}}$
cannot be further improved upon in the minimax sense.} However, the asymptotic theory of the RMS estimator differs significantly
from that for the SMS estimator given the very different forms of
population and sample criterion functions involved.
In particular, the intermediate level of (non-)smoothness in the ReLU
function turns out to be a key driver of the asymptotic behavior of
the RMS estimator. First, the ``kink'' of the ReLU function at 0
(or more precisely, a non-zero first-order derivative from one side)
is essential for the locally quadratic curvature of the population
criterion function around the true parameter $\theta_{0}$. Second, the
Lipschitz continuity of ReLU functions, in contrast with the discontinuous
indicator function, translates small deviations into small deviations,
which is key for a stochastic equicontinuity condition that reduces
the impact of the first-stage nonparametric estimation errors on the
second stage and helps with the convergence rate as well as the asymptotic
normality (instead of a Chernoff-type asymptotic distribution). Third,
the almost-everywhere differentiability of the ReLU function enables
the characterization of the leading term in the asymptotic analysis
as a plug-in estimator of an integration functional of the nonparametric
function $h_{0}\left(x\right)$ over a $\left(d-1\right)$-dimensional
hyperplane (with $d$ being the dimension of $X_{i}$, i.e., the dimension
of the first-stage nonparametric estimation of $h_{0}$). This integral
averages the first-stage estimation error in $\hat{h}$ over a $\left(d-1\right)$-dimensional
space, thus accelerating the convergence to the rate of 1-dimensional
nonparametric estimation, which is the fundamental driver of the final
$n^{-\frac{s}{2s+1}}$ rate of the RMS estimator.\footnote{Relatedly, the asymptotic theory of the SMS estimator \citep{horowitz1992smoothed}
is also driven by the convergence rate of 1-dimensional nonparametric
(kernel) estimation.}
We then (in Section~3) generalize the RMS estimator to an umbrella econometric framework characterized by multi-index single-crossing (MISC) conditions proposed in \cite{gao2020robust}. We show that MISC conditions arise naturally in a wide range of econometric models, and are particularly powerful in multi-index discrete choice and panel multinomial choice settings. In particular, the MISC framework underlies the identification and estimation strategy in \cite{gao2020robust} and \cite*{gao2023logical}, where multi-index single-crossing restrictions are exploited to obtain semiparametric identification in panel multinomial choice models. Our analysis provides a complementary perspective by showing how ReLU-based maximum score ideas can be embedded in the MISC framework and extended to a broad class of models beyond the binary choice benchmark.
Beyond the traditional two-step semiparametric implementation, we also show in Section~\ref{sec:NN} how the RMS/MISC framework can be embedded in a multi-layer neural network architecture. In particular, we construct a special ``RMS layer'' that takes as input a flexible first-stage network $h(x)$ and a low-dimensional direction $\theta$, and applies the composite ReLU transformation that encodes the sign-alignment or MISC restriction. This provides a concrete example of how
economically meaningful low-dimensional parameters can be built into (and estimated within) deep neural networks (DNN) using standard machine learning toolkits. In this way, the paper speaks directly to the broader literatures on interpretable deep learning, by demonstrating how modern neural networks can be used to capture rich nonparametric structure without sacrificing identification for the structural index parameter.
~
Our paper contributes directly to the econometric literature on maximum
score (MS) estimators, dating back to \citet{manski1975maximum,manski1985semiparametric},
and \citet*{kim1990cube}. Of particular relevance is the line of
research on the variants of the MS estimator with different forms
of smoothing. To our best knowledge, our paper is the first to propose
the ReLU-based formulation introduced above, which builds an intermediate
level of smoothness directly into the population criterion. Previously,
\citet*{horowitz1992smoothed} proposes the SMS estimator, where the
indicator function in the MS (sample) criterion is replaced by a smooth
sigmoid function with a bandwidth parameter, and establishes the accelerated
convergence rate and asymptotic normality of the SMS estimator. \citet*{blevins2013local}
works with a local nonlinear least square formulation of the SMS estimator,
and uses debiasing techinques to obtain the SMS convergence rate. \citet*{chen2015binary}
reformulates the sign alignment restriction as a local conditional
moment condition and proposes a corresponding estimator based on local
polynomial smoothing. \citet*{jun2017integrated} considers the integrated
score estimator, a quasi-Bayes estimator where smoothing is achieved
through integration of the MS criterion. Another set of related work
focuses on the inference problem, given that standard bootstrap is
known to be invalid for the MS estimator \citep*{abrevaya2005bootstrap}:
\citet*{horowitz2002bootstrap} establishes bootstrap consistency for
the SMS estimator, \citet*{patra2018consistent} formulates a smoothed
bootstrap procedure for the MS estimator using a semiparametric two-stage
estimator to center the bootstrap samples,\footnote{This semiparametric two-stage estimator in \citet*{patra2018consistent},
defined in their equation (5), utilizes a first-stage nonparametric
estimation of $h_{0}$, which is plugged in along with a nonparametric
density estimator to obtain an integrated estimator of the MS population
criterion function. However, this estimator is then used for the bootstrap
of the original MS estimator, and its properties were not fully developed
in \citet*{patra2018consistent}. } while \citet*{cattaneo2020bootstrap} proposes an alternative approach
to obtain bootstrap consistency by modifying an asymptotically non-random
component of the MS sample criterion. None of the papers cited above
considers our ReLU-based formulation. As discussed above, this new formulation not only leads to a ``more
smooth'' population criterion function that provides both theoretical
and computational advantages, but also greatly generalizes the scope
of applications to which the key idea of maximum score estimation
can be applied.
This paper also builds upon and contributes to the long line of econometric
literature on semiparametric M estimation and inference: see, for
example, \citet*{newey1994large}, \citet*{chen2007sieve}, \citet*{ichimura2007implementing},
and \citet*{kosorok2007introduction} for general surveys on this
topic. In particular, this paper is related to previous work that
analyzes nonsmooth criterion functions, such as \citet*{chen2003estimation},
\citet*{ichimura2010characterization,ichimura2018corrigendum}, \citet*{seo2018local},
and \citet*{delsol2020semiparametric}. A distinct feature of this
paper is the intermediate level of smoothness (``Lipschitz with a
kink'') of the ReLU function leads to the intermediate convergence
rate of the RMS estimator, which is faster than the cubic-root-or-slower
rates obtained in \citet*{kim1990cube}, \citet*{seo2018local} and
the example considered in \citet*{delsol2020semiparametric} (with
``less smooth'' criterion functions), but slower than the root-$n$
rate considered by \citet*{chen2003estimation} and \citet*{ichimura2010characterization,ichimura2018corrigendum}
(with ``more smooth'' criterion functions). More specifically, we show how the ``Lipschitz-with-a-kink'' property of the ReLU function leads to a characterization of the leading term in the RMS asymptotics as a nonparametric plug-in estimator of a lower-dimensional integral functional, and how this lower-dimensional integral becomes the key driver of the final intermediate convergence rate. Our results on the convergence of nonparametric integral functionals over lower-dimensional hyperplanes are of independent interest, which is closely related to the general theory of semiparametric learning of integral functionals on submanifolds developed in \cite{chen2025semiparametric}, which explicitly relates the convergence rate to the dimension of the underlying submanifold. Our contribution also supplements related work in the statistics literature on the estimation of integrals on level sets, which mostly focus on kernel regressions \citep*{dau2020exact} or density estimation \citep*{qiao2021nonparametric}.
Our DNN-based maximum score estimator under the MISC condition framework also speaks directly to the broader machine learning literature on interpretability of deep neural networks (DNN). Surveys such as \cite{fan2021interpretability}
and \cite{zhang2021survey} review a wide range of interpretability tools, which mostly focus on explaining \emph{predictions} or internal
representations, but not on identifying or conducting inference on structural low-dimensional parameters inside a network. In this sense, the DNN-based MISC estimator offer a way to bridge the gap between the interpretability and uncertainty literatures in deep learning and the semiparametric inference literature in econometrics. They allow researchers to use modern DNN to capture rich nonlinearities and heterogeneity in
the data, while still retaining (i) an interpretable, low-dimensional parameter $\theta$ that encodes economically meaningful structure, and (ii) a rigorous large-sample theory that supports conventional confidence intervals and hypothesis tests for that parameter.
~
The rest of the paper is organized as follows. Section~\ref{sec:BinChoice}
introduces the RMS estimator in the binary choice model, develops the basic
identification and asymptotic theory, and compares RMS to the original and
smoothed maximum score estimators. Section~\ref{sec:MISC} embeds the binary
choice setup into the general multi-index single-crossing framework, extends
the RMS criterion to the $J$-index case, and derives the corresponding
asymptotic results, highlighting the effective one-dimensional nature of the
rate. Section \ref{sec:NN} further reformulates the RMS as a specialized layer in a DNN, which allows the estimation of the index parameter to be subsumed under the training of the DNN, for which state-of-art computing software on DNN become applicable. Section~\ref{sec:Sim} presents simulation evidence on the finite-sample
performance of the RMS estimator in both single-index and multi-index designs.
Section~\ref{sec:Con} concludes. Technical proofs and additional auxiliary
results are collected in the appendix.
\section{\label{sec:BinChoice}Special Case: Binary Choice Model}
In this section, we focus on the binary choice model \eqref{eq:BinChoice}
as described in the introduction, and develops the econometric theory
of our ReLU-based maximum score (RMS) estimator with clear lower-level
conditions on the primitives of the model. The binary choice model
is not only of important interest on its own, but also serves as a
canonical setup where our new ReLU-based estimator can be related
to the original maximum score (MS) estimator and its previous variants
in a clear manner.
\subsection{\label{subsec:SetupResults}Setup and Main Results}
Given the binary choice model \eqref{eq:BinChoice} and the ReLU-based
population criterion function $Q$ in \eqref{eq:def_Q}, we define
the ReLU-based maximum score (RMS) estimator as
\begin{equation}
\hat{\theta}:=\arg\max_{\theta\in\mathbb{\mathbb{S}}^{d-1}}\hat{Q}\left(\theta\right)\label{eq:RMS}
\end{equation}
where the sample criterion function $\hat{Q}$ is given by
\[
\hat{Q}\left(\theta\right):=\frac{1}{n}\sum_{i=1}^{n}\left(g_{+,\theta,\hat{h}}\left(X_{i}\right)+g_{-,\theta,\hat{h}}\left(X_{i}\right)\right)
\]
with $\hat{h}$ being some first-stage nonparametric estimator of
$h_{0}\left(x\right)=\mathbb{E}\left[\rest{y_{i}-\frac{1}{2}}X_{i}=x\right]$.
We seek to characterize the asymptotic behaviors of the RMS estimator
$\hat{\theta}$.
It turns out that the $\hat{\theta}$ is ``non-standard'' semiparametric
two-stage M-estimator, and is different from both the usual ``$\sqrt{n}$-normal''
asymptotics in the ``smooth case'' (such as in \citealp{newey1994asymptotic})
and the ``cubic-rate'' asymptotics in \citet*{kim1990cube}.
As we will show subsequently, the ReLU-based maximum score estimator
will feature ``intermediate'' asymptotics (under appropriate conditions
to be made explicit later): $\hat{\theta}$ will converge at nonparametric
rates slower than $n^{\frac{1}{2}}$ but faster than $n^{\frac{1}{3}}$
with asymptotic normal distribution, which can be viewed as a ``semiparametric
two-stage version'' of the asymptotic results in \citealp*{horowitz1992smoothed}.
In particular, the ``intermediate asymptotics'' of ReLU-based maximum
score estimator $\hat{\theta}$ is critically driven by the ``intermediate
smoothness'' allowed by the formulation of the criterion function
\eqref{eq:def_Q} using the the ReLU function $\left[\cdot\right]_{+}$,
which is Lipschitz continuous and everywhere differentiable \emph{except}
at the single ``kink point'' $0$. Interestingly, both the ``smoothness''
and ``kinkiness'' of the ReLU function turns out to be important:
while the Lipschitz continuity of the ReLU function is key in delivering
a ``stochastic equicontinuity'' condition for asymptotic normality,
and the ``kinkiness'' of the ReLU function at $0$ is key in delivering
locally quadratic identification of $\theta_{0}$, i.e., the quadratic
curvature of the population criterion function $Q$ in a neighborhood
of $\theta_{0}$.
~
We start by imposing a set of lower-level assumptions that guarantees
the point identification of $\theta_{0}$ (under scale normalization)
by the RMS criterion function \eqref{eq:def_Q} and that a variety
of densities are smooth and well-behaved. We note that these assumptions
are stronger than necessary, but lend simplicity to the exposition
of our main results.
\begin{assumption}
\label{assu:Basic} Write ${\cal X}:=\text{Supp}\left(X_{i}\right)\subseteq\mathbb{R}^{d}$.
Suppose $\theta_{0}\in\mathbb{\mathbb{S}}^{d-1}$ and the following:
\begin{itemize}
\item[(a)] $\left(y_{i},X_{i},\epsilon_{i}\right)_{i=1}^{n}$ is i.i.d. and satisfies
model \eqref{eq:BinChoice}.
\item[(b)] The conditional median of $\epsilon_{i}$ given $X_{i}=x$ is zero, i.e.,
\[
F\left(\rest 0x\right)=\frac{1}{2},\quad\forall x\in{\cal X}.
\]
\item[(c)] The (unknown) conditional CDF $F\left(\rest{\epsilon}x\right)$ of $\epsilon_{i}$
given $X_{i}=x$ is $d$ times continuously differentiable w.r.t.
$\left(\epsilon,x\right)\in\mathbb{R}\times{\cal X}$ with uniformly bounded derivatives
(bounded by some positive constant $M<\infty$).
\item[(d)] The conditional probability density function $f\left(\rest{\epsilon}x\right)$
of $\epsilon_{i}$ given $X_{i}=x$ is strictly positive for any $\epsilon\in\mathbb{R}$
and $x\in{\cal X}$.
\item[(e)] Furthermore, there exists a finite $M>0$ such that
\[
0<\frac{1}{M}\leq f\left(\rest 0x\right)\leq M,\quad\text{for all }x\in{\cal X}.
\]
\item[(f)] ${\cal X}$ is compact in $\mathbb{R}^{d}$ and contains ${\bf 0}$ as an
interior point. WLOG assume $\norm{x}\leq 1, \forall x\in {\cal X}$.
\item[(g)] Let $p\left(x\right)$ be the probability density function of $X_{i}$.
There exists a finite $M>0$ such that
\[
0<\frac{1}{M}\leq p\left(x\right)\leq M,\quad\text{for all }x\in{\cal X}.
\]
\end{itemize}
\end{assumption}
Assumption \ref{assu:Basic}(a) and (b) consists of a standard random-sampling
assumption for the binary choice model \eqref{eq:BinChoice} with
a conditional median restriction, which are essentially the same as
those imposed in \citet*{horowitz1992smoothed}. Note, however, we
focus on the binary choice model here as a key illustration, but,
just as maximum score estimator can be applied to many models other
than the binary choice model \eqref{eq:BinChoice}, our proposed method
can also be adapted to other settings. See XXX for a more detailed
discussion.
Assumption \ref{assu:Basic}(c)-(e) are regularity conditions on the
conditional distribution of the error term $\epsilon_{i}$, which correspond
to Assumptions 2(b), 9 and 11 in \citet*{horowitz1992smoothed}. The
assumptions of the existence and boundedness (from above) of conditional
densities and their derivatives impose smoothness conditions on model
\eqref{eq:BinChoice} and the conditional expectation function $h_{0}\left(x\right)$
beyond \citet{manski1985semiparametric} and \citet{kim1990cube}.
As in \citet*{horowitz1992smoothed}, these smoothness conditions
are exploited to deliver faster convergence rates than the cubic rate
as well as asymptotic normality. Note that the ``bounded away from
zero'' assumption $f\left(\rest 0x\right)>\frac{1}{M}$ is a local-identification
assumption that deliver the quadratic curvature in the population
criterion function, which is imposed implicitly in Assumption 11 of
\citet*{horowitz1992smoothed}.
Assumption \ref{assu:Basic}(f) imposes assumption on ${\cal X}$,
the support of the covariates $X_{i}$. In particular, the assumption
of ${\cal X}$ containing ${\bf 0}$ as an interior point guarantees
that $X_{i}$ has full ``directional'' support, i.e. $X_{i}/\norm{X_{i}}$
is supported on the whole $\mathbb{\mathbb{S}}^{d-1}$. As explained in \citet*{manski1985semiparametric},
the identification of $\theta_{0}$ is driven by variations in the ``directions''
$X_{i}/\norm{X_{i}}$, and the full-directional-support condition
ensures that $\theta_{0}$ is point identified on $\mathbb{\mathbb{S}}^{d-1}$. As well-known
in the literature, the assumption of ${\bf 0}$ being in the interior
of ${\cal X}$ is a sufficient, but not necessary, condition for the point identification
of $\theta_{0}$. Alternatively, one could work with a ``special regressor''
as in Assumptions 2(a)(c) \& 4 in \citet*{horowitz1992smoothed},
which assume that $\left|\beta_{01}\right|=1$ and that the conditional
distribution of $X_{i1}$ given $\left(X_{i2},...,X_{id}\right)$
has full support on $\mathbb{R}$. This alternative set of assumptions allows
for discreteness in certain components of $X_{i}$ but rules out compactness
of ${\cal X}$, and thus do not nest Assumptions \ref{assu:Basic}(f)(g)
as special cases, nor vice versa. Furthermore, the scale normalization
$\left|\beta_{01}\right|=1$ is dependent on the assumption that a specific
known component of $X_{i}$ has non-zero coefficient. In this paper,
we focus on the normalization $\beta_{0}\in\mathbb{R}^{d}$, i.e., $\norm{\beta_{0}}=1$
and the support Assumptions \ref{assu:Basic}(f)(g), which lends simpler
notation in our asymptotics. However, the substance of our asymptotic
results is not dependent on this specific choice of point-identifying
assumption and scale normalization, and it should be feasible, though
notationally cumbersome, to adapt our asymptotic results to the set
of assumption and normalization using the ``special regressor''
as in \citet*{horowitz1992smoothed}.
We also, note that the assumption of compactness of ${\cal X}$ in
Assumption \ref{assu:Basic}(f) is not necessary, either. Compactness
of ${\cal X}$ is often assumed in the literature, and assumed here
for simpler exposition of results on the nonparametric estimation
of $h_{0}\left(x\right)$. Hence, our results based on the consistency
and convergence rate of nonparametric estimation of $h_{0}$ on compact
${\cal X}$ can be adapted to the case where ${\cal X}$ is not compact
with standard trimming and/or weighting of ${\cal X}$.
Lastly, Assumption \ref{assu:Basic}(g) corresponds to Assumptions
8 and 11 in \citet*{horowitz1992smoothed}, imposing both smoothness
(in terms of bounded-from-above densities) and local-identification
conditions (in terms of bounded-away-from zero densities). Again,
Assumption \ref{assu:Basic}(g) is stated in a stronger-than-necessary
but expositionally simple form. In particular, for local identifiability
it is not necessary to require that $p\left(x\right)$ is bounded
away from zero at every point in ${\cal X}$, since local identifiability
is only concerned with the hyperplane $\left\{ x:x^{'}\theta_{0}=0\right\} $,
not the entire ${\cal X}$. However, given imposed compactness of
${\cal X}$ in Assumption \ref{assu:Basic}(f), the global ``bounded-away-from-zero''
condition here is not very restrictive anyway, and hence we impose
this stronger-than-necessary condition for simpler notation.
We summarize two important implications of Assumption \ref{assu:Basic}
below:
\begin{prop}
\label{prop:ID_Sobolev}Under Assumption \ref{assu:Basic}:
\begin{itemize}
\item[(i)] $\theta_{0}$ is point identified on $\mathbb{\mathbb{S}}^{d-1}$:
\begin{equation}
\theta_{0}=\arg\max_{\theta\in\mathbb{\mathbb{S}}^{d-1}}Q\left(\theta\right).\label{eq:PointID}
\end{equation}
\item[(ii)] $h_{0}\left(x\right)$ is $\left(d+1\right)$ times differentiable
on ${\cal X}$ with uniformly bounded derivatives.
\end{itemize}
\end{prop}
~
Given mild convergence conditions on the first-stage estimator $\hat{h}$,
it is straightforward to establish the consistency of $\hat{\theta}$
in Theorem \ref{thm:Consistency}.
\begin{thm}[Consistency]
\label{thm:Consistency}Suppose that $\norm{\hat{h}-h_{0}}_{\infty}=o_{p}\left(1\right)$.
Then $\hat{\theta}$ is consistent, i.e., $\hat{\theta}\overset{p}{\longrightarrow}\theta_{0}$.
\end{thm}
~
We now proceed to characterize the convergence rate and asymptotic
distribution of $\hat{\theta}$, which are the main results of this section.
While such results can be obtained under higher-level conditions on
the first-stage nonparametric estimators $\hat{h}$, for concreteness
and clarity, we consider two leading types of nonparametric estimators,
the Nadaraya-Waston kernel estimator and the linear series estimator,
and provide lower-level conditions for both.
\begin{assumption}[Kernel/Linear Series First Stage]
\label{assu:KernelSeries} Assume either of the following:
\begin{itemize}
\item[\emph{(a)}] \begin{flushleft}
$\hat{h}$ is given by the Nadaraya-Watson kernel estimator,
\[
\hat{h}\left(x\right):=\frac{\sum_{i=1}^{n}K\left(\frac{X_{i}-x}{b_{n}}\right)\left(Y_{i}-\frac{1}{2}\right)}{\sum_{i=1}^{n}K\left(\frac{X_{i}-x}{b_{n}}\right)}
\]
where $b_{n}$ is a bandwidth parameter and $K\left(x\right)$ is
a $d$-dimensional kernel function of smoothness order $s$ such that
\par\end{flushleft}
\item[] (a.i) $K\left(x\right)=K\left(-x\right)$, and $\int K\left(x\right)dx=1$.
\item[] (a.ii) $\left|K\left(x\right)\right|\leq M<\infty$ and $\int\prod_{j=1}^{d}\left|x_{j}\right|^{l}K\left(x\right)dx<\infty$
for all $l$.
\item[] (a.iii) $\int k\left(t\right)=0$ for $j=1,...,s-1$, and $\kappa_{s}:=\int x_{j}^{s}K\left(x\right)dx\in\left(0,\infty\right)$.
\item[] (a.iv) $K\left(x_{1},...,x_{d}\right)=K\left(x_{\pi_{1}},...,x_{\pi_{d}}\right)$
for any permutation of coordinates $\pi$.
\item[\emph{(b)}] \begin{flushleft}
\emph{ $\hat{h}$ is given by the linear series estimator,}
\[
\hat{h}\left(x\right):=\ol b^{K_{n}}\left(x\right)^{'}\left(\sum_{i=1}^{n}\ol b^{K_{n}}\left(X_{i}\right)\ol b^{K_{n}}\left(X_{i}\right)^{'}\right)^{-1}\sum_{i=1}^{n}\ol b^{K_{n}}\left(X_{i}\right)\left(Y_{i}-\frac{1}{2}\right)
\]
where $K_{n}:=J_{n}^{d}$ is the sieve dimension parameter and $\ol b^{K_{n}}\left(x\right):=\text{vec}\left(\otimes_{j=1}^{d}\left(b_{1}\left(x_{j}\right),...,b_{J_{n}}\left(x_{j}\right)\right)\right)$
is a vector of multivariate basis functions constructed from tensor
products of some univariate orthonormal basis functions $\left(b_{k}\left(\cdot\right)\right)_{k=1}^{\infty}$
such that:
\par\end{flushleft}
\item[] (b.i)\emph{ $\lambda_{\min}\left(\mathbb{E}\left[\ol b^{K_{n}}\left(X_{i}\right)\ol b^{K_{n}}\left(X_{i}\right)^{'}\right]\right)>0$.}
\item[] (b.ii)\emph{ $\inf_{h\in{\cal B}_{K_{n}}}\norm{h-h_{0}}=J_{n}^{-s}$,
where }${\cal B}_{K_{n}}$ denotes the closed span of $\ol b^{K_{n}}$.
\item[] (b.iii)\emph{ $\norm{\Pi_{K_{n},n}}_{\infty}:=\sup_{h:\norm h_{\infty}\neq0}\frac{\norm{\Pi_{K_{n},n}h}_{\infty}}{\norm h_{\infty}}=O_{p}\left(1\right)$,
where} $\Pi_{K_{n},n}$ denotes the empirical projection operator
onto ${\cal B}_{K_{n}}$, i.e.,
\[
\Pi_{K_{n},n}h\left(x\right):=\ol b^{K_{n}}\left(x\right)^{'}\left(\sum_{i=1}^{n}\ol b^{K_{n}}\left(X_{i}\right)\ol b^{K_{n}}\left(X_{i}\right)^{'}\right)^{-1}\sum_{i=1}^{n}\ol b^{K_{n}}\left(X_{i}\right)h\left(x\right).
\]
\end{itemize}
\end{assumption}
\noindent \begin{flushleft}
Assumption \ref{assu:KernelSeries}(a.i-iv) are standard conditions
on the kernel function that covers both product and radial kernels
constructed under a wide range of univariate kernels. Similarly, Assumption
\ref{assu:KernelSeries} (b.i-iii) are standard conditions and properties
on linear series regressions that are satisfied under a wide variety
of sieve classes. See, for example, \citet{chen2007sieve}, \citet{chen2015optimal}
and \citet{belloni2015some} for results on spline, wavelet, Fourier
and many other sieve classes.
\par\end{flushleft}
We now present our main results about the RMS asymptotics.
\begin{thm}[Convergence Rate]
\label{thm:Thm_Bin_Rate} Under Assumption \eqref{assu:Basic}, and
with $\hat{h}$ being given by the Nadaraya-Watson estimator that
satisfies Assumption \emph{\ref{assu:KernelSeries}(a),} for any $b_{n}\to0$
such that $nb_{n}^{2d+1}/\left(\log n\right)^{2}\to\infty$, we have
\begin{equation}
\norm{\hat{\theta}-\theta_{0}}=O_{p}\left(b_{n}^{s}+\frac{1}{\sqrt{nb_{n}}}\right).\label{eq:Bin_Rate_b}
\end{equation}
If $s>d$, then optimal convergence rate can be attained by setting
$b_{n}\sim n^{-\frac{1}{2s+1}}$, giving
\[
\norm{\hat{\theta}-\theta_{0}}=O_{p}\left(n^{-\frac{s}{2s+1}}\right).
\]
The above also holds for linear series $\hat{h}$ with Assumption
\emph{\ref{assu:KernelSeries}(a)} replaced by Assumption \emph{\ref{assu:KernelSeries}(b)}
and $b_{n}$ replaced by $J_{n}^{-1}$.
\end{thm}
The asymptotic distribution can then be derived based on the linearized
argmax theorem (Theorem 3.2.16) in \citet*{van1996weak}.
\begin{thm}[Asymptotic Normality]
\label{thm:AsymNorm} Suppose that Assumption holds with $s>d$.
With $\hat{h}$ being given by the Nadaraya-Watson estimator as in
Assumption \emph{\ref{assu:KernelSeries}(a) with undersmoothing choice
of bandwidth} $b_{n}$ such that $nb_{n}^{2d+1}/\left(\log n\right)^{2}\to\infty$
and $b_{n}=o_{p}\left(n^{-\frac{1}{2s+1}}\right)$, we have
\[
n^{\frac{s}{2s+1}}\left(\hat{\theta}-\theta_{0}\right)\overset{d}{\longrightarrow}\mathcal{N}\left({\bf 0},V^{-}\Omega V^{-}\right).
\]
The above also holds for linear series $\hat{h}$ with Assumption
\emph{\ref{assu:KernelSeries}(a) }replaced by Assumption \emph{\ref{assu:KernelSeries}(b)}
and $b_{n}$ replaced by $J_{n}^{-1}$.
\end{thm}
\subsection{\label{subsec:AsymOutline}Outline of the RMS Asymptotic Theory}
\subsubsection{Decomposition of the Sample Criterion}
To present our formal asymptotic results, we first set up some notation.
Let $Pg_{\theta,h}:=\int g_{\theta,h}\left(x\right)dP\left(x\right),$ $\mathbb{P}_{n}g_{\theta,h}:=\frac{1}{n}\sum_{i=1}^{n}g_{\theta,h}\left(X_{i}\right),$
and $\mathbb{G}_{n}g_{\theta,h}:=\sqrt{n}\left(\mathbb{P}_{n}g_{\theta,h}-Pg_{\theta,h}\right)$,
with which we can rewrite \eqref{eq:PointID} and \eqref{eq:RMS}
as
\[
\theta_{0}=\arg\max_{\theta\in\mathbb{\mathbb{S}}^{d-1}}Pg_{\theta,h_{0}},\quad\hat{\theta}:=\arg\max_{\theta\in\mathbb{\mathbb{S}}^{d-1}}\mathbb{P}_{n}g_{\theta,\hat{h}}
\]
Since the asymptotic behavior of $\hat{\theta}$ is driven by the asymptotic
behavior of $\mathbb{P}_{n}\left(g_{\hat{\theta},\hat{h}}-g_{\theta_{0},\hat{h}}\right)$,
we analyze it by working with the following decomposition
\begin{align}
\mathbb{P}_{n}\left(g_{\hat{\theta},\hat{h}}-g_{\theta_{0},\hat{h}}\right)= & \underset{T_{1}}{\underbrace{\frac{1}{\sqrt{n}}\mathbb{G}_{n}\left(g_{\hat{\theta},h_{0}}-g_{\theta_{0},h_{0}}\right)}}+\underset{T_{2}}{\underbrace{\frac{1}{\sqrt{n}}\mathbb{G}_{n}\left(g_{\hat{\theta},\hat{h}}-g_{\theta_{0},\hat{h}}-g_{\hat{\theta},h_{0}}+g_{\theta_{0},h_{0}}\right)}}\nonumber \\
& +\underset{T_{3}}{\underbrace{P\left(g_{\hat{\theta},h_{0}}-g_{\theta_{0},h_{0}}\right)}}+\underset{T_{4}}{\underbrace{P\left(g_{\hat{\theta},\hat{h}}-g_{\theta_{0},\hat{h}}-g_{\hat{\theta},h_{0}}+g_{\theta_{0},h_{0}}\right)}}\label{eq:EP_decom}
\end{align}
and studying each of the four terms $T_{1},T_{2},T_{3}$ and $T_{4}$.
It turns out that each of the four terms is somewhat ``nonstandard''
relative to the usual case of semiparametric two-stage estimation
theory that delivers $\sqrt{n}$ asymptotic normality under standard
smoothness conditions. Furthermore, the analysis of the four terms
$T_{1},T_{2},T_{3},T_{4}$ reveals some of the key insights in the
asymptotics of our proposed ReLU-based maximum score estimator $\hat{\theta}$.
Hence, we provide an explicit account of the four terms below, where
we show that the terms $T_{1}$ and $T_{2}$ will be of smaller stochastic
orders than $\norm{\hat{\theta}-\theta_{0}}^{2}$ and thus become asymptotically
negligible, while terms $T_{3}$ and $T_{4}$ will be the asymptotically
leading terms of the order $\norm{\hat{\theta}-\theta_{0}}^{2}$. We then
combine the results about the four terms to establish the convergence
rate and asymptotic normality.
\subsubsection{\label{subsec:Term1}Analysis of Term $T_{1}=\frac{1}{\sqrt{n}}\protect\mathbb{G}_{n}\left(g_{\hat{\protect\theta},h_{0}}-g_{\protect\theta_{0},h_{0}}\right)$}
We start with term $T_{1}$, which captures the stochastic variation,
or loosely ``variance'', in the sample criterion function $\mathbb{P}_{n}g_{\hat{\theta},h_{0}}$
when the nonparametric first stage is set to the the true function
$h_{0}$. Lemma \ref{lem:Term1} below presents a maximal inequality
about $T_{1}$ with respect to $\theta$ in a small neighborhood of $\theta_{0}$:
\begin{lem}
\label{lem:Term1} For some constant $M>0$,
\begin{equation}
P\sup_{\norm{\theta-\theta_{0}}\leq\delta}\left|\mathbb{G}_{n}\left(g_{\theta,h_{0}}-g_{\theta_{0},h_{0}}\right)\right|\leq M\delta^{\frac{3}{2}}.\label{eq:Max1}
\end{equation}
\end{lem}
Loosely speaking, the result above in Lemma \ref{lem:Term1} translates
to the following stochastic bounds on $T_{1}$:
\[
T_{1}=O_{p}\left(\frac{1}{\sqrt{n}}\norm{\hat{\theta}-\theta_{0}}^{\frac{3}{2}}\right),
\]
which is $o_{p}\left(\norm{\hat{\theta}-\theta_{0}}^{2}\right)$ since $\norm{\hat{\theta}-\theta_{0}}$
converges no faster than $\frac{1}{\sqrt{n}}$ rate to zero. This
would imply that $T_{1}$ will become asymptotically negligible, which
is ``nonstandard'' in the literature.
Technically, the asymptotic negligibility of $T_{1}$ is directly
driven by the $\delta^{\frac{3}{2}}$-rate bound on the right hand side
of \eqref{eq:Max1}. To see why $\delta^{\frac{3}{2}}$ arises, notice
that
\[
\left|g_{\theta,h_{0}}\left(x\right)-g_{\theta_{0},h_{0}}\left(x\right)\right|=\left|g_{+,\theta,h_{0}}\left(x\right)-g_{+,\theta,h_{0}}\left(x\right)\right|+\left|g_{-,\theta,h_{0}}\left(x\right)-g_{-,\theta,h_{0}}\left(x\right)\right|
\]
and thus, for any $\theta$ close to $\theta_{0}$ in the sense of $\norm{\theta-\theta_{0}}\leq\delta$,
we have
\begin{align}
\left|g_{+,\theta,h_{0}}\left(x\right)-g_{+,\theta_{0},h_{0}}\left(x\right)\right| & =\left|\left[h_{0}\left(x\right)-\left[-x^{'}\theta\right]_{+}\right]_{+}-\left[h_{0}\left(x\right)\right]_{+}\right|\nonumber \\
& \leq\mathbf{\mathbbm1}\left\{ h_{0}\left(x\right)>0\right\} \cdot\mathbf{\mathbbm1}\left\{ x^{'}\theta<0\right\} \cdot\left|x^{'}\theta\right|\nonumber \\
& =\mathbf{\mathbbm1}\left\{ x^{'}\theta_{0}>0>x^{'}\theta\right\} \cdot\left|x^{'}\theta\right|\nonumber \\
& =\mathbf{\mathbbm1}\left\{ x^{'}\theta_{0}>0>x^{'}\theta_{0}+x^{'}\left(\theta-\theta_{0}\right)\right\} \cdot\left|x^{'}\theta_{0}+x^{'}\left(\theta-\theta_{0}\right)\right|\nonumber \\
& \leq\mathbf{\mathbbm1}\left\{ 0<x^{'}\theta_{0}<-x^{'}\left(\theta-\theta_{0}\right)\right\} \cdot\left(\left|x^{'}\theta_{0}\right|+\left|x^{'}\left(\theta-\theta_{0}\right)\right|\right)\nonumber \\
& \leq\mathbf{\mathbbm1}\left\{ 0<x^{'}\theta_{0}<M\norm x\delta\right\} \cdot2M\norm x\delta\label{eq:Bd_gtheta}
\end{align}
In words, the derivation above exploits the observation that $g_{+,\theta,h_{0}}\left(x\right)-g_{+,\theta_{0},h_{0}}\left(x\right)$
is nonzero only if $x^{'}\theta_{0}$ and $x^{'}\theta$ lie on different
sides of $0$, which, given the restriction $\norm{\theta-\theta_{0}}\leq\delta$,
implies that both $\left|x^{'}\left(\theta-\theta_{0}\right)\right|$ and
$x^{'}\theta_{0}$ must be bounded by $M\norm x\delta$. As a result, the
magnitude of $\left|g_{+,\theta,h_{0}}\left(x\right)-g_{+,\theta_{0},h_{0}}\left(x\right)\right|$,
which is at most $\left|x^{'}\theta\right|$, is also bounded above by
a term linear in $\delta$. Furthermore, since $\norm x$ is bounded by
the compactness of ${\cal X}$,\footnote{The compactness of ${\cal X}$ and the boundedness of $\norm x$ allow
for simpler exposition here but are not necessary. If $\norm{X_{i}}$
has unbounded support, the result in Lemma \ref{lem:Term1} will continue
to hold under mild tail-decay condition, or finite--fourth-moment
condition, on $\norm{X_{i}}$.} we have
\[
\left|g_{+,\theta,h_{0}}\left(x\right)-g_{+,\theta_{0},h_{0}}\left(x\right)\right|\leq\ol g_{\delta}\left(x\right):=\mathbf{\mathbbm1}\left\{ \left|x^{'}\theta_{0}\right|<M\norm x\delta\right\} \cdot2M\delta,
\]
and similarly for $\left|g_{-,\theta,h_{0}}\left(x\right)-g_{-,\theta,h_{0}}\left(x\right)\right|$.
Hence, $\ol g_{\delta}\left(x\right)$ is a so-called ``envelope function''
for the function class $\left\{ g_{\theta,h_{0}}\left(x\right)-g_{\theta_{0},h_{0}}\left(x\right):\norm{\theta-\theta_{0}}\leq\delta\right\} $
in the sense of
\[
\sup_{\theta:\norm{\theta-\theta_{0}}\leq\delta}\left|g_{\theta,h_{0}}\left(x\right)-g_{\theta_{0},h_{0}}\left(x\right)\right|\leq\ol g_{\delta}\left(x\right),\quad\forall x\in{\cal X}.
\]
By standard empirical process theory, such as in \citet*{van1996weak},
the magnitude of $\sqrt{\mathbb{E}\left[\ol g_{\delta}\left(X_{i}\right)^{2}\right]}$
is key for the maximal inequality in the style of \eqref{eq:Max1},
which in the current setting is given by
\[
\sqrt{\mathbb{E}\left[\ol g_{\delta}\left(X_{i}\right)^{2}\right]}=\sqrt{\mathbb{P}\left(\left|\frac{X_{i}^{'}}{\norm{X_{i}}}\theta_{0}\right|<M\delta\right)\cdot M\delta^{2}}\ =\sqrt{O\left(\delta\right)\cdot M\delta^{2}}=O\left(\delta^{\frac{3}{2}}\right),
\]
where $\mathbb{P}\left(\frac{X_{i}^{'}}{\norm{X_{i}}}\theta_0\leq M\delta\right)=O\left(\delta\right)$
follows from the observation that $\mathbb{P}\left(\frac{X_{i}^{'}}{\norm{X_{i}}}\leq M\delta\right)$
is the probability of random angle between $\frac{X_{i}^{'}}{\norm{X_{i}}}$
and $\theta_{0}$ being no more than $M\delta$ away from $\pi/2$, which scales linearly with
$\delta$ under the assumption that $p\left(x\right)$ is bounded from
above and away from zero for all $x\in{\cal X}$ in Assumption \ref{assu:Basic}(g).
In summary, $\ol g_{\delta}\left(x\right)^{2}$is at most $M\delta^{2}$ and
nonzero in a region of probability measure at most $M\delta$, and hence
$\mathbb{E}\left[\ol g_{\delta}\left(X_{i}\right)^{2}\right]$ is bounded by $M\delta^{3}$.
Importantly, $\left|x^{'}\theta\right|$ interacts multiplicatively with
the indicator function $\mathbf{\mathbbm1}\left\{ 0<x^{'}\theta_{0}<M\norm x\delta\right\} $
in \eqref{eq:Bd_gtheta}, and hence, even though indicators functions
are invariant under squaring $\mathbf{\mathbbm1}\left\{ \cdot\right\} ^{2}\equiv\mathbf{\mathbbm1}\left\{ \cdot\right\} $,
the magnitude of $\left|x^{'}\theta\right|^{2}\leq M\delta^{2}$ becomes smaller
in the order of magnitude after squaring, leading to the overall $M\delta^{3}$
on $\mathbb{E}\left[\ol g_{\delta}\left(X_{i}\right)^{2}\right]$.
To contrast this with the case of cubic-root asymptotics, say, in
\citet*{kim1990cube}, write the original maximum-score estimand $g_{MS,\theta}\left(y,x\right):=\left(y-\frac{1}{2}\right)\mathbf{\mathbbm1}\left\{ x^{'}\theta>0\right\} $,
and observe that
\begin{align*}
\left|g_{MS,\theta}\left(y,x\right)-g_{MS,\theta_{0}}\left(y,x\right)\right| & =\frac{1}{2}\cdot\left(\mathbf{\mathbbm1}\left\{ x^{'}\theta>0\geq x^{'}\theta_{0}\right\} +\mathbf{\mathbbm1}\left\{ x^{'}\theta_{0}>0\geq x^{'}\theta\right\} \right)\\
& \leq\frac{1}{2}\cdot\left\{ 0<x^{'}\theta_{0}<M\norm x\delta\right\} :=\ol g_{MS,\delta}\left(x\right)
\end{align*}
where the envelope function $\ol g_{MS,\delta}\left(x\right)$ remains
as a discrete function with
\begin{align*}
\sqrt{\mathbb{E}\left[\ol g_{MS,\delta}\left(x\right)^{2}\right]} & =\sqrt{\frac{1}{4}\mathbb{P}\left(\left|\frac{X_{i}^{'}}{\norm{X_{i}}}\theta_{0}\right|<M\delta\right)}\leq M\delta^{\frac{1}{2}},
\end{align*}
leading to a much larger bound than $\delta^{\frac{3}{2}}$ (with $\delta$
thought to be close to $0$). As discussed in \citet*{kim1990cube},
the $\delta^{\frac{1}{2}}$ bound above is the key driver for the cubic-root
asymptotics, and it arises both from the discreteness of the indicator
function $\mathbf{\mathbbm1}\left\{ x^{'}\theta>0\right\} $ as well as the discreteness
of the binary outcome $y_{i}-\frac{1}{2}$. In contrast, in our current
setting, the discrete outcome $y_{i}-\frac{1}{2}$ is replaced by
its conditional expectation, $h_{0}\left(x\right)=\mathbb{E}\left[\rest{y_{i}-\frac{1}{2}}X_{i}=x\right]$,
which is a smooth object, and furthermore the estimand $g_{+,\theta,h_{0}}\left(x\right)-g_{+,\theta_{0},h_{0}}\left(x\right)$
is constructed to be Lipschitz continuous in $x^{'}\theta$.
\subsubsection{Analysis of Term $T_{2}=\frac{1}{\sqrt{n}}\protect\mathbb{G}_{n}\left(g_{\hat{\protect\theta},\hat{h}}-g_{\protect\theta_{0},\hat{h}}-g_{\hat{\protect\theta},h_{0}}+g_{\protect\theta_{0},h_{0}}\right)$}
\noindent We now turn to the second term $T_{2}$, which involves
the first-stage nonparametric estimator $\hat{h}$ of $h$. The asymptotic
negligibility of term $T_{2}$ corresponds to the usual ``stochastic
equicontinuity'' condition, which we will seek to establish here.
To do so, we impose the following standard sup-norm convergence of
the first-stage estimator $\hat{h}$ . First, notice that given Proposition
\ref{prop:ID_Sobolev}(b), $h_{0}\in{\cal H}$ with ${\cal H}$ denoting
the space of functions mapping from ${\cal X}$ to $\left[-\frac{1}{2},\frac{1}{2}\right]$
that possess uniformly bounded derivatives up to order $d+1$. See,
for example, \citet{hansen2008uniform}, \citet{belloni2015some}
and \citet*{chen2015optimal} for results on the sup-norm convergence
of kernel and sieve nonparametric estimators.
\begin{assumption}
\label{assu:FirstStageRate}(i) $\hat{h}\in{\cal H}$ with probability
approaching 1, and (ii) $\norm{\hat{h}-h_{0}}_{\infty}=O_{p}\left(a_{n}\right)$.
\end{assumption}
\begin{lem}
\label{lem:Term2} Under Assumptions \ref{assu:Basic}-\ref{assu:FirstStageRate},
for some constant $M>0$,
\begin{equation}
P\sup_{\theta\in\Theta,h\in{\cal H}:\norm{\theta-\theta_{0}}\leq\delta,\norm{h-h_{0}}_{\infty}\leq Ka_{n}}\left|\mathbb{G}_{n}\left(g_{\theta,h}-g_{\theta_{0},h}-g_{\theta,h_{0}}+g_{\theta_{0},h_{0}}\right)\right|\leq M\delta.\label{eq:Max2}
\end{equation}
\end{lem}
Loosely speaking, Lemma \eqref{lem:Term2} implies that, whenever
$\norm{\hat{\theta}-\theta_{0}}$ converges slower than the $\sqrt{n}$ rate,
\[
T_{2}=O_{p}\left(\frac{1}{\sqrt{n}}\norm{\hat{\theta}-\theta_{0}}\right)=o_{p}\left(\norm{\hat{\theta}-\theta_{0}}^{2}\right),
\]
which will become asymptotically negligible, delivering a ``stochastic
equicontinuity'' condition that is essential for the asymptotic normality
of $\hat{\theta}$. The key model ingredient underlying this result is
the encoding of the sign restrictions via compositions of the Lipschitz-continuous
ReLU-function instead of using the discrete indicator functions as
in the formulation of the original maximum score estimator. The Lipschitz
continuity of ReLU functions, and consequently the Lipschitz continuity
of the function $g_{\theta,h}\left(x\right)=g_{+,\theta,h}\left(x\right)+g_{-,\theta,h}\left(x\right)$,
ensure that small deviations in $\theta,h$ and $x$ translate into small
deviations in $g_{\theta,h}\left(x\right)$, providing the level of smoothness
for the stochastic equicontintuity condition.
\subsubsection{Analysis of Term $T_{3}=P\left(g_{\hat{\protect\theta},h_{0}}-g_{\protect\theta_{0},h_{0}}\right)$}
\noindent Now, we turn to the third term $T_{3}=P\left(g_{\hat{\theta},h_{0}}-g_{\theta_{0},h_{0}}\right)$,
which is a familiar term that captures the quadratic curvature of
the population criterion for $\theta$ in a neighborhood of $\theta_{0}$.
Technically, the characterization of $T_{3}$ boils down to the following
second-order Taylor expansion of $Pg_{\theta,h_{0}}$ around $\theta_{0}$:
\[
P\left(g_{\theta,h_{0}}-g_{\theta_{0},h_{0}}\right)=\nabla_{\theta}Pg_{\theta_{0},h_{0}}\left(\theta-\theta_{0}\right)+\frac{1}{2}\left(\theta-\theta_{0}\right)^{'}\nabla_{\theta\theta}Pg_{\theta_{0},h_{0}}\left(\theta-\theta_{0}\right)+o\left(\norm{\theta-\theta_{0}}^{2}\right)
\]
where the gradient $\nabla_{\theta}Pg_{\theta_{0},h_{0}}$ and the Hessian $\nabla_{\theta\theta}Pg_{\theta_{0},h_{0}}$
are well-defined since $Pg_{\theta,h}$ is differentiable even though
$g_{\theta,h}$ has kinks. Moreover, since $g_{\theta,h}$ is Lipschitz-continuous
and almost surely differentiable, the gradient can be calculated easily
via $\nabla_{\theta}Pg_{\theta,h}=P\nabla_{\theta}g_{\theta,h}$ However, $\nabla_{\theta}g_{\theta,h}$
will no longer be Lipschitz-continuous and in fact involve indicator
functions, and thus the Hessian $\nabla_{\theta\theta}Pg_{\theta,h}=\nabla_{\theta}P\nabla_{\theta}g_{\theta,h}$
involves differentiation with respect to integral boundaries. As a
result, $\nabla_{\theta\theta}Pg_{\theta,h}$ becomes a ``surface integral'',
or formally, an integral over a lower-dimensional manifolds with respect
to a lower-dimensional Hausdorff measure.
Specifically, the $k$-dimensional Hausdorff measure in $\mathbb{R}^{d}$,
denoted by ${\cal H}^{k}$ for some $k\leq d$, is a ``lower-dimensional''
measure that allows us to define nontrivial integrals over lower-dimensional
subsets in $\mathbb{R}^{d}$ that has measure $0$ with respect to ${\cal L}^{d},$
the Lebesgue measure on $\mathbb{R}^{d}$. See, for example, Chapter 2 of
\citet*{evans2015measure} for the formal definition of the Hausdorff
measure. An important feature of the Hausdorff measure is the equivalence
between ${\cal H}^{k}$ and ${\cal L}^{k}$ on $\mathbb{R}^{k}$ for any $k$,
i.e., the $k$-dimensional Hausdorff measure is in some sense the
same as the Lebesgue measure on $\mathbb{R}^{k}$. On the other hand, while
a lower-dimensional space, such as a hyperplane $\left\{ x\in\mathbb{R}^{d}:x^{'}\theta_{0}=0\right\} $
in $\mathbb{R}^{d}$, is a measure-0 set with respect to ${\cal L}^{d}$ and
thus the integral $\int_{\left\{ x\in\mathbb{R}^{d}:x^{'}\theta_{0}=0\right\} }m\left(x\right)d{\cal L}^{d}\left(x\right)$
is trivially 0 for any function $m$, integrals with respect to the
$\left(d-1\right)$-dimensional Hausdorff measure of the form
\[
\int_{\left\{ x\in\mathbb{R}^{d}:x^{'}\theta_{0}=0\right\} }m\left(x\right)d{\cal H}^{d-1}\left(x\right)
\]
is nontrivial (i.e., may take values other than $0$).
\begin{lem}
\label{lem:Term3}For some positive semidefinite matrix of rank $d-1$,
we have
\[
P\left(g_{\theta,h_{0}}-g_{\theta_{0},h_{0}}\right)=-\left(\theta-\theta_{0}\right)^{'}V\left(\theta-\theta_{0}\right)+o\left(\norm{\theta-\theta_{0}}^{2}\right)
\]
with
\begin{equation}
V:=\int_{x^{'}\theta_{0}=0}\frac{f\left(\rest 0x\right)}{f\left(\rest 0x\right)+1}xx^{'}p\left(x\right)d{\cal H}^{d-1}\left(x\right)\label{eq:V_def}
\end{equation}
where ${\cal H}^{d-1}$ denotes the $\left(d-1\right)$-dimensional
Hausdorff measure in $\mathbb{R}^{d}$.
\end{lem}
Lemma \ref{lem:Term3} can be viewed as a local-identification condition,
which says that $Pg_{\theta,h_{0}}$ becomes smaller than $Pg_{\theta_{0},h_{0}}$
locally with quadratic curvature as $\theta$ moves away from the true
$\theta_{0}$. Essentially, \eqref{eq:V_def} can be viewed as a ``surface
integral'' over the $\left(d-1\right)$-dimensional hyperplane $\left\{ x\in\mathbb{R}^{d}:x^{'}\theta_{0}=0\right\} $.
Note that, even though $V$ has rank $d-1$ instead of $d$, $V$
should still be regarded to have ``full rank'' with respect to the
parameter space $\Theta=\mathbb{\mathbb{S}}^{d-1}$, which also has dimension $d-1$ instead
of $d$. This is similar to the corresponding result in \citet*{kim1990cube}.
Note that the formula of the Hessian matrix $V$ features the probability
density $f\left(\rest 0x\right)$ in the integrand, which reflects
the observation that the sign-restriction identification \ref{eq:Mono_Equiv}
is driven by the conditional median restriction and thus local in
nature. If, for example, $f\left(\rest 0x\right)=0$ for all $x\in{\cal X}$,
then the conditional median restriction is vacuous and thus identification
will fail. The dependence of the identification on $f\left(\rest 0x\right)$,
i.e., the ``conditional median density'', here is also featured
in \citet*{kim1990cube} and \citet*{horowitz1992smoothed}, as well
as more broadly in quantile regression settings. Hence, we assume
in Assumption \ref{assu:Basic} that $f\left(\rest 0x\right)$ is
bounded away from $0$.
\subsubsection{Analysis of Term $T_{4}=P\left(g_{\hat{\protect\theta},\hat{h}}-g_{\protect\theta_{0},\hat{h}}-g_{\hat{\protect\theta},h_{0}}+g_{\protect\theta_{0},h_{0}}\right)$}
The last term, $T_{4}$, reflects the influence of the first-stage
nonparametric estimation on the second-stage M-estimation criterion
function, i.e., how $P\left(g_{\hat{\theta},\hat{h}}-g_{\theta_{0},\hat{h}}\right)$
differs from $P\left(g_{\hat{\theta},h_{0}}-g_{\theta_{0},h_{0}}\right).$
This term corresponds to the derivation of the influence function
through functional differentiation in standard semiparametric two-stage
asymptotic theory.
We work with the following second-order Taylor expansion of $T_{4}$
w.r.t. $\theta$ around $\theta_{0}$:
\begin{align*}
& P\left(g_{\theta,h}-g_{\theta_{0},h}-g_{\theta,h_{0}}+g_{\theta_{0},h_{0}}\right)=P\left(g_{\theta,h}-g_{\theta,h_{0}}\right)-P\left(g_{\theta_{0},h}-g_{\theta_{0},h_{0}}\right)\\
=\ & \nabla_{\theta}P\left(g_{\theta_{0},h}-g_{\theta_{0},h_{0}}\right)\left(\theta-\theta_{0}\right)+\left(\theta-\theta_{0}\right)\nabla_{\theta\theta}P\left(g_{\theta_{0},h}-g_{\theta_{0},h_{0}}\right)\left(\theta-\theta_{0}\right)+o\left(\norm{\theta-\theta_{0}}^{2}\right).
\end{align*}
The leading term $\nabla_{\theta}P\left(g_{\theta_{0},h}-g_{\theta_{0},h_{0}}\right)$
can be linearized through pathwise functional differentiation as
\begin{equation}
\nabla_{\theta}P\left(g_{\theta_{0},h}-g_{\theta_{0},h_{0}}\right)=D_{h}\left[\nabla_{\theta}Pg_{\theta_{0},h_{0}},h-h_{0}\right]+O\left(\norm{h-h_{0}}_{\infty}\norm{\nabla_{x}\left(h-h_{0}\right)}_{\infty}\right),\label{eq:DPdg_lin}
\end{equation}
where the formula of $D_{h}\left[\nabla_{\theta}Pg_{\theta_{0},h_{0}},h-h_{0}\right]$
is derived in Lemma \ref{lem:Term4} below. With $\hat{\theta}$ and $\hat{h}$
plugged in, the term $\left(\hat{\theta}-\theta_{0}\right)\nabla_{\theta\theta}P\left(g_{\theta_{0},\hat{h}}-g_{\theta_{0},h_{0}}\right)\left(\hat{\theta}-\theta_{0}\right)$
will become asymptotically negligible provided that $\nabla_{\theta\theta}P\left(g_{\theta_{0},\hat{h}}-g_{\theta_{0},h_{0}}\right)\overset{p}{\longrightarrow}0$
holds, which can be guaranteed by the convergence of $\nabla_{x}\hat{h}$
to $\nabla_{x}h_{0}$.
\begin{assumption}
\label{assu:FirstStageDeriv} $\norm{\nabla_{x}\hat{h}-\nabla_{x}h_{0}}_{\infty}=O_{p}\left(c_{n}\right)$
with $c_{n}\searrow0$.
\end{assumption}
\begin{lem}
\label{lem:Term4}Under Assumption \ref{assu:FirstStageDeriv}\emph{,
we have}
\begin{align*}
& P\left(g_{\theta,\hat{h}}-g_{\theta_{0},\hat{h}}-g_{\theta,h_{0}}+g_{\theta_{0},h_{0}}\right)\\
= & D_{h}\left[P\nabla_{\theta}g_{\theta_{0},h_{0}},\hat{h}-h_{0}\right]^{'}\left(\theta-\theta_{0}\right)+O_{p}\left(\norm{\theta-\theta_{0}}a_{n}c_{n}\right)+o_{p}\left(\norm{\theta-\theta_{0}}^{2}\right)
\end{align*}
where
\begin{align}
D_{h}\left[\nabla_{\theta}Pg_{\theta_{0},h_{0}},h-h_{0}\right] & :=\int_{x^{'}\theta_{0}=0}\left[h\left(x\right)-h_{0}\left(x\right)\right]\frac{1}{f\left(\rest 0x\right)+1}xp\left(x\right)d{\cal H}^{d-1}\left(x\right).\label{eq:Dh_dh}
\end{align}
\end{lem}
The term $O_{p}\left(\norm{\theta-\theta_{0}}a_{n}c_{n}\right)$ will become
asymptotically negligible if $a_{n}c_{n}=o_{p}\left(\norm{\hat{\theta}-\theta_{0}}\right)$,
which can be viewed as a generalization/adaption of the usual ``$o_{p}\left(n^{-1/4}\right)$''
rate requirement on the first-stage convergence in standard semiparametric
two-stage asymptotic theory that features $n^{-1/2}$ convergence
rate for the final estimator $\hat{\theta}$. As we will show in Theorem
\ref{thm:Thm_Bin_Rate} later, the requirement $\norm{\hat{h}-h_{0}}_{\infty}=o_{p}\left(\sqrt{\norm{\hat{\theta}-\theta_{0}}}\right)$
can be satisfied under proper smoothness condition on $h_{0}$.
Note that $D_{h}\left[\nabla_{\theta}Pg_{\theta_{0},h_{0}},\hat{h}-h_{0}\right]$
can be viewed as the convergence of a plug-in estimator of lower-dimensional
integral over the nonparametric function $h_{0}$ over the hyperplane
$\left\{ x:x^{'}\theta_{0}=0\right\} $. Specifically, we can write
\[
D_{h}\left[\nabla_{\theta}Pg_{\theta_{0},h_{0}},h-h_{0}\right]=L\left(\hat{h}\right)-L\left(h_{0}\right)
\]
with
\begin{equation}
L\left(h\right):=\int_{x^{'}\theta_{0}=0}h\left(x\right)\frac{1}{f\left(\rest 0x\right)+1}xp\left(x\right)d{\cal H}^{d-1}\left(x\right).\label{eq:Lh_func}
\end{equation}
Note that $L\left(h\right)$ is a linear functional of $h$, and the
asymptotic behavior of the plugged-in estimator for linear functionals
has been widely studied in the literature on nonparametric and semiparametric
inference. While there are many results available for ``point evaluation
functionals'' and ``full-dimensional integration functionals'',
there are relatively few results for ``lower-dimensional integration
functionals'' like \eqref{eq:Lh_func}. Hence we develop results
for the asymptotic behavior of plug-in estimators of \eqref{eq:Lh_func}
in this paper.
So far we have not restricted the form of the first-stage nonparametric
estimator $\hat{h}$, and thus all our results above hold for any
form of $\hat{h}$ that satisfies Assumption \ref{assu:FirstStageRate}.
However, now we will need to be more explicit about $\hat{h}$, and
focus our attention on the Nadaraya-Watson kernel estimators and linear
series estimators, which are two leading classes of nonparametric
estimators. We provide the required conditions and results for both
classes separately below.
\subsubsection*{First Stage by Nadaraya-Watson Kernel Regression}
\noindent \textbf{Lemma 5a}
\protected@write \@auxout {}{\string \newlabel {lem:T4_Kern}{{5a}{\thepage}{5a}{lem:T4_Kern}{}} }
\hypertarget{lem:T4_Kern}{}
\emph{Under
Assumption \ref{assu:KernelSeries}(a),
\begin{equation}
D\left[P\nabla_{\theta}g_{\theta_{0},h_{0}},\hat{h}-h_{0}\right]=O_{p}\left(\frac{1}{\sqrt{nb_{n}}}+b_{n}^{s}\right)\label{eq:DPDh_rate_Kern}
\end{equation}
Setting $b_{n}\sim n^{-\frac{1}{2s+1}}$ leads to the optimal rate
of convergence $n^{-\frac{s}{2s+1}}$. With undersmoothing bandwidth
$b_{n}=o\left(n^{-\frac{1}{2s+1}}\right)$, we have
\[
\sqrt{nb_{n}}D\left(P\nabla_{\theta}g_{\theta_{0},h_{0}},\hat{h}-h_{0}\right)\overset{d}{\longrightarrow}\mathcal{N}\left({\bf 0},\Omega\right),
\]
with
\begin{align*}
\Omega & :=\int G^{2}\left(t\right)dt\cdot\int_{x^{'}\theta_{0}=0}\frac{\sigma_{0}^{2}\left(x\right)}{(f\left(\rest 0x\right)+1)^2}xx^{'}p\left(x\right)d{\cal H}^{d-1}\left(x\right),\\
G\left(t\right) & :=\int_{x^{'}\theta_{0}=0}K\left(x\right)d{\cal H}^{d-1}\left(x\right)\\
\sigma_{0}^{2}\left(x\right) & :=\text{Var}\left(\rest{Y_{i}}X_{i}=x\right)=\frac{1}{4}-h_{0}^{2}\left(x\right)
\end{align*}
}
Lemma \ref{lem:T4_Kern} shows that the asymptotics of $L\left(\hat{h}\right)$
is similar to the asymptotics of univariate nonparametric (kernel)
regressions. Specifically, the magnitude of the (square root of) variance
term in \eqref{eq:DPDh_rate_Kern} is $\left(nb_{n}\right)^{-1/2}$,
and consequently the optimal rate of convergence $n^{-\frac{s}{2s+1}}$,
do \textbf{not} depend on the dimension $d$ of the first-stage nonparametric
estimation of $h_{0}$.
This is a highly intuitive result. It is well-known from the literature
that plug-in estimators of point evaluation functionals converge at
``fully nonparametric rate'' no faster than $n^{-\frac{s}{2s+d}}$,
while plug-in estimators of (regular) ``full-dimensional integral
functionals'' converge at ``parametric rate'' $n^{-\frac{1}{2}}$,
since the ``full-dimensional integration'' effectively reduces the
dimensionality of the estimation problem by aggregating information
(and errors) over the whole $d$-dimensional support of ${\cal X}$.
Here, we are dealing a ``$\left(d-1\right)$-dimensional integral'',
which can be viewed as an intermediate case between ``point evaluation''
and ``full-dimensional integral'' functionals, and as expected our
result shows that plug-in estimators of our $\left(d-1\right)$-dimensional
integral also features an ``intermediate'' convergence rate. This
result is also consistent to the one in \citet*{newey1994kernel},
who also demonstrates accelerated convergence rates for kernel estimation
of ``partial means'', which are defined as integrals over a subvector
of $x$.\footnote{The result on partial means in \citet*{newey1994kernel} requires
that the partial means are defined with respect to a given subvector
of $x$, while our result here covers linear combinations of the whole
vector of $x$ in the form of $x^{'}\theta_{0}$.}
Lemma \ref{lem:T4_Kern} can be established by an adaption of the
proof in \citet*{newey1994kernel}. The key idea is the observation
that $G\left(t\right)$, defined as a lower-dimensional
integral of the multivariate kernel function $K$ over the $\left(d-1\right)$-dimensional
hyperplane $\left\{ x:x^{'}\theta_{0}=0\right\} $, itself qualifies as
a univariate kernel function. Furthermore, $G\left(t\right)$ is also
of smoothness order $s$. Hence, intuitively the $\left(d-1\right)$-dimensional
integral over $\left\{ x:x^{'}\theta_{0}=0\right\} $ reduces the underlying
dimensionality of the kernel nonparametric regression, thus delivering
accelerated rate of convergence for $L\left(\hat{h}\right)$ relative
to $\hat{h}$.
\subsubsection*{First Stage by Linear Series Regression}
\noindent \textbf{Lemma 5b }
\protected@write \@auxout {}{\string \newlabel {lem:T4_Series}{{5b}{\thepage}{5b}{lem:T4_Series}{}} }
\hypertarget{lem:T4_Series}{}
\emph{Under
Assumptions \eqref{assu:Basic}, and \ref{assu:KernelSeries}(b),
\[
D\left[P\nabla_{\theta}g_{\theta_{0},h_{0}},\hat{h}-h_{0}\right]=O_{p}\left(\sqrt{\frac{J_{n}}{n}}+J_{n}^{-s}\right).
\]
With $J_{n}^{-1}=o\left(n^{-\frac{1}{2s+1}}\right)$, we have
\[
\sqrt{nJ_{n}^{-1}}D\left(P\nabla_{\theta}g_{\theta_{0},h_{0}},\hat{h}-h_{0}\right)\overset{d}{\longrightarrow}\mathcal{N}\left({\bf 0},\Omega\right)
\]
for some positive semidefinite matrix with rank $d-1$ and $\theta_{0}^{'}\Omega\theta_{0}=0$.}
\subsubsection{\label{subsec:Rate}Convergence Rate and Asymptotic Normality of
$\hat{\protect\theta}$}
\noindent Now, we combine the results from Lemmas \ref{lem:Term1},
\ref{lem:Term2}, \ref{lem:Term3}, and \ref{lem:Term4} to obtain
the convergence rate of the ReLU-based estimator. In the following
we use the notation of kernel bandwidth $b_{n}$ as if the first-stage
estimator $\hat{h}$ is given by the Nadaraya-Waston kernel regression.
However, note that the arguments also apply to the setting with linear
series first stages simply with $b_{n}$ replaced by $1/J_{n}$, where
$J_{n}$ is the univariate sieve dimension (with the multivariate
sieve dimension given by $K_{n}=J_{n}^{d}$).
Plugging the implications of Lemmas \ref{lem:Term1}, \ref{lem:Term2},
\ref{lem:Term3}, and \ref{lem:Term4} into the decomposition \eqref{eq:EP_decom},
we have
\begin{align*}
0\leq\ & \mathbb{P}_{n}\left(g_{\hat{\theta},\hat{h}}-g_{\theta_{0},\hat{h}}\right)\\
\asymp\ & o_{p}\left(\norm{\hat{\theta}-\theta_{0}}^{2}\right) & T_{1}+T_{2}\\
& -\left(\hat{\theta}-\theta_{0}\right)^{'}V\left(\hat{\theta}-\theta_{0}\right)+o_{p}\left(\norm{\hat{\theta}-\theta_{0}}^{2}\right) & T_{3}\\
& +Z_{n}^{'}\left(\hat{\theta}-\theta_{0}\right)+O_{p}\left(a_{n}c_{n}\norm{\hat{\theta}-\theta_{0}}\right)+o_{p}\left(\norm{\hat{\theta}-\theta_{0}}^{2}\right) & T_{4}
\end{align*}
where
\[
Z_{n}:=D\left[P\nabla_{\theta}g_{\theta_{0},h_{0}},\hat{h}-h_{0}\right]=O_{p}\left(\frac{1}{\sqrt{nb_{n}}}+b_{n}^{s}\right).
\]
It turns out that the convergence rate of $\hat{\theta}$ is driven by
the convergence rate of $Z_{n}$ in $T_{4}$, provided that
\[
O_{p}\left(a_{n}c_{n}\norm{\hat{\theta}-\theta_{0}}\right)=o_{p}\left(\norm{\hat{\theta}-\theta_{0}}^{2}\right),
\]
i.e. $a_{n}c_{n}=o_{p}\left(\norm{\hat{\theta}-\theta_{0}}\right)$, which
can be guaranteed by a condition on $s$. Hence, $T_{1}$ and $T_{2}$
are asymptotically negligible, while $T_{3}$ and $T_{4}$ are the
asymptotically leading terms.
\section{\label{sec:MISC}General Framework: Multi-Index Single-Crossing Condition Models}
\subsection{\label{subsec:MISC_setup}RMS in the Multi-Index Single-Crossing
Framework}
We now introduce the multi-index single-crossing (MISC) condition framework as proposed in \cite{gao2020robust}, which
generalizes the single-index sign-alignment restriction \eqref{eq:Mono_Equiv} to a $J$-dimensional setting.
Formally, consider a random sample $\left(Y_{i},X_{i}\right)_{i=1}^{n}$
where $Y_{i}$ is an outcome with support ${\cal Y}\subseteq\mathbb{R}^{d_{y}}$
and
\[
X_{i}:=\left(X_{i1},...,X_{iJ}\right)\in\mathbb{R}^{d\times J}
\]
is a $d\times J$ random matrix with support ${\cal X}\subseteq\mathbb{R}^{d\times J}$.
Let $h_{0}:{\cal X}\to\mathbb{R}$ be a real-valued functional of the conditional
distribution of $Y_{i}$ given $X_{i}$ that is directly identified
and nonparametrically estimable from the data.\footnote{For example, in the binary choice model in Section \ref{sec:BinChoice},
we take $h_{0}\left(x\right)=\mathbb{E}\left[\left(Y_{i}-\frac{1}{2}\right)\mid X_{i}=x\right]$.
In other applications $h_{0}$ can be a conditional quantile, a conditional
variance, or a difference of such functionals across two states.}
We are interested in a direction parameter $\theta_{0}\in\Theta\subseteq\mathbb{\mathbb{S}}^{d-1}$
that enters the model through the $J$ parametric indexes
\[
X_{ij}^{'}\theta_{0},\qquad j=1,...,J.
\]
\begin{defn}[Multi-Index Single-Crossing Condition]
Given observable $\left(Y_{i},X_{i}\right)$ and a pair $\left(h_{0},\theta_{0}\right)$,
we say that $\left(h_{0},\theta_{0}\right)$ satisfies the (\emph{multi-index
single-crossing}) MISC condition if, for all $x=\left(x_{1},...,x_{J}\right)\in{\cal X}$,
\begin{align}
x_{j}^{'}\theta_{0}>0,\ \forall j=1,...,J & \quad\Rightarrow\quad h_{0}\left(x\right)\geq0,\nonumber \\
x_{j}^{'}\theta_{0}<0,\ \forall j=1,...,J & \quad\Rightarrow\quad h_{0}\left(x\right)\leq0.\label{eq:MISC}
\end{align}
The condition is said to be \emph{strict} if the inequalities on the
right-hand side of \eqref{eq:MISC} are strict, i.e., $h_{0}\left(x\right)>0$ whenever $x_{j}^{'}\theta_{0}>0$
for all $j$, and $h_{0}\left(x\right)<0$ whenever $x_{j}^{'}\theta_{0}<0$
for all $j$.
\end{defn}
When $J=1$, \eqref{eq:MISC} reduces exactly to the sign-alignment
restriction \eqref{eq:Mono_Equiv} used in the binary choice model
in Section \ref{sec:BinChoice}. For $J\ge2$, the MISC condition
requires the sign of $h_{0}\left(x\right)$ to align with the common
sign of the $J$ indexes whenever those indexes all agree. Importantly,
it imposes no restriction on $h_{0}\left(x\right)$ when the $J$
indexes have mixed signs.
In many applications $X_{i}$ arises as a (possibly nonlinear) transformation
of a lower-dimensional regressor $Z_{i}$, so that $X_{i}=\phi\left(Z_{i}\right)$
for some known transformation $\phi$. In that case it is convenient
to state MISC in terms of such transformed regressors.
\begin{rem}[Weak MISC with transformed regressors]
Let $W_{i}=\phi\left(X_{i}\right)$ for a known measurable map $\phi:{\cal X}\to\mathbb{R}^{d\times J}$,
and write $W_{i}=\left(W_{i1},...,W_{iJ}\right)$. We say that $\left(h_{0},\theta_{0}\right)$
satisfies the (\emph{weak}) MISC condition with respect to $W_{i}$
if, for all $x\in{\cal X}$ and $w=\phi\left(x\right)$,
\begin{align}
w_{j}^{'}\theta_{0}>0,\ \forall j=1,...,J & \quad\Rightarrow\quad h_{0}\left(x\right)\geq0,\nonumber \\
w_{j}^{'}\theta_{0}<0,\ \forall j=1,...,J & \quad\Rightarrow\quad h_{0}\left(x\right)\leq0.\label{eq:MISC-weak}
\end{align}
The strict version is defined analogously. In what follows, we suppress
the distinction when it is clear from context whether $X_{i}$ denotes
the original regressors or a transformed version.
\end{rem}
~
The RMS estimator extends naturally to the MISC framework. Given a
candidate direction $\theta\in\Theta$ and a function $h:{\cal X}\to\mathbb{R}$,
define
\begin{align}
g_{+,\theta,h}\left(x\right) & :=\left[h\left(x\right)-\min_{1\le j\le J}\left(-x_{j}^{'}\theta\right)_{+}\right]_{+},\label{eq:def_gplus_MISC}\\
g_{-,\theta,h}\left(x\right) & :=\left[-h\left(x\right)-\min_{1\le j\le J}\left(x_{j}^{'}\theta\right)_{+}\right]_{+},\label{eq:def_gminus_MISC}
\end{align}
and the population criterion
\[
Q\left(\theta\right):=Q_{+}\left(\theta\right)+Q_{-}\left(\theta\right),\qquad Q_{\pm}\left(\theta\right):=\mathbb{E}\left[g_{\pm,\theta,h_{0}}\left(X_{i}\right)\right].
\]
Then clearly, $$\theta_{0}\in\arg\max_{\theta\in\Theta}Q\left(\theta\right).$$
Intuitively, $g_{+,\theta,h}$ penalizes violations of the ``positive
sign'' restriction in \eqref{eq:MISC} when $h\left(x\right)$ is
positive but some index $x_{j}^{'}\theta$ is nonpositive, while $g_{-,\theta,h}$
penalizes violations of the ``negative sign'' restriction when $h\left(x\right)$
is negative but some index $x_{j}^{'}\theta$ is nonnegative. The inner
$\min$ and ReLU terms ensure that, for each realization $x$, only
the index that is closest to the kink at zero contributes to the loss.
Given a first-stage nonparametric estimator $\hat{h}$ of $h_{0}$,
we define the sample criterion
\[
\hat{Q}\left(\theta\right):=\frac{1}{n}\sum_{i=1}^{n}\left\{ g_{+,\theta,\hat{h}}\left(X_{i}\right)+g_{-,\theta,\hat{h}}\left(X_{i}\right)\right\}
\]
and the RMS estimator under the MISC framework as
\[
\hat{\theta}:=\arg\max_{\theta\in\Theta}\hat{Q}\left(\theta\right).
\]
The binary choice model in Section \ref{sec:BinChoice} is a strict
special case of this framework with $J=1$ and $h_{0}\left(x\right)=\mathbb{E}\left[\left(Y_{i}-\frac{1}{2}\right)\mid X_{i}=x\right]$.
In that case $g_{+,\theta,h},g_{-,\theta,h}$ reduce to the composite ReLU
functions in \eqref{eq:def_g} and the RMS estimator coincides with
the estimator studied in Section \ref{subsec:SetupResults}. When
$J\ge2$ or when $h_{0}$ is a functional other than a conditional
expectation, the traditional MS estimator cannot be applied, but the
RMS estimator remains well-defined under MISC.
\medskip{}
\noindent The MISC framework nests a large class of models, including
binary choice with awareness, selection models with multiple latent
thresholds, and panel models with multiple time indices; detailed
examples can be provided depending on the application. The key common
feature is that $h_{0}\left(x\right)$ is monotone in a common direction
$\theta_{0}$ whenever the $J$ indexes share the same sign.
To further explain the economic relevance of the MISC condition framework and the general applicability of the RMS estimator, we now provide some concrete examples\footnote{Section 4 \cite{gao2020robust} also discusses some of the examples below, as well as other examples under the MISC condition framework with endogeneity.} below along with a discussion about the related work in each specific application.
\begin{example}[Binary Choice with Awareness]
\label{exa:Bin_Aware} Consider the following modification of the
binary choice model above
\[
y_{i}=\mathbf{\mathbbm1}\left\{ X_{i1}^{'}\theta_{01}\geq u_{i}\right\} \cdot\mathbf{\mathbbm1}\left\{ X_{2i}^{'}\theta_{0}\geq v_{i}\right\}
\]
where $y_{i}$ denotes whether consumer $i$ purchases a certain,
$X_{i1}$ denotes a vector of covariates that influence the consumer's
utility from a product, and $X_{i2}$ denotes a vector of covariates
that influence the consumer's awareness of the product (such as advertising).
Here $J=2$, $X_{i}:=\left(X_{i1},X_{i2}\right)$, $W_{i1}:=X_{i1}$,
and $W_{i2}:=X_{i2}$. Let the functional $h_{0}$ be defined by $h_{0}\left(x\right):=\mathbb{E}\left[\rest{y_{i}}X_{i}=x\right]-\frac{1}{4}.$
Then, under the conditional median restrictions $\text{med}\left(\rest{u_{i}}X_{i}\right)=\text{med}\left(\rest{v_{i}}X_{i}\right)=0$
and the conditional independence restriction $\rest{u_{i}\perp v_{i}}X_{i}$,
it can be shown that
\begin{align*}
X_{i1}^{'}\theta_{01}>0,\ X_{i2}^{'}\theta_{02}>0 & \Rightarrow\quad h_{0}\left(X_{i}\right)>0,\\
X_{i1}^{'}\theta_{01}<0,\ X_{i2}^{'}\theta_{02}<0 & \Rightarrow\quad h_{0}\left(X_{i}\right)<0,
\end{align*}
again satisfying the MISC condition.
\end{example}
\begin{example}[Panel Multinomial Choice]
\label{exa:PMC} Consider the following panel multinomial choice
model studied in one of the PI's working papers \citet*{gao2020robust},
\[
y_{ijt}=\mathbf{\mathbbm1}\left\{ u\left(X_{ijt}^{'}\beta_{0},\,A_{ij},\,\epsilon_{ijt}\right)=\max_{k\in\left\{ 1,...,J\right\} }u\left(X_{ikt}^{'}\beta_{0},\,A_{ik},\,\epsilon_{ikt}\right)\right\}
\]
where $y_{ijt}$ is a binary variable indicating whether consumer
$i$ chooses product $j$ at time $t$, $X_{ijt}$ is a vector of
observable covariates, $A_{ij}$ is an unobserved fixed effect that
can be infinite dimensional, $\epsilon_{ijt}$ is an unobserved idiosyncratic
taste shock, and the utility function $u$ is assumed to be unknown
but increasing in its first argument. \citet*{gao2020robust} proposes
a novel strategy to identify and estimate the finite-dimensional parameter
$\beta_{0}$ , and the key idea is to leverage the monotonicity of $u$
to obtain a MISC condition through a sequence of intertemporal differencing
and cross-sectional averaging. Specifically, focusing on a pair of
time periods $\left(t,s\right)$ and a particular product $j_{0}$
for illustration, define $\theta_{0j}:=\beta_{0}$, $X_{i}:=\left(\left(X_{ijt}\right)_{j=1}^{J},\left(X_{ijs}\right)_{j=1}^{J}\right)$,
$h_{0}\left(X_{i}\right):=\mathbb{E}\left[\rest{y_{ij_{0}t}-y_{ij_{0}s}}X_{i}\right]$
and
\[
W_{ij}:=\begin{cases}
X_{ijt}-X_{ijs}, & j=j_{0},\\
-\left(X_{ijt}-X_{ijs}\right) & j\neq j_{0}.
\end{cases}
\]
\citet*{gao2020robust} then shows that, under quite general conditions,
the following MISC condition holds
\begin{align*}
W_{ij}^{'}\theta_{0j}>0,\ \forall j=1,...,J\quad & \Rightarrow\quad h_{0}\left(X_{i}\right)>0,\\
W_{ij}^{'}\theta_{0j}<0,\ \forall j=1,...,J\quad & \Rightarrow\quad h_{0}\left(X_{i}\right)<0.
\end{align*}
\end{example}
\begin{example}[Dyadic Network Formation]
\label{exa:NetForm} Consider the following dyadic network formation
model studied in \citet*{gao2023logical}, which is a generalization
of the one studied in \citet{graham2017econometric}:
\begin{align*}
\mathbb{E}\left[\rest{y_{ij}}X_{i},X_{j},A_{i},A_{j}\right]= & \psi\left(w\left(X_{i},X_{j}\right)^{'}\theta_{0},A_{i},A_{j}\right)
\end{align*}
Here $y_{ij}$ is a binary outcome indicating whether individuals
$i$ and $j$ are linked in an undirected network, $X_{i}$ and $X_{j}$
are the individuals' observable covariates, $w\left(X_{i},X_{j}\right)$
is a known pairwise transformation of individual covariates (with
the leading example being $w_{h}\left(X_{i},X_{j}\right):=\left|X_{i,h}-X_{j,h}\right|$
for each coordinate $h=1,...,d_{x}$), $A_{i}$ and $A_{j}$ are unobserved
individual degree heterogeneity terms, and $\psi:\mathbb{R}^{3}\to\mathbb{R}$ is
an unknown function assumed to be multivariate increasing in all its
three arguments. \citet*{gao2023logical} proposes a method, called
``logical differencing'', to cancel out the unobserved heterogeneity
terms $A_{i}$ despite the lack of additive separability in the model,
a technical complication that arises naturally under nontransferable
utility settings. Specifically, fixing a particular pair of individuals
$\ol i$ and $\ol j$ and two generic realizations $\ol x,\ul x$
of $X_{i}$, it can be shown that, with
\[
\ol w:=w\left(x_{\ol j},\ol x\right)-w\left(x_{\ol i},\ol x\right),\quad\ul w:=w\left(x_{\ol i},\ul x\right)-w\left(x_{j},\ul x\right),
\]
and
\begin{align*}
h_{0}\left(\ol x,\ul x\right):= & \left(\mathbb{E}\left[\rest{y_{\ol ik}-y_{\ol jk}}X_{k}=\ol x\right]\right)_{+}\mathbb{E}\left[\rest{y_{\ol ik}-y_{\ol jk}}X_{k}=\ul x\right],\\
& -\left(\mathbb{E}\left[\rest{y_{\ol jk}-y_{\ol ik}}X_{k}=\ol x\right]\right)_{+}\mathbb{E}\left[\rest{y_{\ol jk}-y_{\ol ik}}X_{k}=\ul x\right]
\end{align*}
the weak MISC condition is satisfied (under quite mild additional
conditions):
\begin{align*}
\ol w^{'}\theta_{0}>0,\ul w^{'}\theta_{0}>0\quad & \Rightarrow\quad h_{0}\left(\ol x,\ul x\right)\geq0,\\
\ol w^{'}\theta_{0}<0,\ul w^{'}\theta_{0}<0\quad & \Rightarrow\quad h_{0}\left(\ol x,\ul x\right)\leq0.
\end{align*}
\end{example}
\begin{example}[Conditional Quantile Model for Continuous Outcomes]
Consider the following model
\[
y_{i}=\phi\left(X_{i}^{'}\theta+\epsilon_{i}\right),\quad\text{med}\left(\rest{\epsilon_{i}}X_{i}\right)=0,
\]
where $\phi$ is some unknown strictly increasing function. If ${\bf 0}\in Supp\left(X_{i}\right)$,
we can take $h_{0}$ to be the difference in conditional median functions
\[
h_{0}\left(x\right):=\text{med}\left(\rest{y_{i}}X_{i}=x\right)-\text{med}\left(\rest{y_{i}}X_{i}=0\right),
\]
so that \eqref{eq:MISC} holds since
\begin{align*}
\text{med}\left(\rest{y_{i}}X_{i}=x\right) & =\phi\left(\text{med}\left(\rest{X_{i}^{'}\theta+\epsilon_{i}}X_{i}=x\right)\right)\\
& =\phi\left(x_{i}^{'}\theta+\text{med}\left(\rest{\epsilon_{i}}X_{i}=x\right)\right)=\phi\left(x_{i}^{'}\theta\right).
\end{align*}
Alternatively, we could also state the single-crossing condition in
terms of pairwise differences by
\[
h_{0}\left(\ol x,\ul x\right):=\text{med}\left(\rest{y_{i}}X_{i}=\ol x\right)-\text{med}\left(\rest{y_{i}}X_{i}=\ul x\right)
\]
so that
\[
h_{0}\left(\ol x,\ul x\right)\lessgtr0\quad\Leftrightarrow\quad\left(\ol x-\ul x\right)^{'}\theta\lessgtr0,
\]
which is a special case of \eqref{eq:MISC} with $J=2$ and $g\left(\ol x,\ul x\right)=\ol x-\ul x.$
\end{example}
\begin{example}[Stochastic Volatility for Continuous Outcomes]
Consider the following simple ``stochastic volatility'' model of
some centered (mean-zero) variable $y_{t}$:
\begin{align*}
y_{t} & =\sigma\left(X_{t}^{'}\theta+\epsilon_{t}\right)\cdot u_{t}
\end{align*}
where $\sigma$ is some unknown strictly increasing function and $u_{t}$
is mean-zero exogenous error with $\mathbb{E}\left[\rest{u_{t}^{2}}X_{t}\right]=1$.
Suppose that $\epsilon_{t}\perp\left(X_{t},u_{t}\right)$. Then we can
set
\begin{align*}
h_{0}\left(\ol x,\ul x\right) & :=\mathbb{E}\left[\rest{y_{t}^{2}}X_{t}=\ol x\right]-\mathbb{E}\left[\rest{y_{t}^{2}}X_{t}=\ul x\right]\\
& =\mathbb{E}\left[\sigma^{2}\left(\ol x^{'}\theta+\epsilon_{t}\right)-\sigma^{2}\left(\ul x^{'}\theta+\epsilon_{t}\right)\right]
\end{align*}
so that
\[
h_{0}\left(\ol x,\ul x\right)\lessgtr0\quad\Leftrightarrow\quad\left(\ol x-\ul x\right)^{'}\theta\lessgtr0.
\]
\end{example}
\noindent It should be pointed out that the above are just a few illustrations
of many plausible econometric models nested under the MISC condition
framework. Given that the exact specifications of $y,X,\phi,h_{0}$
are left largely unrestricted, they can be user-configured in very
flexibly manners depending on the economic contexts: for example,
$X$ can be decomposed into an ``endogenous/structural'' part and
an ``exogenous/IV'' part, while $W=\phi\left(X\right)$ and $h_{0}\left(X\right)$
may involve a subvector or the whole of $X$ with potentially nonlinear
transformations.
One main advantage of the MISC framework lies in its ability to identify
and estimate index parameters in models with rich forms of unobserved
heterogeneity and additively nonseparable interactions between modeling
ingredients.
\subsection{RMS Asymptotic Theory under MISC}
\label{subsec:MISC-asymptotics}
We now derive the convergence rate and asymptotic distribution of the RMS
estimator $\hat\theta$ in the multi-index single-crossing (MISC) framework of
Section~3.1. As in the single-index case, the key ingredients are: (i) a
linearization of the effect of first-stage estimation errors through a
lower-dimensional submanifold integral functional, and (ii) a local quadratic
expansion of the population criterion $Q(\theta)$ around $\theta_0$. In the
MISC case, both objects have a particularly transparent form.
Recall that
\[
Q(\theta;h)
:= P g_{\theta,h}, \qquad
Q(\theta) := Q(\theta;h_0) = P g_{\theta,h_0},
\]
and define the (vector-valued) directional derivative functional
\[
L(h)
:= D_h\big(P \nabla_\theta g_{\theta_0,h_0}\big)\big[h-h_0\big]
\in\mathbb{R}^d.
\]
Throughout this subsection, we view $L(h)$ as a map on a suitable function
space $H$ containing $h_0$ and the first-stage estimator $\hat h$.
The next lemma collects the two structural properties that drive the
asymptotics: a submanifold-integral representation of the linear functional
$L(h)$ and a quadratic expansion of $Q(\theta)$ around $\theta_0$.
\begin{lem}[Asymptotics via Submanifold Integrals]
\label{lem:MISC-curv-L}
Under the strict MISC condition~\textup{(20)} hold,
\begin{enumerate}
\item[(a)]
For any $c\in\mathbb{R}^d$, define the scalar functional
\[
\Gamma_c(h)
:= c' P\nabla_\theta g_{\theta_0,h}.
\]
Then $\Gamma_c$ is Fr\'echet differentiable at $h_0$ and its derivative
satisfies
\begin{equation}
D_h\Gamma_c(h_0)[v]
= \sum_{j=1}^J \int_{\{x\in\mathcal{X}: x_j'\theta_0=0\}}
v(x)\, w_{c,j}(x)\, d\mathcal{H}^{d-1}(x),
\qquad \forall v\in H,
\label{eq:MISC-submanifold-deriv}
\end{equation}
for some uniformly bounded weight functions $w_{c,j}:\mathcal{X}\to\mathbb{R}$.
In particular, each component of $L(h)$ can be written as a finite sum of
integrals of $(h-h_0)$ over the hyperplanes $\{x:x_j'\theta_0=0\}$.
\item[(b)]
There exists a symmetric positive semidefinite $d\times d$ matrix $V$ of rank
$d-1$ such that, for all $\theta$ in a neighborhood of $\theta_0$ with
$\|\theta\|=1$,
\begin{equation}
Q(\theta) - Q(\theta_0)
= -(\theta-\theta_0)'V(\theta-\theta_0)
+ o\big(\|\theta-\theta_0\|^2\big),
\label{eq:MISC-quadratic}
\end{equation}
and $V\theta_0=0$. Moreover, $V$ admits the representation
\begin{equation}
V
= \sum_{j=1}^J \int_{\{x\in\mathcal{X}: x_j'\theta_0=0\}}
m_j(x,\theta_0)\, x_j x_j' p(x)\, d\mathcal{H}^{d-1}(x),
\label{eq:MISC-V-rep}
\end{equation}
for some nonnegative Lipschitz functions $m_j(\cdot,\theta_0)$,
$j=1,\dots,J$, and $(d-1)$-dimensional Hausdorff measure
$\mathcal{H}^{d-1}$.
\end{enumerate}
\end{lem}
Lemma~\ref{lem:MISC-curv-L} shows that the second-stage curvature is governed by a $(d-1)$-dimensional matrix $V$ and that the first-stage impact enters only through submanifold integrals of $h-h_0$ over those $(d-1)$-dimensional hyperplanes. This is precisely the setting analyzed in \cite{chen2025semiparametric}, with submanifold dimension $m=d-1$ (codimension $d-m=1$).
\begin{assumption}
\label{ass:MISC-CG}
Suppose that:
\begin{enumerate}
\item[(i)] The true function $h_0$ belongs to a H\"older (or Sobolev) ball of
smoothness order $s>1$ on a compact support $\mathcal{X}\subset\mathbb{R}^d$.
\item[(ii)] The first-stage estimator $\hat h$ is either a kernel or linear series
(sieve) estimator constructed as in Section~2, with smoothing parameter
(bandwidth or sieve dimension) chosen so that the conditions of Assumptions~6–8
in Chen and Gao (2025) hold for the regressors $X_i$ and the basis. In
particular, if $K_n$ denotes the sieve dimension, then
\[
K_n\log K_n / n \to 0
\quad\text{and}\quad
K_n^{-s/d} = o\Big(\sqrt{K_n^{(d-1)/d}/n}\Big).
\]
\item[(iii)] For each $c\in\mathbb{S}^{d-1}$, the scalar functional
$\Gamma_c(h)=c'P\nabla_\theta g_{\theta_0,h}$ satisfies the linearization
and regularity conditions in Assumptions~9–11 of \cite{chen2025semiparametric} with
submanifold dimension $m=d-1$ and level-set function
$g_j(x)=x_j'\theta_0$, $j=1,\dots,J$.
\end{enumerate}
\end{assumption}
Assumption~\ref{ass:MISC-CG}(c) is essentially a restatement, in our notation,
of the high-level conditions required to apply Theorems~2 and~3 of \cite{chen2025semiparametric} to the functionals $c'L(h)$. Under these conditions, those
theorems yield both the convergence rate and the asymptotic normality of
$L(\hat h)$ as an estimator of $L(h_0)$.
We can now state the main result of this subsection.
\begin{thm}[RMS Asymptotics under MISC]
\label{thm:Asymp-MISC}
Suppose the MISC condition \eqref{eq:MISC}, and Assumption~\ref{ass:MISC-CG} hold.
\begin{enumerate}
\item[(a)]
For the linear functional $L(h)$ defined above, under undersmoothing,
\begin{equation}
c_n\big(L(\hat h) - L(h_0)\big)
\;\xrightarrow{d}\;
\mathcal{N}(0,\Omega),
\label{eq:Lhat-CLT}
\end{equation}
with $c_n$ can be taken to be slower than but arbitrarily close to $n^{-s/(2s+1)}$.
\item[(b)]
Let $\hat\theta$ denote the RMS estimator under MISC,
Then
\[
c_n(\hat\theta - \theta_0)
= - V^{-} c_n L(\hat h) + o_p(1),
\]
where $V$ is the Hessian in \eqref{eq:MISC-quadratic} and $V^{-}$
its Moore–Penrose inverse restricted to the tangent space orthogonal to $\theta_0$. Consequently, with undersmoothing,
\begin{equation}
c_n(\hat\theta - \theta_0)
\;\xrightarrow{d}\;
\mathcal{N}\big(0,\,V^{-}\Omega V^{-}\big).
\label{eq:thetahat-CLT}
\end{equation}
\end{enumerate}
\end{thm}
\begin{rem}[Effective one-dimensional rate in the $J$-index case]
\label{rem:MISC-1D-rate}
By Lemma~\ref{lem:MISC-curv-L}(ii), the submanifold functional $L(h)$ depends
on $h$ only through its restriction to the $(d-1)$-dimensional hyperplanes
$\{x:x_j'\theta_0=0\}$, $j=1,\dots,J$. The analysis in \cite{chen2025semiparametric}
shows that, for kernel or sieve estimators of $h_0$ on a $d$-dimensional
support, the minimax-optimal rate for such submanifold integrals is
$n^{-s/(2s+1)}$, independent of $J$. Thus $c_n=n^{s/(2s+1)}$ in
Theorem~\ref{thm:Asymp-MISC}, and the RMS estimator under MISC achieves the
same “one-dimensional” nonparametric rate as in the single-index binary choice
model. Increasing $J$ affects only the constants and the asymptotic variance
matrix $V^{-}\Omega V^{-}$, not the convergence rate.
\end{rem}
\section{DNN-Based Maximum Score Estimator}
\label{sec:NN}
In this section we show how the RMS estimator can be further adapted to be implemented within a neural network architecture. The key observation is that the RMS criterion is itself a composition of ReLU units with a simple, interpretable structure. This allows us to view the RMS estimator as a special multi-layer network with a dedicated ``RMS layer'' that extracts the sign information of the index parameter $\theta$, and to estimate $\theta$ using standard machine learning software.
\subsection{RMS as a Special Neural Network Layer}
\label{subsec:NN-layer}
We first describe the single-index binary choice model. Let $x\in\mathbb{R}^d$ denote
the covariate and recall that in Section~\ref{sec:BinChoice} we defined, for a
generic function $h$ and direction $\theta\in\Theta\subset\mathbb{\mathbb{S}}^{d-1}$,
\[
g_{+,\theta,h}(x)
:= \bigl[h(x) - [-x'\theta]_+\bigr]_+,
\qquad
g_{-,\theta,h}(x)
:= \bigl[-h(x) - [x'\theta]_+\bigr]_+,
\]
and the RMS population criterion $Q(\theta)=E[g_{+,\theta,h_0}(X_i) +
g_{-,\theta,h_0}(X_i)]$. These maps are compositions of three elementary
operations:
\begin{enumerate}
\item a \emph{directional projection} $s(x;\theta)=x'\theta$;
\item a \emph{sign-extracting pair} of ReLU units $[s(x;\theta)]_+$ and
$[-s(x;\theta)]_+$; and
\item a final \emph{RMS transform} that compares $h(x)$ to the
ReLU-transformed index via an outer ReLU.
\end{enumerate}
This structure can be encoded as a small neural network module
$R_\theta(h)(x)$ that takes as input the scalar $h(x)$ and the vector $x$,
computes $x'\theta$, passes it through ReLUs, and outputs $g_{+,\theta,h}(x)$
and $g_{-,\theta,h}(x)$ (or their difference). In particular, for any fixed
$\theta$, $h\mapsto R_\theta(h)$ is a Lipschitz, piecewise linear operator.
A convenient way to embed RMS into a network is to treat $h$ as the output of a
generic multi-layer perceptron $f_\beta:\mathbb{R}^d\to\mathbb{R}$ with parameters
$\beta\in\mathbb{R}^p$, and then apply the RMS layer to $(x,f_\beta(x))$. In notation,
set
\[
g_{+}(x;\theta,\beta)
:= \bigl[f_\beta(x) - [-x'\theta]_+\bigr]_+,
\qquad
g_{-}(x;\theta,\beta)
:= \bigl[-f_\beta(x) - [x'\theta]_+\bigr]_+,
\]
and define
\[
h_{\theta,\beta}(x) := g_{+}(x;\theta,\beta) - g_{-}(x;\theta,\beta).
\]
The map $x\mapsto h_{\theta,\beta}(x)$ is then a neural network with one
special ``RMS layer'' on top of a generic (deep) regression network $f_\beta$. When
$\theta=\theta_0$ and $f_\beta$ approximates $h_0$, the outputs
$(g_{+},g_{-})$ implement the same sign-alignment structure as in the
population RMS criterion, and the resulting $h_{\theta,\beta}$ inherits the
economic interpretation of $h_0$.
\subsection{DNN-Based MISC Estimation}
\label{subsec:NN-MISC}
In the $J$-index MISC setting of Section~\ref{sec:MISC}, the relevant
population criterion is defined based on the following: for each $x=(x_1,\dots,x_J)$,
\[
g_{+,\theta,h_0}(x)
= \Bigl[h_0(x) - \bigl(\min_{1\le j\le J}(-x_j'\theta)_+\bigr)\Bigr]_+,
\qquad
g_{-,\theta,h_0}(x)
= \Bigl[-h_0(x) - \bigl(\min_{1\le j\le J}(x_j'\theta)_+\bigr)\Bigr]_+.
\]
which can be encoded in a neural network with the following special architecture:
\begin{enumerate}
\item A MLP neural network to approximate $h_0$.
\item A \emph{multi-index generation layer} that computes the $J$ scalar indexes
$s_j(x;\theta)=x_j'\theta$ and their ReLU transforms
$[s_j(x;\theta)]_+$, $[-s_j(x;\theta)]_+$.
\item A \emph{MISC aggregation layer} that takes the elementwise minimum
\[
u(x;\theta) := \min_j [-s_j(x;\theta)]_+,
\qquad
v(x;\theta) := \min_j [s_j(x;\theta)]_+,
\]
and passes them, together with $h(x)$, through outer ReLUs as above.
\end{enumerate}
The resulting multi-layer neural network encodes exactly the MISC conditions as in Section~\ref{sec:MISC}. The MISC parameter $\theta$ appears only in the linear projections $x_j'\theta$ inside this special layer, while the possibly high-dimensional parameters $\beta$ govern flexible, nonparametric features through $h(x)=f_\beta(x)$.
From an applied perspective, one of the main appeals of the DNN-based MISC formulation
is that it provides a principled way to extract an economically meaningful index parameter $\theta$ from a high-dimensional black-box DNN.
\subsection{Implementation using Machine Learning Packages}
\label{subsec:NN-implementation}
The network architectures described above are straightforward to implement in
standard machine learning frameworks such as \texttt{PyTorch} or
\texttt{TensorFlow}. The main ingredients are:
\begin{itemize}
\item a base MLP $f_\beta$ with ReLU activation (possibly deep),
\item a directional parameter $\theta$ constrained to lie on the unit sphere,
implemented via explicit normalization or a reparameterization, and
\item a custom ``RMS layer'' that takes $(x,f_\beta(x),\theta)$ as input and
outputs $g_{+}(x;\theta,\beta)$ and $g_{-}(x;\theta,\beta)$.
\end{itemize}
Since all components are compositions of affine maps and ReLU activations,
the network is differentiable almost everywhere and compatible with automatic
differentiation. Training can therefore be carried out using standard
gradient-based optimizers (e.g.\ ADAM) with GPU acceleration.
There are two natural training strategies:
\begin{enumerate}
\item \emph{Two-step RMS:} First estimate $h_0$ by training $f_\beta$ to
minimize a standard loss (e.g.\ squared error between $Y_i$ and
$f_\beta(X_i)$). Then plug in $\hat h(x)=f_{\hat\beta}(x)$ and optimize
$\hat Q(\theta)$ over $\theta$ only, using the RMS layer as in
Sections~\ref{sec:BinChoice} and~\ref{sec:MISC}.
\item \emph{Joint DNN:} Parameterize the outcome as
$Y_i\approx h_{\theta,\beta}(X_i)$ via the RMS layer and estimate both
$\theta$ and $\beta$ jointly by minimizing a loss such as
$\frac{1}{n}\sum_i (Y_i-h_{\theta,\beta}(X_i))^2$ subject to $\|\theta\|=1$.
This corresponds to embedding the MISC structure directly into a deep network
and training it with standard backpropagation.
\end{enumerate}
The two-step approach falls directly under our existing asymptotic theory, once
$\hat h$ is shown to satisfy the first-stage conditions. The joint-DNN approach is more demanding theoretically but conceptually attractive, as it treats $\theta$ as a low-dimensional ``interpretable head'' on top of a deep, flexible feature extractor.
Formally establishing the asymptotic properties of $\hat{\theta}$ in the joint DNN estimation approach is an interesting direction for future research. One natural route would be to show that, under suitable conditions on the loss, architecture and regularization, the joint estimator of $\theta$ is asymptotically equivalent to the two-step/profile RMS estimator studied here, given that the MISC parameter $\theta$ only shows up in the ``outmost'' hidden layer of the DNN. An alternative route would be to use sample-splitting or cross-fitting to obtain valid inference for $\theta$ directly from the joint optimization problem.
\section{Simulation}
\label{sec:Sim}
Our goal in this section is to investigate the finite-sample performance of the
RMS estimator $\hat\theta$ for $\theta_0$ in both the single-index binary
choice model and the two-index MISC setting. We first describe the common
simulation design and implementation, and then discuss an alternative neural
network implementation that embeds the MISC structure directly into the network
architecture.
\subsection{Simulation Design and Implementation}
\label{subsec:SimDesign}
Each Monte Carlo experiment follows the same basic four-step procedure:
\begin{enumerate}
\item Generate a random sample from a given data-generating process (DGP).
\item Obtain an estimate $\hat{\theta}$ either using a two-step plug-in procedure or the joint DNN procedure.
\item Evaluate the performance of $\hat\theta$ across $B$ Monte Carlo
replications.
\end{enumerate}
\subsubsection{DGP Specification}
\paragraph{Single-index DGP.}
In the baseline design we consider the binary choice model
\[
y_i = 1\{X_i'\theta_0 > \varepsilon_i\},
\]
with
\[
\theta_0
= \Bigl(\tfrac{\sqrt{3}}{3}, -\tfrac{\sqrt{3}}{3}, \tfrac{\sqrt{3}}{3}\Bigr)',
\qquad \| \theta_0 \| = 1.
\]
The regressors are drawn independently as
$X_{i1}, X_{i2}, X_{i3} \sim \mathrm{Unif}[-2,2]$, and the error terms
$\varepsilon_i$ are i.i.d.\ logistic. Denoting by $F$ the logistic CDF, the
true first-stage function is
\[
h_0(x)
:= E\bigl[y_i - \tfrac{1}{2} \mid X_i=x\bigr]
= F(x'\theta_0) - \tfrac{1}{2}
= \frac{1}{1 + \exp(-x'\theta_0)} - \tfrac{1}{2},
\]
which is known in closed form but treated as unknown in the estimation
procedure.
\paragraph{Two-index (MISC) DGP.}
To illustrate the multi-index setting, we also consider a two-index model
($J=2$) that satisfies the MISC condition. For each $i$, we generate
\[
y_i
= 1\{X_{i1}'\theta_0 > \varepsilon_{i1}\}\,
1\{X_{i2}'\theta_0 > \varepsilon_{i2}\},
\]
where $\varepsilon_{i1},\varepsilon_{i2}$ are i.i.d.\ logistic and each
component of $X_{i1}$ and $X_{i2}$ is i.i.d.\ $\mathrm{Unif}[-2,2]$. Writing
$X_i = (X_{i1},X_{i2})$ and using the same $\theta_0$ as above, we have
\[
P(y_i=1 \mid X_i=(x_1,x_2))
= F(x_1'\theta_0)F(x_2'\theta_0),
\]
so that
\[
h_0(x_1,x_2)
:= E\bigl[y_i - \tfrac{1}{4} \mid X_{i1}=x_1,X_{i2}=x_2\bigr]
= F(x_1'\theta_0)F(x_2'\theta_0) - \tfrac{1}{4}.
\]
This DGP satisfies the strict MISC condition: $h_0(x_1,x_2)>0$ whenever both
$x_1'\theta_0$ and $x_2'\theta_0$ are positive, and $h_0(x_1,x_2)<0$ whenever
both are negative.
\subsubsection{Two-Stage Implementation}
\paragraph{First-Stage Nonparametric Regression}
Given simulated data, we estimate $h_0$ nonparametrically by regressing
$y_i - \tfrac{1}{2}$ on $X_i$ in the single-index design, and
$y_i - \tfrac{1}{4}$ on $(X_{i1},X_{i2})$ in the two-index design. We consider
two main classes of estimators (implemented using standard R packages):
\begin{itemize}
\item \emph{Kernel regression}, with a polynomial kernel and bandwidth
selected over a small grid (e.g.\ using a rule of thumb or simple
cross-validation). In the reported simulations we use a polynomial kernel with
tuning parameters $\alpha = 0.1$ and $\gamma = 0.0001$.
\item \emph{Series (sieve) regression}, based on tensor-product spline bases,
with the number of basis functions playing the role of the smoothing parameter.
\item \emph{Neural Network regression}: a standard multi-layer perceptron (MLP) with
ReLU activation, where the main tuning parameters are the number of hidden
units and layers. In the experiments reported below, a typical configuration
uses a hidden size of 10, 2 hidden layers, a learning rate of 0.01, and 100
epochs of training with the ADAM optimizer.
\end{itemize}
\paragraph{Second-stage optimization of the RMS criterion.}
Given $\hat h$, we form the sample analogue of the RMS criterion,
\[
\hat Q(\theta)
:= \frac{1}{n}\sum_{i=1}^n
\Bigl\{ g_{+,\theta,\hat h}(X_i) + g_{-,\theta,\hat h}(X_i) \Bigr\},
\]
and maximize $\hat Q(\theta)$ over
$\theta$ on the unit sphere $\{\theta:\theta'\theta=1\}$. We use a
gradient-based algorithm (ADAM) together with a simple projection step to
enforce the unit-norm constraint. In practice, this amounts to running ADAM
updates on the unconstrained parameter vector and renormalizing $\theta$ to
unit length after each update. The learning rate is set to 0.01 and we run 500
epochs for each replication. The use of ReLU functions makes the objective
continuous and Lipschitz in $\theta$, so gradients are well defined almost
everywhere and standard optimization routines are stable in these simulations.
\subsubsection{Joint Implementation via Neural Networks}
We also consider the DNN-based joint estimation of $h_0$ and $\theta_0$ as described in Section \ref{sec:NN}. Specifically, we use a three-stage training strategy:
\begin{itemize}
\item \textbf{Stage 1}: Freeze $\theta$ parameters (initialized to zero vectors), and train only the MLP component parameters to learn basic function approximation.
\item \textbf{Stage 2}: Freeze the MLP component, reinitialize and train only the directional parameter $\theta$.
\item \textbf{Stage 3}: Jointly train all parameters for fine-tuning.
\end{itemize}
\subsubsection{Performance Measures}
For each design, we consider two sample sizes $N \in \{1000,5000\}$, and re report summary measures of the distribution of $\hat\theta$ across $B=1000$ Monte Carlo replications. The basic componentwise diagnostics are the Monte Carlo mean squared error (MSE), bias, and standard deviation (SD) of each
coordinate of $\hat\theta$. To capture overall performance in a
rotation-invariant way, we also report: the $\ell^1$ error of each coordinate, the $\ell^2$ norm of the bias vector, and mean/median ``angular similarity'', defined as one minus the cosine of the angle between $\hat\theta$ and $\theta_0$.
\subsection{Results}
\subsubsection{Single-Index DGP}
For the single-index design, Tables 1–3 report the performance of the RMS
estimator with three different first-stage implementations: kernel regression
(Table 1), a separate neural network nonparametric estimator (Table 2), and an
“all-in-one’’ neural network that jointly estimates the first stage and
$\theta$ (Table 3). In all cases, increasing the sample size from $N=1000$ to
$N=5000$ substantially reduces MSE, standard deviations, and angular errors:
$1-\text{mean angular similarity}$ falls from roughly $6\times 10^{-3}$ to
$4\times 10^{-3}$ for the kernel, and from about $1.0\times 10^{-2}$ to
$3\text{–}4\times 10^{-3}$ for the neural network implementations. The kernel
first stage is slightly more accurate than the neural network alternatives at
$N=1000$, but by $N=5000$ all three approaches deliver very similar accuracy,
with small biases and tight angular concentration around $\theta_0$.
\begin{table}
\centering
\caption{Two-Stage RMS with Kernel First Stage}
\begin{tabular}{lcc}
\toprule
\textbf{Metric} & \textbf{N=1000} & \textbf{N=5000} \\
\midrule
MSE of $\theta_1$ & 0.00375 & 0.00249 \\
MSE of $\theta_2$ & 0.00405 & 0.00228 \\
MSE of $\theta_3$ & 0.00395 & 0.00257 \\
Bias of $\theta_1$ & -0.00279 & -0.00157 \\
Bias of $\theta_2$ & 0.00380 & 0.00154 \\
Bias of $\theta_3$ & -0.00358 & -0.00325 \\
SD of $\theta_1$ & 0.06114 & 0.04986 \\
SD of $\theta_2$ & 0.06356 & 0.04771 \\
SD of $\theta_3$ & 0.06274 & 0.05058 \\
L1 Error of $\theta_1$ & 0.04643 & 0.03657 \\
L1 Error of $\theta_2$ & 0.04886 & 0.03562 \\
L1 Error of $\theta_3$ & 0.04620 & 0.03615 \\
\midrule
L2 Norm of Bias & 0.005923 & 0.003920 \\
1-- Mean Angular Similarity & 0.005874 & 0.003668 \\
1-- Median Angular Similarity & 0.003250 & 0.001833 \\
\bottomrule
\end{tabular}
\end{table}
\begin{table}
\centering
\caption{Two-Stage RMS with Neural-Net First Stage}
\begin{tabular}{lcc}
\toprule
\textbf{Metric} & \textbf{N=1000} & \textbf{N=5000} \\
\midrule
MSE of $\theta_1$ & 0.00703 & 0.00260 \\
MSE of $\theta_2$ & 0.00728 & 0.00241 \\
MSE of $\theta_3$ & 0.00678 & 0.00242 \\
Bias of $\theta_1$ & -0.00410 & -0.00328 \\
Bias of $\theta_2$ & 0.00640 & 0.00312 \\
Bias of $\theta_3$ & -0.00776 & -0.00004 \\
SD of $\theta_1$ & 0.08374 & 0.05085 \\
SD of $\theta_2$ & 0.08508 & 0.04903 \\
SD of $\theta_3$ & 0.08195 & 0.04924 \\
L1 Error of $\theta_1$ & 0.06505 & 0.03846 \\
L1 Error of $\theta_2$ & 0.06705 & 0.03713 \\
L1 Error of $\theta_3$ & 0.06458 & 0.03759 \\
\midrule
L2 Norm of Bias & 0.010862 & 0.004529 \\
1-- Mean Angular Similarity & 0.010542 & 0.003717 \\
1-- Median Angular Similarity & 0.006961 & 0.002064 \\
\bottomrule
\end{tabular}
\end{table}
\begin{table}
\centering
\caption{Joint DNN-Based Estimation}
\begin{tabular}{lcc}
\toprule
\textbf{Metric} & \textbf{N=1000} & \textbf{N=5000} \\
\midrule
MSE of $\theta_1$ & 0.01047 & 0.00275 \\
MSE of $\theta_2$ & 0.00997 & 0.00260 \\
MSE of $\theta_3$ & 0.00990 & 0.00285 \\
Bias of $\theta_1$ & -0.00856 & -0.00174 \\
Bias of $\theta_2$ & 0.01187 & 0.00184 \\
Bias of $\theta_3$ & -0.00584 & -0.00352 \\
SD of $\theta_1$ & 0.10198 & 0.05238 \\
SD of $\theta_2$ & 0.09916 & 0.05097 \\
SD of $\theta_3$ & 0.09933 & 0.05327 \\
L1 Error of $\theta_1$ & 0.07945 & 0.04148 \\
L1 Error of $\theta_2$ & 0.07802 & 0.04006 \\
L1 Error of $\theta_3$ & 0.07886 & 0.04231 \\
\midrule
L2 Norm of Bias & 0.015764 & 0.004335 \\
1-- Mean Angular Similarity & 0.015174 & 0.004099 \\
1-- Median Angular Similarity & 0.009594 & 0.002806 \\
\bottomrule
\end{tabular}
\end{table}
\subsubsection{Two-Index Design}
For the two-index MISC design, Tables 4–6 show the same three
implementations. The problem is clearly harder: MSEs and angular errors are larger than in the single-index case, though they still improve remarkably with sample size. Here the choice of first-stage method matters more. The kernel version (Table 4) achieves reasonable performance, but the two-step neural network first stage (Table 5) delivers substantially smaller MSE and angular error, especially at $N=5000$ (where MSEs drop from about $10^{-2}$ to roughly
$3\times 10^{-3}$, and $1-\text{mean angular similarity}$ from about $1.5\times 10^{-2}$ to around $4.6\times 10^{-3}$). The all-in-one neural network (Table 6) performs similarly to the kernel in this two-index setting and does\ not match the accuracy of the two-step neural network. Overall, the tables confirm that (i) the RMS estimator behaves in line with the theory as $N$ grows, (ii) the two-step architecture is robust and competitive in the single-index case, and (iii) in more complex multi-index designs, flexible neural network first stages can yield clear gains over standard kernel
smoothing.
\begin{table}
\centering
\caption{Two-Stage RMS with Kernel First Stage: $J=2$}
\begin{tabular}{lcc}
\toprule
\textbf{Metric} & \textbf{N=1000} & \textbf{N=5000} \\
\midrule
MSE of $\theta_1$ & 0.02901 & 0.01047 \\
MSE of $\theta_2$ & 0.02810 & 0.01003 \\
MSE of $\theta_3$ & 0.02805 & 0.00995 \\
Bias of $\theta_1$ & -0.02186 & -0.01029 \\
Bias of $\theta_2$ & 0.02827 & 0.00720 \\
Bias of $\theta_3$ & -0.02361 & -0.00887 \\
SD of $\theta_1$ & 0.16890 & 0.10178 \\
SD of $\theta_2$ & 0.16522 & 0.09989 \\
SD of $\theta_3$ & 0.16580 & 0.09937 \\
L1 Error of $\theta_1$ & 0.12860 & 0.07589 \\
L1 Error of $\theta_2$ & 0.12873 & 0.07545 \\
L1 Error of $\theta_3$ & 0.12967 & 0.07591 \\
\midrule
L2 Norm of Bias & 0.042832 & 0.015379 \\
1-- Mean Angular Similarity & 0.042575 & 0.015223 \\
1-- Median Angular Similarity & 0.029377 & 0.008477 \\
\bottomrule
\end{tabular}
\end{table}
\begin{table}
\centering
\caption{Two-Stage RMS with Neural-Net First Stage: $J=2$}
\begin{tabular}{lcc}
\toprule
\textbf{Metric} & \textbf{N=1000} & \textbf{N=5000} \\
\midrule
MSE of $\theta_1$ & 0.01890 & 0.00303 \\
MSE of $\theta_2$ & 0.02165 & 0.00314 \\
MSE of $\theta_3$ & 0.01637 & 0.00302 \\
Bias of $\theta_1$ & -0.01421 & -0.00304 \\
Bias of $\theta_2$ & 0.02280 & 0.00113 \\
Bias of $\theta_3$ & -0.01229 & -0.00380 \\
SD of $\theta_1$ & 0.13674 & 0.05498 \\
SD of $\theta_2$ & 0.14536 & 0.05606 \\
SD of $\theta_3$ & 0.12736 & 0.05485 \\
L1 Error of $\theta_1$ & 0.09536 & 0.04228 \\
L1 Error of $\theta_2$ & 0.10210 & 0.04332 \\
L1 Error of $\theta_3$ & 0.09336 & 0.04255 \\
\midrule
L2 Norm of Bias & 0.029541 & 0.004994 \\
1-- Mean Angular Similarity & 0.028461 & 0.004600 \\
1-- Median Angular Similarity & 0.013980 & 0.002799 \\
\bottomrule
\end{tabular}
\end{table}
\begin{table}
\centering
\caption{Joint DNN-Based Estimation: $J=2$}
\begin{tabular}{lcc}
\toprule
\textbf{Metric} & \textbf{N=1000} & \textbf{N=5000} \\
\midrule
MSE of $\theta_1$ & 0.02751 & 0.01160 \\
MSE of $\theta_2$ & 0.02881 & 0.01188 \\
MSE of $\theta_3$ & 0.02689 & 0.01183 \\
Bias of $\theta_1$ & -0.02517 & -0.00710 \\
Bias of $\theta_2$ & 0.01779 & 0.01367 \\
Bias of $\theta_3$ & -0.02910 & -0.00981 \\
SD of $\theta_1$ & 0.16394 & 0.10745 \\
SD of $\theta_2$ & 0.16881 & 0.10816 \\
SD of $\theta_3$ & 0.16137 & 0.10834 \\
L1 Error of $\theta_1$ & 0.12984 & 0.08602 \\
L1 Error of $\theta_2$ & 0.13342 & 0.08658 \\
L1 Error of $\theta_3$ & 0.13054 & 0.08681 \\
\midrule
L2 Norm of Bias & 0.042388 & 0.018265 \\
1-- Mean Angular Similarity & 0.041604 & 0.017658 \\
1-- Median Angular Similarity & 0.026760 & 0.012385 \\
\bottomrule
\end{tabular}
\end{table}
\section{\label{sec:Con}Conclusion}
We have proposed a rectified-linear-unit-based maximum score (RMS) estimator for models characterized by sign-alignment restrictions. By replacing the
discontinuous indicator in Manski’s maximum score with composite ReLU functions,
the population criterion becomes piecewise smooth with quadratic curvature,
while preserving the underlying identification logic. This structure delivers
an intermediate, ``one-dimensional'' rate $n^{-s/(2s+1)}$ and asymptotic
normality, but also yields a sample objective that is much more amenable to
modern gradient-based optimization methods. In practice, RMS can be optimized
using off-the-shelf routines from machine learning, avoiding the fragile,
combinatorial searches often required for discontinuous maximum score criteria.
We also embed the binary choice model in a general multi-index single-crossing
(MISC) framework, where several indexes enter through a common direction
parameter. Even in this multi-index setting, the leading term in the
asymptotic expansion depends on the nonparametric component only through
its restriction to a finite union of $(d-1)$-dimensional hyperplanes, so the
effective nonparametric dimension remains one and the convergence rate is
unchanged. Taken together, these results show that ReLU-based formulations can
retain the robustness and partial identification features of maximum score,
while offering significant computational advantages and extending naturally to
richer multi-index environments.
\bibliographystyle{ecca}
\bibliography{ReMS}