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.
37,390 characters
Gaussian Approximation for Maximum Score and Non-Smooth M-Estimators with Multiway Dependence
\begin{abstract}
The maximum score estimator of \citet{manski1975maximum} provides an elegant approach to estimate slope coefficient in binary choice models without requiring parametric assumptions on the error distribution. However, under i.i.d. sampling, it admits a non-Gaussian limiting distribution and exhibits cube-root asymptotics, which complicates statistical inference.
We show that, under multiway dependence, the maximum score estimator attains asymptotic normality at a parametric rate. We obtain this surprising result through the development of a general M-estimation theory that accommodates non-smooth objective functions under multiway dependence. We further propose and establish the validity of a bootstrap procedure for inference.
\end{abstract}
\maketitle
\section{Introduction}
The maximum score estimator of \citet{manski1975maximum} offers a robust and conceptually elegant method for estimating parametric coefficients in binary response models under minimal assumptions. In contrast to conventional parametric approaches, such as maximum likelihood estimation, it does not require specification of the error distribution or independence between covariates and unobservables, relying instead solely on the ordinal information encoded in the sign of the latent index. This distribution-free feature renders the estimator particularly attractive from a theoretical standpoint, as it is robust to misspecification of the disturbance term and remains valid in settings where standard parametric assumptions are untenable.
Despite these appealing properties, the estimator exhibits fundamentally non-standard asymptotic behaviour under i.i.d.\ sampling. In a seminal contribution, \citet{kim1990cube} show that the maximum score estimator converges at the cube-root rate $n^{1/3}$, rather than the usual parametric rate, and possesses a non-Gaussian Chernoff-type limiting distribution characterised by the argmax of a stochastic process. This irregular asymptotic structure arises from the non-smooth, discontinuous nature of the objective function, which invalidates the classical central limit theory and the standard bootstrap is known to fail \citep{abrevaya2005bootstrap}.
These features pose substantial challenges for statistical inference. In particular, the absence of asymptotic normality complicates the construction of confidence intervals and hypothesis tests, while the slower convergence rate implies reduced precision in finite samples. As a consequence, practical implementation often requires nonstandard inference procedures tailored to cube-root asymptotics-- see, for example, \cite{delgado2001subsampling,patra2018consistent,cattaneo2020bootstrap}--which further limits the estimator’s accessibility and widespread empirical use.
In this paper, we show that, under multiway dependence (also known as multiway clustering)—a complex dependence structure commonly encountered in empirical research \citep{miglioretti2007marginal,petersen2008estimating,cameron2011robust,thompson2011simple,cameron2015practitioner,mackinnon2023cluster}—the original maximum score estimator of \citet{manski1975maximum} attains asymptotic normality with a parametric convergence rate. This finding stands in sharp contrast to the classical i.i.d.\ setting and suggests that, in this context, dependence can be beneficial rather than detrimental for inference.
Our maximum score results are built upon a general asymptotic theory for non-smooth M-estimators under multiway dependence, which may be of independent interest.
A substantial body of work has developed asymptotic theory under multiway dependence; see, for example, \citet{menzel2021bootstrap}, \citet{2021daveziesEmpiricalProcessResults,davezies2025analytic}, \citet{chiang2023inference}, and \citet{graham2024sparse}. Despite these advances, a general asymptotic framework for $M$-estimators—whether based on smooth or non-smooth objective functions—remains largely unexplored in this setting. The present paper contributes to this literature by developing a unified theory for $M$-estimation under multiway dependence, thereby filling this gap.
Despite asymptotic normality (equivalently, asymptotic Gaussianity), analytical inference is challenging because the asymptotic variance depends on derivatives of smooth population objects. These derivatives are difficult to estimate because their sample counterparts are non-smooth and do not provide direct analogues. We circumvent these challenges by proposing a procedure based on a multiplier bootstrap. Building upon our Gaussian approximation theory for M-estimators, this approach provides a robust and practical framework for implementation.
The intuition behind our Gaussian approximation result is related to the asymptotic theory for
$U$-process-based estimators, such as simplicial depth \citep{liu1990notion} and Oja's spatial medians \citep{oja1983descriptive}, established in \cite{arcones1994estimators}, as well as the notion of data-generating-process (DGP)–induced smoothing studied in panel quantile regression with common shocks by \citet{chiang2026panel}. However, the present setting is substantially more complex. Unlike the convex and subgradient-friendly quantile check loss, the maximum score objective is inherently discontinuous, which renders standard arguments inapplicable. In particular, even when a smooth approximation exists, it is non-trivial to show that the maximiser of the original objective is well approximated by that of the smoothed problem. Addressing this issue requires a careful quadratic approximation argument, drawing on
$U$-process M-estimation ideas and the stochastic differentiability framework of \cite{pollard1985new}.
A further complication arises in establishing stochastic equicontinuity for higher-order projections. In the classical
$U$-process settings of \cite{arcones1994estimators}, such results rely on weak convergence for degenerate processes \citep{arcones1993limit}, which is not available under multiway dependence. To overcome this, we develop an alternative approach combining localised empirical process methods with an iterated argument that first delivers a preliminary rate and then sharpens it, using maximal inequalities to control the resulting entropy integrals \citep{chen2026cross}.
Building on these ingredients, we establish a general asymptotic theory for non-smooth M-estimators under multiway dependence. By carefully localising the effect of DGP-induced smoothing, we control the discrepancy between the original objective and its smooth approximation, which in turn yields asymptotic linearity, a Gaussian approximation, and a valid bootstrap procedure for the maximum score estimator at the parametric rate.
A number of approaches have been proposed to address the non-Gaussian limiting distribution of the maximum score estimator and the attendant challenges for statistical inference. In a seminal contribution, \citet{horowitz1992smoothed} proposes the smoothed maximum score estimator, whereby the discontinuous indicator function in the objective is replaced with a smooth kernel approximation. This modification renders the objective function differentiable and delivers asymptotic normality under i.i.d.\ sampling, albeit at a nonparametric rate determined by the bandwidth choice. While this restores inferential tractability, it introduces an additional tuning parameter and entails smoothing bias, both of which require careful calibration in practice. More recently, \citet{chen2025relu} extend this line of work by employing the ReLU functions to encode the sign alignment restrictions, obtaining asymptotic normality at a nonparametric rate faster than $n^{-1/3}$, though similar issues concerning tuning and approximation bias persist.
An alternative strand of the literature exploits the maximum rank correlation approach. For example, \citet{lee1999root} considers a two-period panel model and show that, following an appropriate pairwise differencing transformation, the estimation problem can be recast as a maximum rank correlation estimator of the type analysed by \citet{sherman1993limiting}. This reformulation allows them to establish asymptotic normality under suitable regularity conditions. Closely related is \citet{honore2000panel}, which studies a conditional maximum score estimator for panel data and similarly relies on smoothing.
Nevertheless, all the aforementioned modifications alter the original estimator. For Manski's maximum score estimator, Gaussian asymptotics remain unavailable, and the inherent non-standard behaviour continues to pose a fundamental challenge for inference.
Our contribution demonstrates that the unmodified maximum score estimator of \cite{manski1975maximum} can, in fact, be asymptotically Gaussian under conditions of complex multiway dependence.
\subsection*{Notation}
Let $\mathbb{N}$ denote the set of positive integers and $\mathbb{R}$ denote the real line. For $a,b\in \mathbb{R}$, let $a\vee b=\max\{a,b\}$ and $a\wedge b = \min\{a,b\}$. Throughout the paper, the asymptotics are taken with respect to $n\to\infty$ for $n$ that will be defined in \Cref{sec:maximum_score}. For vectors $\bm{a},\bm{b}\in\mathbb{R}^K$, write $\bm{a}\le\bm{b}$ for coordinate-wise inequality.
\(\|\cdot\|\) is the Euclidean norm unless otherwise specified. Henceforth, $\leadsto$ denotes weak convergence and $\overset{\mathrm{p}}{\to}$ denotes convergence in probability.
\section{Maximum Score Estimator}\label{sec:maximum_score}
In this section, we study the asymptotic properties of the maximum score estimator under multiway dependence.
Let \(K \in \mathbb{N}\), \({\bm{N}} = \left( N_{1}, \dots, N_{K} \right) \in
\mathbb{N}^{K}\) and \(n = \min \left\{ N_{1}, \dots, N_{K} \right\}\).
Each \(N_{k}\) is the sample size in dimension \(k\).
Let \(I_{{\bm{N}}} = \prod_{k = 1}^{K} \left\{ 1, \dots, N_{k}
\right\}\) be the set of $K$-dimensional indices.
The econometrician observes the data \(\mathcal{D}_{{\bm{N}}} = \left\{ W_{{\bm{i}}} \right\}_{{\bm{i}} \in I_{{\bm{N}}}} =
\left\{ \left( Y_{{\bm{i}}}, X_{{\bm{i}}}^{\prime} \right)^{\prime} \right\}_{{\bm{i}} \in
I_{{\bm{N}}}}\), where \(Y_{{\bm{i}}}\) is a binary outcome with values in \(\{- 1, 1\}\)
and \(X_{{\bm{i}}}\) is a \(d\)-dimensional vector of covariates.
Let \(\mathbb{\mathbb E}_{{\bm{N}}}\) denote sample averaging over \(\mathcal{D}_{{\bm{N}}}\), so that
\(\mathbb{\mathbb E}_{{\bm{N}}} [f (W)] = \left| I_{{\bm{N}}} \right|^{-1} \sum_{{\bm{i}} \in
{\mathcal{I}}_{{\bm{N}}}} f \left( W_{{\bm{i}}} \right)\).
Denote the surface of the unit sphere by
\(\mathbbm{S}^{d - 1} = \left\{ v :
\|v\| = 1 \right\}\) and
\(\mathcal{B} \subseteq \mathbbm{S}^{d - 1}\).
The maximum score estimator solves
\begin{equation}
\widehat{Q}_{{\bm{N}}} \left( \widehat{\beta}_{{\bm{N}}} \right) = \sup_{b \in \mathcal{B}}
\widehat{Q}_{{\bm{N}}} (b) - o_P (1)
= \sup_{b \in \mathcal{B}} \mathbb{\mathbb E}_{{\bm{N}}} \left[ Y \cdot \perp\!\!\!\perp \left\{ X^{\prime}
b \geq 0 \right\} \right] - o_P (1).
\label{eqn--max-score-estimator-def}
\end{equation}
For each $k=1,\dots,K$, define
\(
\mathcal{E}_k=\{{\bm{e}}\in\{0,1\}^K:\|{\bm{e}}\|_0=k\},
\)
so that $\{0,1\}^K=\bigcup_{k=0}^K\mathcal{E}_k$ and let ${\bm{e}}_k$ denote the $k$-th unit vector. In other words, $\mathcal{E}_1=\{{\bm{e}}_k:k=1,...,K\}$ and $\mathcal{E}_2=\{{\bm{e}}=(e_1,...,e_K)\in \{0,1\}^K: \sum_{k=1}^K e_k=2\}$.
We impose the following standard assumption for multiway dependence, which is also used in \cite{menzel2021bootstrap} and \cite{2021daveziesEmpiricalProcessResults,davezies2025analytic}, among others.
\begin{assumption}
\label{asm--exchangeable}
\(\left\{ W_{{\bm{i}}} \right\}_{{\bm{i}} \in\mathbb N^K}\) are separately exchangeable and
dissociated.
Equivalently, each \(W_{{\bm{i}}}\) admits the Aldous--Hoover--Kallenberg
representation:
there is a Borel measurable \(\tau (\cdot)\) and
mutually independent latent \(\mathrm{Uniform} [0, 1]\) variables
\(\left\{ U_{\bm{i} \odot \bm{e}} : \bm{i} \in
\mathbb{N}^{K}, \bm{e} \in \{0, 1\}^{K} \setminus \{\bm{0}\} \right\}\) such that
\begin{align}
W_{{\bm{i}}} \overset{d}{=} \tau \left( \left\{ U_{\bm{i} \odot \bm{e}} \right\}_{\bm{e} \in \{0,
1\}^{K} \setminus \{\bm{0}\}} \right),
\label{eqn--AHK-representation}
\end{align}
where $\odot$ denotes the Hadamard product.
\end{assumption}
Under \Cref{asm--exchangeable}, the \(W_{{\bm{i}}}\)'s are typically not mutually independent,
but identically distributed.
Let \(W = (Y, X^{\prime})^{\prime}\) denote a random variable distributed
according to the common marginal distribution of \(W_{{\bm{i}}}\), and denote
\begin{equation}
m_{0} (X) = \mathbb{\mathbb E} [Y | X].
\label{eqn--m0-def}
\end{equation}
We impose the following modelling and regularity conditions on the joint distributions of $(Y,X')$.
\begin{assumption}
\label{asm--max-score-consistency}
\hfill
\begin{enumerate}[(i)]
\item \label{asm--latent-linear-median-regression}
\(Y = 2 \perp\!\!\!\perp \left\{ X^{\prime} \beta_{0} - \epsilon \geq 0 \right\} -
1\) for some \(\beta_{0} \in \mathbbm{S}^{d - 1}\), and the conditional median of the error satisfies \(\mathrm{Med}
[\epsilon | X] = 0\) almost surely.
\item \label{asm--index-sgn-disagree-pos-pr}
For any \(b \in \mathbbm{S}^{d - 1} \setminus \beta_{0}\), \(\text{Pr} \left\{ \left|
m_{0} (X) \right| > 0, \,\mathrm{sgn} \left( X^{\prime} \beta_{0} \right) \neq \mathrm{sgn}
\left( X^{\prime} b \right) \right\} > 0\).
\item \label{asm--continuity-under-constraint}
The estimation set \(\mathcal{B}\) in \eqref{eqn--max-score-estimator-def}
is a closed subset of \(\mathbbm{S}^{d - 1}\).
Furthermore for every \(b \in \mathcal{B}\), \(\text{Pr} \left\{ X^{\prime} b =
0 \right\} = 0\).
\end{enumerate}
\end{assumption}
Lower-level sufficient conditions for \Cref{asm--max-score-consistency}
\eqref{asm--index-sgn-disagree-pos-pr} and
\eqref{asm--continuity-under-constraint} are well known in the literature.
For example \citet{1985manskiSemiparametricAnalysisDiscrete} assumes an index
restriction coupled with a special regressor.
In particular, Assumption 2 of \citet{1985manskiSemiparametricAnalysisDiscrete}
requires there to exist an index \(l \in \{1, \dots, k\}\) such that \(\beta_{0,
l} \neq 0\) and the corresponding regressor, \(X [l]\), is assumed to have full
support on \(\mathbb{R}\) (almost surely) conditional on the other regressors.
In proving consistency, Theorem 1 of
\citet{1985manskiSemiparametricAnalysisDiscrete} further restricts the effective
parameter space to be \(\mathcal{B} = \left\{ b \in \mathbbm{S}^{d - 1} : \left|
b_{l} \right| \geq \eta \right\}\), for some \(0 < \eta < \beta_{0, l}\).
Together with some additional regularity conditions, these all imply
\Cref{asm--max-score-consistency}
\eqref{asm--index-sgn-disagree-pos-pr} and
\eqref{asm--continuity-under-constraint}.
Whilst the assumptions of \citet{1985manskiSemiparametricAnalysisDiscrete} allow
for some discrete (or constant) regressors, \citet{kim1990cube}
directly assume absolute continuity of \(X\) with respect to Lebesgue measure
with a differentiable density, and further assume that \(X / \|X\|\) has a
density with respect to surface measure on the sphere.
These imply \Cref{asm--max-score-consistency}
\eqref{asm--continuity-under-constraint} with \(\mathcal{B} = \mathbbm{S}^{d - 1}\).
The following results provides consistency of the maximum score estimator under multiway dependence.
\begin{lemma}
\label{lem--max-score-consistency}
Suppose \Cref{asm--exchangeable,asm--max-score-consistency} hold, and further assume
the restricted estimation set \(\mathcal{B}\) in
\eqref{eqn--max-score-estimator-def} satisfies \(\beta_{0} \in \mathcal{B}\) for
\(\beta_{0}\) in
\Cref{asm--max-score-consistency} \eqref{asm--latent-linear-median-regression}.
Then \(\widehat{\beta}_{n}\) in \eqref{eqn--max-score-estimator-def} satisfies
\(\widehat{\beta}_{n} \overset{\mathrm{p}}{\to} \beta_{0}\).
\end{lemma}
\begin{proof}[Proof of \Cref{lem--max-score-consistency}]
See \Cref{sec--prf--lem--max-score-consistency}.
\end{proof}
We now consider asymptotic Gaussianity.
Consistency of \(\widehat{\beta}_{{\bm{N}}}\) allows us to focus on values of \(\beta \in
\mathcal{B}\) close to \(\beta_{0}\).
The following local parameter space will be particularly useful for our purpose:
\begin{equation*}
\mathcal{B}_{0} = \{\beta \in \mathbbm{S}^{d - 1} : \|\beta - \beta_{0}\| \leq
\sqrt{2}\} = \{\beta \in \mathbbm{S}^{d - 1} : \beta_{0}' \beta \geq 0\}.
\end{equation*}
\(\mathcal{B}_{0}\) admits a smooth parametrisation by
\(\Theta = \left\{ \theta \in \mathbb{R}^{d - 1} : \|\theta\| \leq 1 \right\}\).
Let \(B_{0}\) be an orthonormal \(d \times (d - 1)\) basis matrix
for the subspace \(\{\beta : \beta^{\prime} \beta_{0} = 0\}\),
i.e. \(B_{0}^{\prime} B_{0} = \mathbb{I}_{d - 1}\), the $(d-1)$-dimensional identity
matrix, and \(\beta_{0}^{\prime} B_{0}\) is the \((d - 1)\)-dimensional zero
vector.
Define
\begin{equation}
\beta (\theta) = B_{0} \theta + \sqrt{1 - \|\theta\|^{2}} \cdot \beta_{0}.
\label{eqn--beta-local}
\end{equation}
The parametrisation in \eqref{eqn--beta-local} is a diffeomorphism from
\(\Theta\) to \(\mathcal{B}_{0}\) and satisfies \(\beta_{0} = \beta (0)\), \(\|\beta
(\theta)\| = 1\) for every \(\theta\), and \(\beta_{0}^{\prime} \beta (\theta) =
\sqrt{1 - \|\theta\|^{2}}\).
Thus as an added bonus, \(\beta_{0}\) is ``interior'' since it is identified
with \(\theta = 0 \in \mathrm{int} (\Theta)\).
Finally, let
\begin{equation*}
\widehat{\theta}_{{\bm{N}}} = B_{0}^{\prime} \widehat{\beta}_{{\bm{N}}} = B_{0}^{\prime} \left(
\widehat{\beta}_{{\bm{N}}} - \beta_{0} \right).
\end{equation*}
Then \(\widehat{\theta}_{{\bm{N}}} \overset{\mathrm{p}}{\to} 0\) by \(\widehat{\beta}_{{\bm{N}}} \overset{\mathrm{p}}{\to} \beta_{0}\).
It can be shown that
\begin{equation*}
\widehat{\beta}_{{\bm{N}}} = B_{0} \widehat{\theta}_{{\bm{N}}} + \mathrm{sgn} \left( \widehat{\beta}_{{\bm{N}}}'
\beta_{0} \right) \sqrt{1 - \|\widehat{\theta}_{{\bm{N}}}\|^{2}} \cdot \beta_{0}.
\end{equation*}
Hence, \(\widehat{\beta}_{{\bm{N}}} = \beta \left( \widehat{\theta}_{{\bm{N}}} \right)\) on the
event \(\left\{ \widehat{\beta}_{{\bm{N}}} \in \mathcal{B}_{0} \right\}\).
We further impose the following regularity conditions on the joint distributions
of $Y$, $X$, and the latent shocks $(U_{{\bm{e}}'})_{{\bm{e}}'\le {\bm{e}}}$ arising from
\eqref{eqn--AHK-representation}, for a collection of indices ${\bm{e}}\in \mathcal{E}_1 \cup
\mathcal{E}_2$ that will be specified below.
Notably, these assumptions do not impose topological structure on the latent
shocks.
\begin{assumption}
\label{asm--X-ac-smooth-m-diff}
The following approximate maximiser condition holds at rate $n^{-1}$:
\[\sup_{\theta \in \Theta}\widehat{Q}_{{\bm{N}}} ( \beta (\theta)
) - \widehat{Q}_{{\bm{N}}} ( \widehat{\beta}_{{\bm{N}}}) = o_P (n^{-1}).\] Furthermore,
denote $\mathcal{U}_{\bm{e}}=\{U_{{\bm{e}}'}\}_{{\bm{e}}'\le {\bm{e}}}$ and $\mathbf{u}_{\bm{e}}=\{u_{{\bm{e}}'}\}_{{\bm{e}}'\le {\bm{e}}}$. The following conditions hold.
\begin{enumerate}[(i)]
\item \label{asm--X-ac-smooth-m-diff-p-smooth}
For each $\bm{e} \in \mathcal{E}_2$, the conditional distribution \(X \mid \mathcal{U}_{\bm{e}}\) is absolutely continuous
against Lebesgue measure, and the associated conditional density, \(p_{\bm{e}} (x | \bm{u}_{\bm{e}})\), is
continuously differentiable with respect to \(x\) on the interior of its
support. Furthermore, $\mathbb{\mathbb E}[\|X\|^3]<\infty$.
\item \label{asm--X-ac-smooth-m-diff-m-smooth}
For each $\bm{e} \in \mathcal{E}_2$, let \(m_{0, \bm{e}} (x, \bm{u}_{\bm{e}}) = \mathbb{\mathbb E} \left[ Y \middle| X = x, \mathcal{U}_{\bm{e}}=\bm{u}_{\bm{e}}
\right]\).
Then \(m_{0, \bm{e}} (x, \bm{u}_{\bm{e}})\) is continuously differentiable with respect to
\(x\) for each \(\bm{u}_{\bm{e}}\) with derivative \(\dot{m}_{0, \bm{e}} (x, \bm{u}_{\bm{e}}) =
\frac{\partial}{\partial x} m_{0, \bm{e}} (x, \bm{u}_{\bm{e}})\).
\item \label{asm--X-ac-smooth-m-diff-m-dominance}
For each $\bm{e} \in \mathcal{E}_2$, there exists \(m_{\ast} (x, \bm{u}_{\bm{e}})\) that is square-integrable
against \(\left( X, \mathcal{U}_{\bm{e}} \right)\) and
\begin{equation*}
\max \left\{
p_{\bm{e}} (x | \bm{u}_{\bm{e}}), m_{0, \bm{e}} (x, \bm{u}_{\bm{e}}), |\dot{m}_{0, \bm{e}} (x, \bm{u}_{\bm{e}})| \right\} \leq
m_{\ast} (x, \bm{u}_{\bm{e}}).
\end{equation*}
\item \label{asm--X-ac-smooth-m-diff-m-increasing}
For each \({\bm{e}}_k \in \mathcal{E}_1\), \(m_{0, {\bm{e}}_k} (x, u) = 0\) if \(x^{\prime} \beta_{0} = 0\) and
the map \(t \mapsto m_{0, {\bm{e}}_k} (x + t \beta_{0}, u)\) is strictly
increasing in \(t\) almost surely with respect to the joint distribution of
\(\left( X, U_{\bm{e}_{k}} \right)\).
\end{enumerate}
\end{assumption}
The approximate maximiser condition at rate $n^{-1}$ is mild and can be ensured
by the econometrician during numerical optimisation.
Condition \eqref{asm--X-ac-smooth-m-diff-p-smooth} imposes regularity on the conditional distribution of the covariates by requiring absolute continuity of \(X \mid \mathcal{U}_{\bm{e}}\) and continuous differentiability of the associated density \(p_{\bm{e}}(x \mid \bm{u}_{\bm{e}})\) in \(x\), which ensures the absence of point masses and enables local smooth approximations of the population objective. Condition \eqref{asm--X-ac-smooth-m-diff-m-smooth} complements this by requiring that the conditional regression function \(m_{0,\bm{e}}(x,\bm{u}_{\bm{e}})\) is continuously differentiable in \(x\), with the derivative \(\dot{m}_{0,\bm{e}}(x,\bm{u}_{\bm{e}})\) characterising the local curvature of the population objective. Condition \eqref{asm--X-ac-smooth-m-diff-m-dominance} introduces a square-integrable envelope \(m_\ast(x,\bm{u}_{\bm{e}})\) that uniformly dominates the density, regression function, and its derivative, thereby ensuring the validity of dominated convergence arguments. Finally, Assumption \eqref{asm--X-ac-smooth-m-diff-m-increasing} imposes a strong median identification condition: the normalisation \(m_{0,{\bm{e}}_k}(x,u)=0\) when \(x'\beta_0=0\) centres the decision boundary, while strict monotonicity of \(t \mapsto m_{0,{\bm{e}}_k}(x+t\beta_0,u)\) guarantees a unique sign change across the hyperplane \(x'\beta_0=0\), ruling out flat regions and ensuring point identification together with the local curvature required for quadratic expansion.
The following result establishes asymptotic distributional theory for the
maximum score estimator under the local parametrisation \eqref{eqn--beta-local}.
\begin{theorem}
\label{thm--max-score-asymp-normal}
Let \Cref{asm--exchangeable,asm--max-score-consistency,asm--X-ac-smooth-m-diff}
hold, and suppose that \(n / N_{k} \to \lambda_{k} \in [0, \infty)\) for each
\(k \in \{1, \dots, K\}\).
Then there is a positive semi-definite matrix \(V\) such that
\begin{align*}
n^{1 / 2} \widehat \theta_{{\bm{N}}} \leadsto \mathrm{N} (0, V).
\end{align*}
\end{theorem}
\begin{proof}[Proof of \Cref{thm--max-score-asymp-normal}]
A proof can be found in \Cref{sec:proof--thm-max-score-asymp-normal} of the
appendix.
\end{proof}
The proof proceeds by verifying the conditions of \Cref{thm--non-smooth-M-Estimation} in \Cref{sec:M-estimation} for the maximum score estimator. In particular, establishing stochastic differentiability and local stochastic equicontinuity is highly non-trivial. The former draws on analytic techniques for maximum score developed in \cite{kim1990cube}, while the latter relies on a localisation argument for empirical processes based on a local maximal inequality established in \cite{chen2026cross}.
An inspection of the proof of \Cref{thm--max-score-asymp-normal} reveals that the asymptotic variance possesses a complex structure and depends on various unknown quantities; consequently, obtaining a consistent variance estimator is challenging. Fortunately,
the bootstrap procedure described in \Cref{thm:bootstrap} of \Cref{sec:M-estimation} provides a simple and valid framework for statistical inference regarding the maximum score estimator under multiway dependence. This validity is guaranteed as the assumptions of \Cref{thm--max-score-asymp-normal} sufficiently satisfy all underlying requirements.
\section{M-Estimation with Multiway Dependence}\label{sec:M-estimation}
In this section, we develop a general asymptotic theory for non-smooth M-estimators under multiway dependence as well as a bootstrap procedure for statistical inference.
Let $\mathcal{F}=\{f_\theta:\mathcal{S}\to \mathbb{R}:\theta\in \Theta\subset\mathbb{R}^d,\; \mathbb{\mathbb E}|f_\theta (W)|<\infty\}$, where $\Theta$ contains a neighbourhood of the origin and $\theta\mapsto f_\theta $ may be non-smooth.
We define the population objective function for the M-estimation problem by
\(
Q(\theta)=\mathbb{\mathbb E}[f_\theta(W)], \; \theta\in\Theta,
\)
and, without loss of generality, normalise the true parameter to $\theta_0$ to the origin.
Suppose we observe random variables
\(
\mathcal D_{\bm{N}}=\{ W_{\bm{i}} :{\bm{i}} \in I_{\bm{N}} \},
\)
an estimator based on $\mathcal D_{\bm{N}}$ can now be defined as
\begin{align*}
\widehat \theta_{{\bm{N}}}=\underset{\theta\in\Theta}{\text{argmax}}\:\mathbb{\mathbb E}_{{\bm{N}}} f_\theta,
\end{align*}
or more generally
\begin{align*}
\sup_{\theta\in \Theta}\mathbb{\mathbb E}_{{\bm{N}}} f_\theta-\mathbb{\mathbb E}_{{\bm{N}}} f_{\widehat\theta_{{\bm{N}}}}=o_P(n^{-1}).
\end{align*}
\begin{comment}
We assume that the data random variables $\{W_{\bm{i}}\}_{\bm{i}}$ are defined on a probability space $(S,\mathcal{S},P)$ and admits the Aldous--Hoover--Kallenberg (AHK) representation
\begin{align}
W_{\bm{i}} = \tau\Bigl( \{ U_{\bm{i} \odot \bm{e}} \}_{\bm{e} \in \{0,1\}^K \setminus \{\bm{0}\}} \Bigr), \label{eq:AHK_representation}
\end{align}
where $\odot$ denotes the Hadamard product, the latent variables
\(
\{ U_{\bm{i} \odot \bm{e}} \}
\)
are mutually independent and identically distributed on $(0,1)$, and $\tau$ is Borel measurable.
\end{comment}
Before stating our assumptions, let us introduce the Hoeffding-type decomposition \citep{hoeffding1948class,chiang2023inference} for multiway dependence.
For an $f:\mathcal{S}\to \mathbb{R}$ with $\mathbb{\mathbb E} [f(W_{\bm{i}})]=0$\footnote{Otherwise one may replace $ f$ with $\widetilde f(W_{\bm{i}})=f(W_{\bm{i}})-\mathbb{\mathbb E}[f(W_{\bm{i}})]$.} and any $\bm{i}\in I_{\bm{N}}$, define the conditional expectation
\[
(P_{\bm{e}}f)(\{U_{\bm{i}\odot \bm{e}'}\}_{\bm{e}'\le \bm{e}})
=\mathbb{E}\!\left[f(W_{\bm{i}})\mid\{U_{\bm{i}\odot \bm{e}'}\}_{\bm{e}'\le \bm{e}}\right].
\]
The orthogonal projections $\pi_{\bm{e}}$ are then defined recursively: for ${\bm{e}}_k\in \mathcal{E}_1$,
\[
(\pi_{\bm{e}_k}f)(U_{\bm{i}\odot \bm{e}_k})
=(P_{\bm{e}_k}f)(U_{\bm{i}\odot \bm{e}_k}),
\]
and for $\bm{e}\in\cup_{k=2}^K\mathcal{E}_k$,
\[
(\pi_{\bm{e}}f)\bigl(\{U_{\bm{i}\odot \bm{e}'}\}_{\bm{e}'\le \bm{e}}\bigr)
=(P_{\bm{e}}f)\bigl(\{U_{\bm{i}\odot \bm{e}'}\}_{\bm{e}'\le \bm{e}}\bigr)
-\sum_{\bm{e}'\le \bm{e},\, \bm{e}'\neq \bm{e}}
(\pi_{\bm{e}'}f)\bigl(\{U_{\bm{i}\odot \bm{e}''}\}_{\bm{e}''\le \bm{e}'}\bigr).
\]
Following Lemma~1 in \cite{chiang2023inference}, for any \(\ell\in \mathrm{supp}(\bm{e})\)\footnote{That is, $\mathrm{supp}({\bm{e}}) = \{ j=1,\dots,K : e_j \ne 0\}$.} the random variable
\(
(\pi_{\bm{e}}f)\bigl(\{U_{\bm{i}\odot \bm{e}'}\}_{\bm{e}'\le \bm{e}}\bigr)
\)
is centred conditionally on
\(
\{U_{\bm{i}\odot \bm{e}'}\}_{\bm{e}'\le \bm{e}-\bm{e}_{\ell}}.
\)
Define
\(
I_{\bm{N},\bm{e}} = \{\bm{i}\odot \bm{e} : \bm{i}\in I_{\bm{N}}\},
\)
so that
\(
|I_{\bm{N},\bm{e}}| = \prod_{k'\in\mathrm{supp}(\bm{e})} N_{k'}.
\)
Accordingly, define the ${\bm{e}}$-specific Hoeffding-type projection by
\[
H_{\bm{N}}^{\bm{e}}(f)
=\frac{1}{|I_{\bm{N},\bm{e}}|}
\sum_{\bm{i}\in I_{\bm{N},\bm{e}}}
(\pi_{\bm{e}}f)\bigl(\{U_{\bm{i}\odot \bm{e}'}\}_{\bm{e}'\le \bm{e}}\bigr).
\]
We then obtain the Hoeffding-type decomposition
\begin{align}
\mathbb{E}_N f
= \sum_{k=1}^K \sum_{\bm{e}\in\mathcal{E}_k} H_{\bm{N}}^{\bm{e}}(f).
\label{eq:hoeffding}
\end{align}
We impose the following conditions for asymptotic theory of general non-smooth M-estimators.
\begin{assumption}[Non-smooth M-estimation]
\label{asm--M-estimation-general}
Assume the following hold:
\begin{enumerate}[(i)]
\item \label{asm--M-estimation-general-id}
(Identification) $Q$ has a unique maximum attained at $\theta=0$
and the following identity holds
\begin{align}
Q(\theta)=Q(0)-\frac{1}{2}\theta'H\theta+o(\|\theta\|^2),
\label{eq:a_quadratic_approximation}
\end{align}
where $H$ is a symmetric and positive definite $d\times d$ matrix.
\item \label{asm--M-estimation-general-consistency}
(Consistency) Let $\{\widehat \theta_{{\bm{N}}} \}$ be such that
\begin{align}
\widehat \theta_{{\bm{N}}} =o_P(1)\label{eq:a_consistency}
\end{align}
and it approximates the maximiser of the sample counterpart of $Q$ in the sense that
\begin{align}
\sup_{\theta\in \Theta}\mathbb{\mathbb E}_{{\bm{N}}} f_\theta-\mathbb{\mathbb E}_{{\bm{N}}} f_{\widehat\theta_{{\bm{N}}}}=o_P(n^{-1}).\label{eq:a_solution}
\end{align}
\item \label{asm--M-estimation-general-stoch-diff}
(Stochastic Differentiability)
For each unit vector ${\bm{e}}_k\in \mathcal{E}_1$, there
exists a measurable $\mathbb{R}^d$-valued function $\Delta_k=\Delta_k (U_{{\bm{e}}_k})$ satisfying $\mathbb{\mathbb E} \Delta_k=0$ and $\mathbb{\mathbb E} \|\Delta_k\|^2<\infty$, and $r_k=r_k(U_{{\bm{e}}_k};\theta)$ such that
$r_k(u;0)=0$ for all $u\in (0,1)$,
\begin{align*}
r_k(u;\theta)=\frac{\pi_{{\bm{e}}_k}(f_\theta-\mathbb{\mathbb E}[f_\theta])
(u)-\pi_{{\bm{e}}_k}(f_0-\mathbb{\mathbb E}[f_0] )(u)-\theta'\Delta_k(u)}{ \|\theta\| }
\end{align*}
for all $\theta\in \Theta\setminus \{0\}$, $u\in (0,1)$, and
\begin{align}
\sup_{\|\theta\|\le
\delta}\frac{|\sqrt{N_{k}} (\mathbb{\mathbb E}_{k, {\bm{N}}}-\mathbb{\mathbb E}) r_k(\cdot;\theta)|}{1 +
n^{1/2} \|\theta\|}=o_P(1),
\label{eq:a_differentiablity}
\end{align}
where $\mathbb \mathbb{\mathbb E}_{k, {\bm{N}}}=N_k^{-1}\sum_{i_k=1}^{N_k}\delta_{i_k}$ is the sample average over $i_k=1,...,N_k$.
\item \label{asm--M-estimation-general-stoch-eq}
(Local Stochastic Equicontinuity)
The class of functions $\mathcal{F}$ is of VC-type with characteristics $A\ge (e^{2(K-1)}/16)\vee e$ and $v\ge 1$ and an envelope $F$ with $\mathbb{\mathbb E}[F^2]<\infty$. Furthermore, for ${\bm{e}}\in\mathcal{E}_2$, and $\delta_n=O(n^{-1/2})$, it holds that
\begin{align}
\sup_{\|\theta\|\le \delta_n}\left|H_{{\bm{N}}}^{\bm{e}}(f_\theta - f_0) \right|=o_P(n^{-1})\label{eq:a_stochastic_equicontinuity}
\end{align}
\end{enumerate}
\end{assumption}
Assumption \ref{asm--M-estimation-general} (i) imposes a standard identification condition, together with a quadratic expansion of the population objective, and normalises the unique maximiser to the origin. Part (ii) requires the existence and consistency of an estimator $\widehat\theta_{{\bm{N}}}$, which can typically be established using conventional M-estimation arguments; see Section 2 of \citet{newey1994large}. Part (iii) imposes a mild stochastic differentiability condition on the first-order projections. Importantly, it accommodates a broad class of non-smooth objective functions; see \citet{pollard1985new} for details and sufficient conditions. Finally, part (iv) imposes a local stochastic equicontinuity condition on the localised second-order projections of the objective function, which can typically be verified using an appropriate maximal inequality. This requirement is weaker than the comparable stochastic equicontinuity condition employed in $U$-process results, such as Theorem 1 of \citet{arcones1994estimators}, as it only requires \ref{eq:a_stochastic_equicontinuity} to hold for $\delta_n = O(n^{-1/2})$, rather than for all sequences $\delta_n \to 0$.
The following provides asymptotic theory for general non-smooth M-estimators with multiway dependence.
\begin{theorem}[Non-smooth M-estimation]\label{thm--non-smooth-M-Estimation}
Suppose \Cref{asm--exchangeable,asm--M-estimation-general} hold, and $n/N_k\to \lambda_k\in[0,\infty)$ for each $k=1,...,K$, then
\begin{align*}
n^{1/2}\widehat \theta_{{\bm{N}}} = n^{1/2}H^{-1}\overline \psi_{\bm{N}}+o_P(1)\leadsto N(0,V),
\end{align*}
where $\overline \psi_{\bm{N}}=\sum_{k=1}^K \mathbb \mathbb \mathbb E_{k,{\bm{N}}}\Delta_k$, $V=H^{-1}\Omega H^{-1}$, and $\Omega=\sum_{k=1}^K \lambda_k\mathbb{\mathbb E}[\Delta_k(U_{\bm{1}\odot {\bm{e}}_k})\Delta_k(U_{\bm{1}\odot {\bm{e}}_k})']$.
\end{theorem}
\begin{proof}[Proof of \Cref{thm--non-smooth-M-Estimation}]
A proof can be found in \Cref{sec:proof--thm--non-smooth-M-estimation} of the
appendix.
\end{proof}
Note that the variance can be zero -- in this case, the limiting Gaussian
distribution degenerates at zero.
Our proof strategy differs from existing approaches, such as \citet{arcones1994estimators}, in that we do not rely on stochastic equicontinuity. In the present multiway setting, an analogue of the weak convergence result for degenerate
$U$-processes (cf.\ Corollary 5.7 of \citet{arcones1993limit}) is not available, and hence the higher-order projection terms must be controlled by alternative means. To this end, we adopt an approach combining localised empirical process methods with an iterated argument that first delivers an \(O_P(n^{-1})\) rate for the higher-order projection terms and then sharpens it to \(o_P(n^{-1})\), using Condition \eqref{eq:a_stochastic_equicontinuity}.
In non-smooth M-estimation problems, estimation of the asymptotic variance is typically challenging. The components $\Delta_k$ and $H$ involve unknown derivatives of population objects and, because of the lack of smoothness of $\theta\mapsto f_\theta$, generally do not possess obvious feasible sample counterparts. Consequently, direct variance estimation is challenging. Fortunately, the linear representation and Gaussian approximation established in Theorem \ref{thm--non-smooth-M-Estimation} enable us to conduct inference via a bootstrap procedure. The next theorem states this result.
\begin{theorem}[Bootstrapping M-estimators]\label{thm:bootstrap}
Suppose the conditions of Theorem \ref{thm--non-smooth-M-Estimation} hold and the asymptotic variance $V$ is positive definite. For each $k=1,\ldots,K$, let
\(
\xi^k=(\xi_1^k,\ldots,\xi_{N_k}^k)
\)
consist of i.i.d. random variables with $\mathbb{\mathbb E}[\xi_{i_k}^k]=\operatorname{\text{Var}}(\xi_{i_k}^k)=1$, $\mathbb{\mathbb E}[|\xi_{i_k}^k|^3]<\infty$, and $\xi^1,\ldots,\xi^K$ are mutually independent and independent of $\mathcal D_{\bm{N}}$. Define
\(
\xi_{\bm{i}}=\prod_{k=1}^K \xi_{i_k}^k
\)
and \(
\mathbb{\mathbb E}_{{\bm{N}}} \xi f_\theta=|I_{\bm{N}}|^{-1}\sum_{{\bm{i}}\in I_{\bm{N}}}\xi_{\bm{i}} f_\theta(W_{\bm{i}}).
\)
Let $\widehat\theta_{\bm{N}}^*$ denote a bootstrap $M$-estimator
that satisfies
\[\sup_{\theta\in\Theta} \mathbb{\mathbb E}_{{\bm{N}}} \xi f_\theta-\mathbb{\mathbb E}_{{\bm{N}}} \xi f_{\widehat \theta_{{\bm{N}}}^* }=o_P(n^{-1}).\] Then, as $n\to \infty$,
we have
\[
n^{1/2}(\widehat\theta_{\bm{N}}^*-\widehat\theta_{\bm{N}})\overset{*}{\leadsto }N(0,V)
\]
conditionally on $\mathcal D_{\bm{N}}$
with probability approaching one.
\end{theorem}
\begin{proof}[Proof of \Cref{thm:bootstrap}]
A proof can be found in \Cref{sec:proof_thm:bootstrap} of the appendix.
\end{proof}
The restrictions imposed on the weights $\xi^k$ are met, for instance, when $\xi^k$ is drawn from either an $\text{Exponential}(1)$ or a $\text{Poisson}(1)$ distribution. This bootstrap is closely related to the pigeonhole bootstrap considered by \citet{owen2007pigeonhole} and \citet{2021daveziesEmpiricalProcessResults}.
Unlike the pigeonhole bootstrap, which relies on multinomial weights, the proposed method assigns i.i.d. weights to each $k$.
The proof proceeds by introducing an alternative Hoeffding-type decomposition for the bootstrapped process on an expanded probability space that incorporates the bootstrap weights. We employ i.i.d.\ weights for each $k$, which is convenient for deriving this decomposition. We then show that the bootstrap estimator admits an analogous unconditional asymptotic linear representation with the corresponding multiplicative i.i.d.\ weights, a property that underpins the validity of our bootstrap procedure and may be of independent interest.
\section{Conclusion}
In summary, this paper develops a unified asymptotic framework for non-smooth M-estimators under multiway dependence and applies it to the maximum score estimator, establishing asymptotic Gaussianity at the parametric rate together with a valid bootstrap procedure. These results highlight that complex dependence structures can fundamentally alter the inferential properties of non-smooth estimators, rendering standard difficulties under i.i.d. settings tractable in this context.
\bigskip