ReLU-Based and DNN-Based Generalized Maximum Score Estimators

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{#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 {6pt} {6pt}

Introduction

In a sequence of papers, \citet*{manski1975maximum,manski1985semiparametric} proposed and analyzed the properties of the maximum-score estimator in the context of semiparametric discrete choice models. To be specific, consider the following canonical binary choice model

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,

which is a 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,

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

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 (ref) above. However, instead of using multiplication with indicator functions on $X_{i}^{'}\theta_{0}$ as in (ref), our new formulation employs compositions of ReLU functions. Specifically, define

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

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 (ref),

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

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

The corresponding $J$-index RMS population criterion is $Q_J(\theta) := Q_J^+(\theta)+Q_J^-(\theta)$ with

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 (ref) and (ref).

The main focus of this paper is to show how this new ReLU-based population criterion $Q$, as defined by (ref)-(ref), 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), 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 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 kim1990cube, and it coincides with the rate of the 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 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 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 gao2020robust and 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) 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 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 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 fan2021interpretability and zhang2021survey review a wide range of interpretability tools, which mostly focus on explaining 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) 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) 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) 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) presents simulation evidence on the finite-sample performance of the RMS estimator in both single-index and multi-index designs. Section (ref) concludes. Technical proofs and additional auxiliary results are collected in the appendix.

Special Case: Binary Choice Model

In this section, we focus on the binary choice model (ref) 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.

Setup and Main Results

Given the binary choice model (ref) and the ReLU-based population criterion function $Q$ in (ref), we define the ReLU-based maximum score (RMS) estimator as

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 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 (ref) using the the ReLU function $\left[\cdot\right]_{+}$, which is Lipschitz continuous and everywhere differentiable 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 (ref) 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.

Assumption (ref)(a) and (b) consists of a standard random-sampling assumption for the binary choice model (ref) 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 (ref), our proposed method can also be adapted to other settings. See XXX for a more detailed discussion.

Assumption (ref)(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 (ref) and the conditional expectation function $h_{0}\left(x\right)$ beyond manski1985semiparametric and 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)(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)(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)(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)(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)(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)(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)(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) below:

Given mild convergence conditions on the first-stage estimator $\hat{h}$, it is straightforward to establish the consistency of $\hat{\theta}$ in Theorem (ref).

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.

We now present our main results about the RMS asymptotics.

The asymptotic distribution can then be derived based on the linearized argmax theorem (Theorem 3.2.16) in \citet*{van1996weak}.

Outline of the RMS Asymptotic Theory

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 (ref) and (ref) 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

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.

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) below presents a maximal inequality about $T_{1}$ with respect to $\theta$ in a small neighborhood of $\theta_{0}$:

Loosely speaking, the result above in Lemma (ref) 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 (ref). 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

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) 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 (ref), 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)(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 (ref), 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

where the envelope function $\ol g_{MS,\delta}\left(x\right)$ remains as a discrete function with

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$.

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)$

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)(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, hansen2008uniform, belloni2015some and \citet*{chen2015optimal} for results on the sup-norm convergence of kernel and sieve nonparametric estimators.

Loosely speaking, Lemma (ref) 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.

Analysis of Term $T_{3}=P\left(g_{\hat{\protect\theta},h_{0}}-g_{\protect\theta_{0},h_{0}}\right)$

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$).

Lemma (ref) 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, (ref) 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) 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) that $f\left(\rest 0x\right)$ is bounded away from $0$.

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}$:

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

where the formula of $D_{h}\left[\nabla_{\theta}Pg_{\theta_{0},h_{0}},h-h_{0}\right]$ is derived in Lemma (ref) 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}$.

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) 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

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 (ref). Hence we develop results for the asymptotic behavior of plug-in estimators of (ref) 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). 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.

First Stage by Nadaraya-Watson Kernel Regression

Lemma 5a \protected@write \@auxout {\string \newlabel {lem:T4_Kern}{{5a}{\thepage}{5a}{lem:T4_Kern}} } \hypertarget{lem:T4_Kern} Under Assumption (ref)(a),

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

}

Lemma (ref) 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 (ref) is $\left(nb_{n}\right)^{-1/2}$, and consequently the optimal rate of convergence $n^{-\frac{s}{2s+1}}$, do 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) 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}$.

First Stage by Linear Series Regression

Lemma 5b \protected@write \@auxout {\string \newlabel {lem:T4_Series}{{5b}{\thepage}{5b}{lem:T4_Series}} } \hypertarget{lem:T4_Series} Under Assumptions (ref), and (ref)(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$.

Convergence Rate and Asymptotic Normality of $\hat{\protect\theta}$

Now, we combine the results from Lemmas (ref), (ref), (ref), and (ref) 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), (ref), (ref), and (ref) into the decomposition (ref), we have

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.

General Framework: Multi-Index Single-Crossing Condition Models

RMS in the Multi-Index Single-Crossing Framework

We now introduce the multi-index single-crossing (MISC) condition framework as proposed in gao2020robust, which generalizes the single-index sign-alignment restriction (ref) 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), 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. \]

When $J=1$, (ref) reduces exactly to the sign-alignment restriction (ref) used in the binary choice model in Section (ref). 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.

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

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 (ref) 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) 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 (ref) and the RMS estimator coincides with the estimator studied in Section (ref). 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.

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 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.

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.

RMS Asymptotic Theory under MISC

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$.

Lemma (ref) 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 chen2025semiparametric, with submanifold dimension $m=d-1$ (codimension $d-m=1$).

Assumption (ref)(c) is essentially a restatement, in our notation, of the high-level conditions required to apply Theorems 2 and 3 of 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.

DNN-Based Maximum Score Estimator

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.

RMS as a Special Neural Network 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) 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:

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$.

DNN-Based MISC Estimation

In the $J$-index MISC setting of Section (ref), 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:

The resulting multi-layer neural network encodes exactly the MISC conditions as in Section (ref). 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.

Implementation using Machine Learning Packages

The network architectures described above are straightforward to implement in standard machine learning frameworks such as PyTorch or TensorFlow. The main ingredients are:

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:

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.

Simulation

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.

Simulation Design and Implementation

Each Monte Carlo experiment follows the same basic four-step procedure:

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.

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):

\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.

Joint Implementation via Neural Networks

We also consider the DNN-based joint estimation of $h_0$ and $\theta_0$ as described in Section (ref). Specifically, we use a three-stage training strategy:

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$.

Results

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$.

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.

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.