Gaussian approximation for maximum score and non-smooth M-estimators with multiway dependence

Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.

37,390 characters · 5 sections · 43 citation commands

Rendered from LaTeX for readability, not typeset faithfully. Citation keys are highlighted; maths is left as source; figures, tables and equation environments are summarised rather than reproduced; unrecognised commands are greyed out so nothing is silently dropped. Email addresses are removed.

Gaussian Approximation for Maximum Score and Non-Smooth M-Estimators with Multiway Dependence

abstractThe maximum score estimator of 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.

Introduction

The maximum score estimator of 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, 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 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, 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 miglioretti2007marginal,petersen2008estimating,cameron2011robust,thompson2011simple,cameron2015practitioner,mackinnon2023cluster—the original maximum score estimator of 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, menzel2021bootstrap, 2021daveziesEmpiricalProcessResults,davezies2025analytic, chiang2023inference, and 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 liu1990notion and Oja's spatial medians oja1983descriptive, established in arcones1994estimators, as well as the notion of data-generating-process (DGP)–induced smoothing studied in panel quantile regression with common shocks by 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 pollard1985new. A further complication arises in establishing stochastic equicontinuity for higher-order projections. In the classical $U$-process settings of arcones1994estimators, such results rely on weak convergence for degenerate processes 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 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, 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, 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, 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 sherman1993limiting. This reformulation allows them to establish asymptotic normality under suitable regularity conditions. Closely related is 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 manski1975maximum can, in fact, be asymptotically Gaussian under conditions of complex multiway dependence.

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

Maximum Score Estimator

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

equation[equation omitted — 336 chars of source]

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 menzel2021bootstrap and 2021daveziesEmpiricalProcessResults,davezies2025analytic, among others.

assumption$\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), \end{align} where $\odot$ denotes the Hadamard product.

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

equation[equation omitted — 79 chars of source]

We impose the following modelling and regularity conditions on the joint distributions of $(Y,X')$.

assumption\begin{enumerate}[(i)] • $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. • 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$. • The estimation set $\mathcal{B}$ in (ref) 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}

Lower-level sufficient conditions for (ref) (ref) and (ref) are well known in the literature. For example 1985manskiSemiparametricAnalysisDiscrete assumes an index restriction coupled with a special regressor. In particular, Assumption 2 of 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 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 (ref) (ref) and (ref). Whilst the assumptions of 1985manskiSemiparametricAnalysisDiscrete allow for some discrete (or constant) regressors, 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 (ref) (ref) with $\mathcal{B} = \mathbbm{S}^{d - 1}$.

The following results provides consistency of the maximum score estimator under multiway dependence.

lemmaSuppose (ref) hold, and further assume the restricted estimation set $\mathcal{B}$ in (ref) satisfies $\beta_{0} \in \mathcal{B}$ for $\beta_{0}$ in (ref) (ref). Then $\widehat{\beta}_{n}$ in (ref) satisfies $\widehat{\beta}_{n} \overset{\mathrm{p}}{\to} \beta_{0}$.

proof[Proof of (ref)] See (ref).

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:

equation*[equation* omitted — 173 chars of source]

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

equation[equation omitted — 118 chars of source]

The parametrisation in (ref) 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

equation*[equation* omitted — 164 chars of source]

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

equation*[equation* omitted — 215 chars of source]

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

assumptionThe 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)] • 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$. • 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}})$. • 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*} • 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}

The approximate maximiser condition at rate $n^{-1}$ is mild and can be ensured by the econometrician during numerical optimisation. Condition (ref) 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 (ref) 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 (ref) 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 (ref) 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 (ref).

theoremLet (ref) 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*}

proof[Proof of (ref)] A proof can be found in (ref) of the appendix.

The proof proceeds by verifying the conditions of (ref) in (ref) 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 kim1990cube, while the latter relies on a localisation argument for empirical processes based on a local maximal inequality established in chen2026cross.

An inspection of the proof of (ref) 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 (ref) of (ref) 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 (ref) sufficiently satisfy all underlying requirements.

M-Estimation with Multiway Dependence

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

align*[align* omitted — 127 chars of source]

or more generally

align*[align* omitted — 157 chars of source]

commentWe 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), \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.

Before stating our assumptions, let us introduce the Hoeffding-type decomposition 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 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

align[align omitted — 123 chars of source]

We impose the following conditions for asymptotic theory of general non-smooth M-estimators.

assumption[Non-smooth M-estimation] Assume the following hold: \begin{enumerate}[(i)] • (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), \end{align} where $H$ is a symmetric and positive definite $d\times d$ matrix. • (Consistency) Let $\{\widehat \theta_{{\bm{N}}} \}$ be such that \begin{align} \widehat \theta_{{\bm{N}}} =o_P(1) \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}). \end{align} • (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), \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$. • (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}) \end{align} \end{enumerate}

Assumption (ref) (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 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 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 arcones1994estimators, as it only requires (ref) 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.

theorem[Non-smooth M-estimation] Suppose (ref) 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})']$.

proof[Proof of (ref)] A proof can be found in (ref) of the appendix.

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

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) enable us to conduct inference via a bootstrap procedure. The next theorem states this result.

theorem[Bootstrapping M-estimators] Suppose the conditions of Theorem (ref) 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.

proof[Proof of (ref)] A proof can be found in (ref) of the appendix.

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

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.