Cross-Fitting-Free Debiased Machine Learning 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.

57,416 characters · 9 sections · 52 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.

\address{ School of Economics, Shanghai University of Finance and Economics, 111 Wuchuan Road, Yangpu District, Shanghai 200434, China} \email{[email removed]}

\address{Department of Economics, University of Wisconsin-Madison, 1180 Observatory Drive Madison, WI 53706-1393, USA.} \email{[email removed]}

abstractThis paper develops an asymptotic theory for two-step debiased machine learning (DML) estimators in generalised method of moments (GMM) models with general multiway clustered dependence, without relying on cross-fitting. While cross-fitting is commonly employed, it can be statistically inefficient and computationally burdensome when first-stage learners are complex and the effective sample size is governed by the number of independent clusters. We show that valid inference can be achieved without sample splitting by combining Neyman-orthogonal moment conditions with a localisation-based empirical process approach, allowing for an arbitrary number of clustering dimensions. The resulting debiased GMM estimators are shown to be asymptotically linear and asymptotically normal under multiway clustered dependence. A central technical contribution of the paper is the derivation of novel global and local maximal inequalities for general classes of functions of sums of separately exchangeable arrays, which underpin our theoretical arguments and are of independent interest.

Introduction

The debiased machine learning (DML), also known as the double/debiased machine learning framework, has become a leading approach to two-step estimation with high-dimensional or nonparametric nuisance components; see, among many others, chernozhukov2018double, chernozhukov2022, and escanciano2022debiased. A central feature of this literature is cross-fitting, whereby sample splitting is used to separate nuisance estimation from target parameter evaluation. This device is motivated by two considerations. First, the conventional view holds that it mitigates overfitting bias arising from the use of highly flexible first-stage learners. Second, it relaxes empirical process conditions by reducing the dependence between first-stage estimation errors and second-stage score evaluation. As a result, cross-fitting has become close to a default recommendation in both theoretical analyses and empirical implementations of DML procedures.

At the same time, many empirical applications in economics and finance involve multiway clustered dependence; see, for example, petersen2008estimating,cameron2015practitioner. In such settings, observations may be correlated along multiple dimensions, such as firm and time or region and industry. Extensions of DML to multiway clustered environments are developed in chiang2022multiway.

However, extensive sample splitting is not without cost. First, with a finite number of folds, cross-fitting yields a random estimator, which may hinder reproducibility. Increasing the number of folds can mitigate this concern to some extent, but at the expense of substantially greater computational burden, particularly when the first stage is complex and/or requires extensive tuning. Second, cross-fitting effectively reduces the sample size available to each first-stage problem. This is especially consequential when these stages involve high-dimensional or nonparametric estimation, where performance is inherently variance-sensitive. The resulting loss in nuisance estimation precision may, in finite samples, translate into non-negligible efficiency losses for the parameter of interest.

These concerns are further amplified under multiway clustering, where the effective sample size is determined by the number of independent cluster units rather than the total number of observations. Partitioning the data into folds may therefore leave each subsample with only a limited number of independent clusters, making cross-fitting particularly costly, in addition to increasing computational burden. Moreover, in two-step procedures such as double machine learning, overfitting bias arising from first-stage nuisance estimation need not be intrinsically detrimental for inference on the target parameter, in contrast to classical one-step settings. This suggests that the conventional rationale for cross-fitting may be less central than is sometimes presumed. Indeed, even under i.i.d. settings, the literature provides little general theoretical support for the advantages of cross-fitting, apart from a few special cases considered, for example, by newey2018cross.

Motivated by these considerations, a growing literature seeks to weaken or dispense with cross-fitting in general nonlinear estimation problems. One approach to obtaining theoretical guarantees for DML without cross-fitting is the “localisation" method, as employed for example in belloni2015uniform,belloni2018uniformly. The key idea is to analyse the supremum of an empirical process indexed by estimating equations evaluated over deterministic sets that localise the nuisance parameter around its true value. These sets are constructed so that the nuisance estimator lies in them with probability approaching one, thereby replacing a stochastic index with a deterministic one and disentangling the dependence between estimating equations and first-stage estimates. If the localisation sets shrink at appropriate rates—typically verified through maximal inequalities—the associated empirical process remainder is of smaller order than the leading asymptotically linear term and does not affect the limiting distribution. An alternative strategy is based on a “stability" conditions; see chen2022debiased under i.i.d. and a generalisation in cao2025neighborhood under spatial/network $\beta$-mixing. Verifying such conditions, however, is often substantially more involved beyond certain well-understood cases, and we therefore do not pursue this route.

In this paper, we develop a general asymptotic theory for two-step debiased GMM estimators under multiway clustered dependence, without relying on cross-fitting. We consider a broad class of locally robust two-step GMM problems similar to those studied in chernozhukov2022 in which low-dimensional target parameters are identified by orthogonal moment conditions that depend on high-dimensional or nonparametric nuisance components. The data are allowed to exhibit dependence along an arbitrary number of clustering dimensions, accommodating empirical settings in which correlation arises simultaneously across, for example, firms, time periods, locations, or networks. Our results establish asymptotic linearity and asymptotic normality of the resulting estimators under conditions that permit flexible highly first-stage learners while avoiding sample splitting. The analysis explicitly accounts for the reduced effective sample size induced by multiway clustering and provides inference procedures that remain valid when the number of independent cluster units, rather than the total number of observations, governs the stochastic order. By combining orthogonality with a localisation-based empirical process argument tailored to multiway clustered arrays, we show that the impact of first-stage estimation can be controlled without cross-fitting. This yields a unified framework for debiased GMM inference that is well suited to empirically relevant clustered environments where conventional cross-fitting can be statistically and computationally costly.

Maximal inequalities are central to localisation-based arguments in DML, yet for multiway clustered—more precisely, separately exchangeable (SE)—arrays, the available theory remains limited. In particular, no general global maximal inequality accommodates arbitrary numbers of clustering dimensions $K$, arbitrary moments $q$, and infinite pointwise measurable function classes. Existing results address only special cases: Theorem B.2 of chiang2023inference allows general $K$ and $q\in[1,\infty)$ but restricts attention to finite classes, while Lemma C.3 of liu2024estimation covers general classes only for $q=1$ and $K=2$. The situation is even more restrictive for local maximal inequalities, which are essential for sharper convergence rates: unlike in the related $U$-statistics literature (see ChenKato2019b), beyond the $K=1$ (i.i.d.) case, no such result appears to be available for SE arrays. These gaps reflect the intrinsic difficulty posed by multiway dependence, which generates complex interactions across observations and undermines classical tools such as Hoeffding averaging. To overcome this, we develop a new proof strategy based on a transversal partition of the index set that effectively decouples dependence and permits the use of the Hoffmann--J{\o}rgensen inequality, yielding sharp higher-moment bounds. This approach delivers both global and local maximal inequalities for potentially uncountable, pointwise measurable function classes under SE sampling.

The paper is organised as follows. Section 2 introduces the multiway clustered sampling setups and notations. Section 3 develops the cross-fitting-free debiased GMM estimator and presents the main asymptotic results under general high-level conditions, while also providing three examples for the rate and complexity conditions and validity of variance estimation. Section 4 provides the core technical tools in the form of new maximal inequalities for empirical processes for separately exchangeable arrays. Proofs and supplementary arguments are collected in the appendices.

\vskip 0.15in

Setups and Notations

In this section, we introduce the framework of multiway clustered sums that will be used throughout the paper. Let $K$ be a fixed positive integer, and denote a $K$-tuple index by $ \bm{i} = (i_1, i_2, \dots, i_K) \in \mathbb{N}^K. $ Suppose we observe a $K$-way array of data \[ \{ X_i : \bm{i} \in [N_1] \times \cdots \times [N_K] \}, \] where $[N_k] := \{1,\dots,N_k\}$ denotes the index set of dimension $k$ and $N_k$ is the sample size. Let $\bm{N} = (N_1, N_2, \dots, N_K)$ and define $ [\bm{N}] = \prod_{k=1}^{K} \{1,2,\dots,N_k\}. $ We also denote \[ N = \prod_{k=1}^K N_k,\quad n = \min\{N_1, N_2, \dots, N_K\},\quad \text{and}\quad \overline{N} = \max\{N_1,...,N_K\}. \]

Suppose $\{ X_i\}_i$ are random variables defined on a probability space $(S,\mathcal{S}, P)$\footnote{Our framework accommodates high-dimensional regimes in which the array of data-generating processes may depend on the sample sizes; for notational economy, this dependence is left implicit. }, and $\{ X_i\}$ satisfy the separate exchangeability (SE) and dissociation (D) conditions defined below.

enumerate• For any $\pi =(\pi_1,\dots,\pi_K)$, a $K$-tuple of permutations of $ N$, $\{X_{\bm{i}}\}_{\bm i \in [\bm N]}$ and $\{X_{\pi(\bm{i})}\}_{\bm i \in [\bm N]}$ are identically distributed. • For any two disjoint sets of indices $I,I'\subset \mathbb{N}^K$, $\{X_{\bm{i}}\}_{{\bm{i}}\in I}$ and $\{X_{\bm{i}}\}_{{\bm{i}}\in I'}$ are independent.

Under Conditions (SE) and (D), the Aldous--Hoover--Kallenberg (AHK) representation (see Corollary 7.35 in kallenberg2005probabilistic) guarantees the existence of the following representation:

align[align omitted — 149 chars of source]

where $\odot$ denotes the Hadamard (element-wise) product, the collection $ \{ U_{\bm{i} \odot \bm{e}} : \bm{i} \in \mathbb{N}^K,\; \bm{e} \in \{0,1\}^K \setminus \{\bm{0}\} \} $ consists of mutually independent and identically distributed (i.i.d.) random variables, and $\tau$ is a Borel measurable map taking values in $\mathcal{S}$.

We say a class of functions $\mathcal{F}:\mathcal{S}\to \mathbb{R}$ is pointwise measurable if there exists a countable subclass $\mathcal{F}'\subset \mathcal{F}$ such that for each $f\in \mathcal{F}$, there exists a sequence $(f_j)_j\subset \mathcal{F}'$ such that $f_j\to f$ pointwisely. Given the observed set of random variables $\{X_{\bm{i}}:{\bm{i}} \in [\bm{N}]\}$ that satisfy Conditions (SE) and (D), and a pointwise measurable class of functions $\mathcal{F}$ with elements $f: \mathcal{S} \to \mathbb{R}$, define the sample mean process by $ \mathbb{E}_N f = N^{-1}\sum_{\bm{i} \in [\bm{N}]} f(X_{\bm{i}}) $ and, suppose $P|f|:=\int |f| dP<\infty$\footnote{See Section 2.3 in pollard2002user for detailed explanation of empirical process notation for measure and integration.}, the empirical process by \[ \mathbb{G}_{n}(f) = \frac{\sqrt{n}}{N} \sum_{\bm{i} \in [\bm{N}]} \Bigl\{ f(X_{\bm{i}}) - P(f)\bigr] \Bigr\}. \] Throughout, we use the shorthand $\mathbb E[f(X_{\bm{1}})]= P(f)$, where $\bm{1}=(1,...,1)$, for expectations taken with respect to the data-generating distribution.

Notation.

Let $\mathbb{N}$ denote the set of positive integers and $\mathbb{R}$ for the real line. For $a,b\in \mathbb{R}$, let $a\vee b=\max\{a,b\}$ and $a\wedge b = \min\{a,b\}$. Denote for $m\in \mathbb N$ that $ [m] = \{1,2,\ldots,m\}. $ For real vectors $\bm{a}= (a_{1},\dots,a_{K})$ and $\bm{b} = (b_{1},\dots,b_{K})$, we denote $\bm{a} \le \bm{b}$ for $a_{j} \le b_{j}$ for all $1 \le j \le K$. Let $\mathrm{supp}(\bm{a}) = \{ j : a_j \ne 0\}$. We denote by $\odot$ the Hadamard product: for ${\bm{i}} = (i_1,\dots,i_K)$ and $\bm{j} = (j_1,\dots,j_K)$, ${\bm{i}} \odot \bm{j}= (i_1 j_1,\dots,i_K j_K)$. For each $k=0,1,2,...,K$, define $ \mathcal{E}_k=\{{\bm{e}}\in \{0,1\}^K: \Vert{\bm{e}} \Vert_0 =k\}$ and thus $\{0,1\}^K=\cup_{k=0}^K \mathcal{E}_k$. For $q\in [1,\infty]$, let $\|f\|_{Q,q}=( Q|f|^q)^{1/q}$. For a non-empty set $T$ and $f:T\to \mathbb{R}$, denote $\|f\|_T=\sup_{t\in T}|f(t)|$. For a pseudometric space $(T,d)$, let $N(T,d,\varepsilon)$ denote the $\varepsilon$-covering number for $(T,d)$. We say $F:\mathcal{S}\to \mathbb{R}_+$ is an envelope for a class of functions $\mathcal{F}\ni f:\mathcal{S}\to \mathbb{R}$ if $\sup_{f\in\mathcal{F}}|f(x)|\le F(x)$ for all $x\in\mathcal{S}$.

\vskip 0.15in

Debiased Machine Learning for GMM

In this section, we DML type two-step estimation and inference approaches in a setting where data is multiway clustered. Particularly, it is of practical importance to study debiased machine learning without sample-splitting due to the poor usage of samples when splitting the data. While cross-fitting can improve the sample usage by switching the roles of split samples, the improvement is limited in a multiway clustered setting where cross-fitting is done in a way that more data is excluded when estimating the high-dimensional nuisance parameters.

Since most econometric and statistical models can be reduced to moment restrictions, we consider DML without sample splitting in a GMM setup similar to those considered in chernozhukov2022:

enumerate• The parameters $\theta_0 \in \Theta \subset R^d$ of interest are fixed-dimensional and (over-)identified by a set of moment conditions that depends on high-dimensional nuisance parameters: $g: \mathcal{X}\times\Theta\times {\Gamma}\rightarrow\mathbb{R}^q$ such that \begin{equation} \mathbb{E}\left[ g(X,\theta_0,\gamma_0) \right] = 0, \end{equation} • The nuisance parameters $\gamma_0\in \Gamma$ are exactly identified and estimated by machine learners using the full sample. • With the full sample again, solving the GMM problem using an orthogonalized moment function and the corresponding optimal weighting matrix, with plug-in nuisance estimates.

DML for two-step GMM with i.i.d. data is studied in chernozhukov2022. In contrast, our setting involves multiway-clustered sampling, which leads to substantially different asymptotic arguments. Moreover, our theoretical framework eliminates the need for cross-fitting.

Main results for the debiased GMM estimator.

The first component of DML is the orthogonalisation of the moment condition. Specifically, we construct $\psi(X,\theta_0,\eta_0)$ by adding an adjustment term (which may depends on extra nuisance parameters contained in $\eta_0\in \Gamma$) to $g$ such that $\psi$ is mean zero and the path-wise derivative with respect to $\eta$ in the direction $\widetilde\eta\in \Gamma$ is zero (or vanishing) when evaluated at the truth:

align[align omitted — 198 chars of source]

Such adjustment offsets the effect of local perturbation of $\gamma$ on the identifying moment condition. This component of DML is a property with respect to the population moment condition, i.e., irrelevant of multiway clustering or cross-fitting, and it is well-established in the GMM setting due to aforementioned literature, among others. Therefore, we take this condition as given for our analyses.

Let $\widehat\eta$ be some machine learners that are appropriate for multiway clustering data\footnote{For example, cluster-LASSO from belloni2016inference can be used for one-way clustering data. For two-way clustering panels, LASSO in chen2025inference can be employed. For clustering more than two dimensions, the multiplier bootstrap for jointly exchangeable arrays in chiang2023inference can be used for choosing the valid penalty levels.}, and let $\widehat\psi_N$ denote the empirical average with plug-in estimate $\widehat\eta$:

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

With some positive semi-definite weighting matrix $\widehat\Upsilon$ (e.g., the inverse of a multiway cluster-robust variance-covariance estimator of $\psi(X;\widehat\theta^{(0)},\widehat\eta)$ with some initial estimate $\widehat\theta^{(0)}$), the debiased GMM estimator of $\theta_0$ is defined as

align[align omitted — 150 chars of source]

When $\theta_0$ is exactly identified, the debiased GMM estimator reduces to $ \widehat\theta$ as a solution to

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

Our analyses are based on the general case ((ref)).

We denote the population and empirical Jacobian as \[ J_0 := -\partial_\theta \mathbb{E}[\psi(X;\theta,\eta_0)]\big|_{\theta=\theta_0}, \qquad \widehat J_N(\widehat\theta) := -\partial_\theta \widehat \psi_N(\theta)\big|_{\theta=\widehat\theta}. \] Suppose we have an interior minimizer $\widehat\theta$ from ((ref)), then the first-order condition holds as follows:

align[align omitted — 109 chars of source]

Let $f(\eta) = \psi(X,\theta_0,{\eta}) - \psi(X,\theta_0,{\eta}_0)$. By a standard mean-value expansion of $\psi_N(\widehat\theta)$ in ((ref)), we can write $ \sqrt{n}(\widehat \theta - \theta_0)$ as a function of a well-behaved term $ \mathbb{E}_N [\psi(X,\theta_0,{\eta}_0)]$, an empirical process term $\mathbb{G}_{n}\left(f(\widehat\eta)\right)$, an extra error term $ \mathbb{E}[ f(\eta)]_{\eta = \widehat\eta}$ due to the nuisance parameter estimation, as well as the empirical Jacobian $J_N(\widehat\theta)$ and feasible weighting matrix $\widehat\Upsilon$. As in the DML literature, $ \mathbb{E}[ f(\eta)]_{\eta = \widehat\eta}$ can be bounded by the orthogonality condition. It is relatively straightforward to deal with the empirical Jacobian term, and it is standard in the literature to put aside the weighting matrix estimation, as long as $\widehat{\Upsilon}\overset{p}{\to}\Upsilon$ to some positive-definite matrix $\Upsilon$. Now the difficult term left is $\mathbb{G}_{n}(f(\widehat\eta))$, and we bound it using a localisation approach through the maximum inequality under multiway clustering.

The idea of the localisation approach is that under the exchangeability and dissociation conditions, we can utilize the Hoeffding decomposition of the empirical process $\mathbb{G}_{n}\left(f(\widehat\eta)\right)$ and bound each of the decomposed terms by the maximum inequality in a neighborhood of $\eta_0$. As long as the first-step machine learner of the nuisance parameter lies in the neighborhood $\Gamma_n(\eta_0)$ with high probability, we can show $\mathbb{G}_{n}\left(f(\widehat\eta)\right)$ vanishes asymptotically. To formally define the neighborhood of $\eta_0\in\Gamma$, we equip $\Gamma$ with the $L_2(P)$ norm $\Vert.\Vert_{P,2}$.

assumptionLet $\Gamma_n(\eta_0)= \{\eta: \Vert \eta-\eta_0\Vert_{P,2} \le C_1 n^{-1/4}\} $ for a constant $C_1<\infty$ and $\mathcal{N}(\theta_0)$ be a shrinking neighborhood of $\theta_0$. \begin{enumerate}[(i)] • $(X_{\bm{i}})_{{\bm{i}}\in [\bm{N}]}$ satisfy Conditions (SE) and (D). • { $\psi(X,\theta_0,\eta)$ is continuous in $\eta$ a.s.}, $\mathbb{E}[\psi(X,\theta_0,\eta)]$ is Lipschitz-continuously Gateaux-differentiable on $\Gamma$ , and $\mathbb{E} \|\psi(X,\theta_0,{\eta}) - \psi(X,\theta_0,{\eta}_0)\|^2 \leq C_2\Vert \eta-\eta_0\Vert_{P,2}^2$ for a constant $C_2<\infty$. • $\psi(X,\theta,\eta)$ is differentiable in $\theta$, $\partial_\theta\psi(X,\theta_0,\eta)$ is continuous at $\eta_0$ a.s., and \\ $ \mathbb{E} [\sup_{\eta\in\Gamma_n(\eta_0)}\|\partial_\theta\psi(X,\theta_0,\eta)\|]<\infty$; There exists positive $B(X,\eta)$ such that \\ $\Vert\partial_\theta\psi(X,\theta,\eta) - \partial_\theta\psi(X,\theta_0,\eta) \Vert \leq B(X,\eta) \Vert \theta-\theta_0\Vert^\alpha$, $\forall \theta\in\mathcal{N}(\theta_0)$ with some $\alpha>0$; $B(X,\eta)> 0$ is continuous at $\eta_0$ a.s. and $ \mathbb{E}[\sup_{\eta\in\Gamma_n(\eta_0)}B(X,\eta)]<\infty$; . • $ \|\widehat\eta - \eta_0\|_{P,2}= o_P\left(n^{-1/4}\right)$ and $\partial_\eta \mathbb{E}[ \psi(X,\theta_0,\eta_0)](\widetilde \eta) = o\left(n^{-1/2}\right)$ for all $\widetilde\eta \in \Gamma$. • $\mathbb{E}\Vert \psi(X,\theta_0,\eta_0)\Vert^2<\infty$. • ${\rm rank}(J_0) = d$. \end{enumerate}

Assumption (ref)(i) characterizes the multiway clustered data by the exchangeability and dissociation conditions, which are standard in clustering robust inference literature. Assumptions (ref)(ii) and (iii) are score regularity conditions. Assumption (ref)(ii) is satisfied when the scores come from a likelihood function or moment conditions that are smooth in terms of the nuisance parameters. Assumption (ref)(iii) is a nonlinear counterpart of the linear-in-$\theta$ condition common in the DML literature. For a score that is twice differentiable in $\theta$, the existence of the integrable envelope $B(X,\eta)$ reduces to the integrability of the Hessian matrix locally. See Remark (ref) below for more details on the choice of $B(X,\eta)$. The locality in $\Gamma_n(\eta_0)$ ensures that $\eta$ takes values that do not explode up the envelope, e.g., $\eta$ as inverse probability weights. Assumption (ref)(iv) is a high-level condition governing the quality of nuisance parameter estimation. This requirement can be verified using existing theoretical results for a range of machine-learning estimators under multiway clustering; for instance, in the case of LASSO, it follows from Proposition 2 of chiang2023inference. Assumptions (ref)(v) and (vi) are standard and mild finite moment and full rank conditions.

remark[On $B$ in Assumption (ref)(iii)] $B(X,\eta)$ in Assumption (ref)(iii) is a local Hölder modulus term that controls how nonlinear the score is in $\theta$, uniform in $\eta\in\Gamma$. This is generally defined by, for a shrinking neighborhood of $\theta_0$, $\mathcal{N}(\theta_0)$, \begin{align*} B(X,\eta) :=\sup_{\theta\in\mathcal{N}(\theta_0)} \frac{\Vert \partial_\theta\psi(X,\theta,\eta) - \partial_\theta\psi(X,\theta_0,\eta)\Vert}{\Vert \theta-\theta_0\Vert^\alpha}, \quad \alpha>0. \end{align*} For twice differentiable $\psi(X,\theta,\eta)$ in $\theta$, by the mean-value theorem, there exists $\widetilde\theta $ such that $$ \partial_\theta\psi(X,\theta,\eta) - \partial_\theta\psi(X,\theta_0,\eta) = \partial_{\theta\theta}\psi(X,\widetilde\theta,\eta)(\theta-\theta_0) \quad s.t. \ \Vert \widetilde\theta - \theta_0 \Vert \leq \Vert \theta - \theta_0 \Vert.$$ Then, by setting $\alpha= 1$, we can take $ B(X,\eta) = \sup_{\theta\in \mathcal{N}(\theta_0)}\Vert \partial_{\theta\theta}\psi(X,\theta,\eta)\Vert$. If $\psi$ is linear in $\theta$, then this reduces to $ B(X,\eta) =\Vert \partial_{\theta\theta}\psi(X,\theta,\eta)\Vert$ which does not depend on $\theta$. To verify the integrability of $\sup_{\eta\in\Gamma_n(\eta_0)}B(X,\eta)$, one can often look for integrable dominating functions that do not depend on $\eta$. If $B(X,\eta)$ is Hölder modulus at $\eta_0$, then the desired integrability reduces to the integrability of $B(X,\eta_0)$ as well as the integrability of the Hölder modulus term at $\eta_0$. $\lozenge$

Following Chapter 3.6 in gine2016mathematical, a function class $\mathcal{F}$ on $\mathcal{S}$ with a measurable envelope $F$ is called Vapnik–Chervonenkis-type (VC-type) with characteristics $(A,v)$ if

align[align omitted — 171 chars of source]

where the supremum is taken over all finite discrete distributions. It can be shown that a wide range of commonly used models and estimators in econometrics, machine learning, and statistics give rise to sequences of function classes that satisfy this VC-type condition; see Section (ref) below for illustrative examples.

The following theorem presents our first main result, establishing the asymptotic linearity and asymptotic normality of a generic debiased GMM estimator without cross-fitting.

theorem[Asymptotic linearity and normality] Let $\widehat{\theta}$ be a solution to $(\ref{foc})$. Suppose Assumption (ref) holds, and $\widehat{\theta}\overset{p}{\to} \theta_0$, $\widehat\Upsilon\overset{p}{\to} \Upsilon$ for some positive-definite limit $\Upsilon$ as $n\to\infty$. For each $n\in \mathbb N$, let $F_n$ be an envelope for the function class $ \mathcal{F}_n: = \{ f(\eta) - \mathbb{E}[ f(\eta)]:\eta\in \Gamma_n(\eta_0)\}$. Then suppose that, for some $q\geq 2$ and for each $n$, $\|F_n\|_{P,q}<\infty$\footnote{This condition need not hold uniformly in $n$; in particular, $\|F_n\|_{P,q}\to\infty$ is allowed.} and \begin{align} $\mathcal{F}_{n}$ is a VC-type class with characteristics $A_n\ge (e^{2(K-1)/16}\vee e)$ and $v_n\ge1$, \end{align} then we have (i) the following linear representation holds \begin{align*} \sqrt{n}(\widehat{\theta}-\theta_0) = &\left( J_0' \Upsilon J_0\right)^{-1} J_0' \Upsilon \sqrt{n} \mathbb{E}_N [\psi(X,\theta_0,{\eta}_0)] + \sum_{k=1}^K\sum_{\bm{e}\in \mathcal{E}_k} O_P(\rho_{n,k}) +o_P(1) \\ where \ \rho_{n,k}=&\left(\frac{v_n\log(A_n\vee \overline{N} )}{n^{1-1/2k}}\right)^{k/2} \vee \left(\frac{\left\| F_n \right\|_{P,q} \{v_n\log(A_n\vee \overline{N} )\}}{n^{1/2-1/q}}\right)^{k}. \end{align*} (ii) If, additionally, (a) it holds that \begin{align} \frac{v_n\log(A_n\vee \overline{N} )}{n^{1-1/2k}}\vee \frac{\left\| F_n \right\|_{P,q} \{v_n\log(A_n\vee \overline{N} )\}}{n^{1/2-1/q}}=o(1), \end{align} and (b) the smallest eigenvalue of $\Psi_0$, defined below, is bounded from below by some constant $c>0$, \begin{align*} \sqrt{n}(\widehat{\theta}-\theta_0)\overset{d}{\to} N(0,V), \end{align*} where $V :=\left(J_0' \Upsilon J_0\right)^{-1} J_0' \Upsilon \Psi_0 \Upsilon J_0 \left(J_0' \Upsilon J_0\right)^{-1}$ and \begin{align*} \Psi_0 : =\mu_1 {\rm Var}( \mathbb{E}[\psi(X,\theta_0,\eta_0)|U_{1,0,...,0}]) +...+ \mu_K {\rm Var}( \mathbb{E}[\psi(X,\theta_0,\eta_0)|U_{0,0,...,1}]). \end{align*} where $\mu_k =\lim_{n\to\infty} \frac{n}{N_k}$ for $k=1,...,K$; $\{U_{\bm i}\}_{\bm i> \bm 0}$ are as defined in (ref).

A proof can be found in Section (ref) in the appendix. The additional condition (b) in statement (2) of the theorem is a non-degeneracy requirement, ensuring that at least one clustering dimension enters the score in a linear manner. This condition is mild for larger $K $'s, as it only requires that at least a single one latent shock of the $K$ clustering dimensions has a non-trivial effect on $\psi$. This condition was also imposed in, e.g., davezies2021empirical, chiang2022multiway and chiang2023inference. In the case of i.i.d data, the non-degeneracy condition does not hold. In such conventional settings, the asymptotic normality result for the full-sample DML approach has been established in belloni2015uniform, among others, while here we focus on the non-degenerate case.

To build intuition, note that the first-order condition implies the expansion

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

for some invertible matrix $M$, where $ f(\eta)=\psi(X,\theta_0,\eta)-\psi(X,\theta_0,\eta_0). $ Such a decomposition is standard in semiparametric theory; see, for example, andrews1994asymptotics.

The first term is asymptotically normal. The second term is controlled by the orthogonality condition (ref), and is therefore negligible. Consequently, the main technical challenge is to control the localised empirical process $\mathbb{G}_n\bigl(f(\widehat \eta)\bigr),$ which is non-standard due to the dependence of $\widehat \eta$ on the full sample. A standard approach is cross-fitting, which removes this dependence. Conditional on $\widehat \eta$, one may apply Hoeffding-type decomposition and Markov's inequality to obtain, under suitable smoothness conditions, with probability $1-o(1)$ \[ \mathbb{G}_n\bigl(f(\widehat \eta)\bigr) \;\lesssim\; \|\widehat \eta - \eta_0\|^{v}_{P,2} \quad \text{for some } v>0. \]

An alternative is a localisation argument. Suppose there exists a sequence of shrinking function classes $\{\mathcal{F}_n\}$ such that $ \mathbb{P}(\widehat \eta \in \mathcal{F}_n) \to 1. $ Then, with probability $1-o(1)$, \[ \bigl|\mathbb{G}_n(f(\widehat \eta))\bigr| \;\le\; \sup_{\eta\in\mathcal{F}_n} \bigl|\mathbb{G}_n(f(\eta))\bigr|, \] which removes the stochastic dependence on $\widehat \eta$. Such classes $\{\mathcal{F}_n\}$ can often be constructed tightly when the convergence rate of $\|\widehat\eta - \eta_0\|_{P,2}$ is available. This is typically the case, as the same rate is also needed to verify the orthogonality condition regardless of whether cross-fitting is used.

In the classical semiparametric literature, the function class $\mathcal{F}$ is typically fixed, and stochastic equicontinuity follows from standard uniform (functional) CLT arguments. In contrast, with machine-learning first stages, a fixed $\mathcal{F}$ is generally too large to control, necessitating shrinking (localised) classes. This localisation strategy underlies the i.i.d.\ analyses of belloni2015uniform further generalised in belloni2018uniformly, which rely on maximal inequalities from CCK2014AoS to control the supremum. Since comparable results are unavailable under multiway clustering, we develop the required global and local maximal inequalities in Section (ref).

remark[Choosing the envelopes $F_n$] The envelope $F_n$ of the function class $\mathcal{F}_n$ can be taken as $\sup_{\eta\in\Gamma_n(\eta_0)} \left|f(\eta) - \mathbb{E}[f(\eta)]\right| = \sup_{\eta\in\Gamma_n(\eta_0)} \left|\psi(X,\theta_0,\eta) - \mathbb{E}[\psi(X,\theta_0,\eta) ] - \psi(X,\theta_0,\eta_0) \right|$. The $L^q $ integrability of this object can be ensured by $\mathbb{E}\left[\sup_{\eta\in\Gamma_n(\eta_0)}|\psi(X,\theta_0,\eta)|^q\right]<\infty$, which is in turn verifiable in specific GMM models. For example, if $\psi(X,\theta_0,\eta)$ is Hölder continuous in $\eta$ w.p.1, then the $L^q$ integrability of $\psi(X,\theta_0,\eta_0)$ and the shrinking neighborhood $\Gamma_n(\eta_0)$ delivers the $L^q $ integrability of $F_n$. $\lozenge$

remark[Alternative asymptotics for semiparametric estimation] Making the empirical process component in the asymptotic expansion negligible is not the only route to asymptotic linearity and normality in DML and related two-step estimation problems. An alternative arises when the nuisance parameter $\eta$ itself admits an asymptotically linear representation and satisfies a central limit theorem. A canonical example is the density-weighted average derivative estimator studied by powell1989semiparametric. In such settings, valid inference can instead be conducted under small-bandwidth asymptotics, as developed by cattaneo2014small. Under this regime, the empirical process term $\mathbb{G}_n(f(\widehat\eta))$ need not vanish asymptotically; rather, its limiting distribution can be explicitly characterised using a CLT for quadratic forms\footnote{See, e.g. de1987central.} and may be of the same order as, or even dominate, the usual asymptotic linear component. This framework is particularly relevant when the nuisance parameter is estimated using kernel-based methods. It is also worth noting that (leave-one-out) cross-fitting continues to play a role in this literature. $\lozenge$

remark[Multiway-clustering stability] One may alternatively pursue asymptotic normality of the debiased GMM estimator without using a maximal inequality by directly controlling the empirical process component $\mathbb{G}_n(f(\widehat\eta))$ appearing in Remark (ref), through a multiway-clustering stability condition analogous to that of chen2022debiased. To define such a condition, for $k=1,\ldots,K$ and $i_k'=1,\ldots,N_k$, let $X^{(i_k=i_k')}=\{X_{\bm i}:\bm i\in[\bm N],\, i_k=i_k'\}$ and let $\widetilde X^{(i_k=i_k')}$ denote an independent copy, and for $\bm i'=(i_1',\ldots,i_K')\in[\bm N]$ define $X^{\neg(\bm i')}$ as the dataset obtained by replacing $\bigcup_{k=1}^K X^{(i_k=i_k')}$ with $\bigcup_{k=1}^K \widetilde X^{(i_k=i_k')}$; by the dissociation condition, $X_{\bm i'}$ is independent of $X^{\neg(\bm i')}$, and we denote by $\widehat\eta^{\neg(\bm i')}$ the corresponding nuisance estimator. An analogue of the stability condition in chen2022debiased then requires that, for all $j=1,\ldots,d$, \begin{align*} \max_{\bm{i} \in[\bm N]} \mathbb{E}\left\vert \psi_j(X_{\bm i},\theta_0,\widehat{\eta}) -\psi_j(X_{\bm i},\theta_0,\widehat{\eta}^{\neg (\bm{i})}) \right\vert &= o(n^{-1/2}) \\ \max_{\bm{i} \in[\bm N]} \left( \mathbb{E}\left\vert \psi_j(X_{\bm i},\theta_0,\widehat{\eta}) -\psi_j(X_{\bm i},\theta_0,\widehat{\eta}^{\neg (\bm{i})}) \right\vert^2\right)^{1/2} &= o(n^{-1/2}). \end{align*} which would render the empirical process term asymptotically negligible without localisation or sample splitting; however, since no existing first-stage machine learning estimators are known to satisfy such rate conditions under multiway clustered dependence, we leave this approach for future research. $\lozenge$

Verification of complexity and rate conditions: three examples

Under Assumption (ref), to apply Theorem (ref) it suffices to verify the high-level VC-type condition in (ref) and the rate condition in (ref). Verifying these conditions is not entirely straightforward in general, owing to their abstract nature. Below we discuss several examples of machine learning estimators for the first-stage nuisance function $\eta$ that satisfy these two conditions.

To isolate the role of the first-stage nuisance parameter learner, we consider $\psi(\cdot,\eta)$ as a map in $\eta$ that preserves the VC-type properties of the underlying function class $\mathcal{G}_n$ to which $\eta$ belongs. For instance, $\psi(\cdot,\eta)$ may be a monotone or Lipschitz transformation of $\eta$, or a finite combination such as sums, products, minima, or maxima; see Section 3.6 of gine2016mathematical.

theorem[Complexity/rate for machine learners] Suppose $\psi(.,\eta)$ is a map that preserves the VC-type characteristics of $\mathcal{G}_n $ and $\Vert F_n\Vert_{P,q}<\infty$ with some $q>4$. Then under each of the following three cases, Conditions (ref) and (ref) in Theorem (ref) hold. \begin{enumerate} • (Generalised linear models with $\ell^1$-regularisation) If for a monotonic link $g$ $$\mathcal{G}_n = \{x\mapsto g(x^T\beta): \beta\in \mathbb{R}^p,\Vert\beta \Vert_0<s, \Vert\beta \Vert_2<\infty \},$$ then $\mathcal{G}_n$ is VC-type with $v_n=s$ and $A_n = \frac{C \ ep}{s}\vee (e^{2(K-1)/16}\vee e)$ for some $C<\infty$. Furthermore, ${s\log(p/s)} \vee {s\log(\overline{N})} = o(n^{1/4})$. • (Regression trees) Let $\{R_l\}_{l=1}^L$ denote random partitions of a regression tree with axis-aligned threshold splits based on features in $\mathbb{R}^p$, and the output on each leaf is a constant. The corresponding function class can be written as $$\mathcal{G}_n = \left\{x\mapsto g(x): g(x) = \sum_{l=1}^L \mu_{l} 1\{x\in R_l\},\ \bigcup_{l=1}^LR_l = \mathbb{R}^p, \ |\mu_l|<\infty\right\}.$$ Then, $\mathcal{G}_n$ is a VC-subgraph class with pseudo VC-dimension of order $O(L\log(Lp))$. For some $C<\infty$, $v_n=2CL\log(2Lp)$ and $A_n= C\vee(e^{2(K-1)/16}\vee e)$. Further assume that $ {L\log(2Lp)\log(A_n \vee \overline{N})} = o(n^{1/4})$. • (Deep neural networks) Consider a feed-forward neural network with the ReLU activation function, $p$ features, $L-2$ hidden layers, $U$ total hidden units, and $W-1$ total parameters. Let $\mathcal{G}_n$ be the class of functions generated by such a neural network with a fixed structure. Then, $\mathcal{G}_n$ is a VC-subgraph class with pseudo VC-dimension of order $O(LW\log(pU))$. For some $C<\infty$, $v_n=2CLW\log(pU)$ and $A_n= C\vee(e^{2(K-1)/16}\vee e)$. Further assume that ${LW\log(pU)\log(A_n \vee \overline{N})} = o(n^{1/4})$. \end{enumerate}

A proof can be found in Section (ref) in the appendix.

Case (i) admits generalized linear sparse models, including sparse linear, logit, and exponential models as special cases. These models correspond to LASSO-type ($\ell^1$-penalty) machine learners for a sparse generalised linear model. Cases (ii) and (iii) correspond to the regression tree and deep neural networks, respectively. More details are given in the appendix on how the rate conditions are obtained. Basic definitions and textbook treatments of these methods can be found in e.g. chernozhukov2024applied.

Variance estimation.

For hypothesis testing using the results given in Theorem (ref), a missing piece is the unknown asymptotic variance $V$. In this section, we propose a full-sample variance estimator that takes into account (1) multiway clustering dependence and (2) estimation errors from both high-dimensional nuisance estimation and the GMM estimation. To account for the multiway clustering dependence, we follow the formulation of the multiway clustering-robust variance estimator in DDG2018, except that the empirical scores here are replaced by the Neyman orthogonalised scores with estimated nuisance parameters.

For any $\bm i, \bm j \in [\bm N]$, let $i_k$ and $j_k$ denote their $k$-th elements, and let $\mathbbm 1_{k}\{\bm i, \bm j\}$ indicate whether the two observations $\bm i, \bm j$ share the same cluster at $k$-th dimension, i.e., $\mathbbm 1_{k}\{\bm i, \bm j\} = \mathbbm 1\{ i_k = j_k \}.$ We define the estimator for the middle term $\Psi_0$ as follows:

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

Then the estimator for $V$ is given as follows:

align[align omitted — 377 chars of source]

In practice, if there are more than one observation in some cells $\bm i$, we simply aggregate within each cell by replacing $\psi(X_{\bm i},\widehat{\theta},\widehat{\eta})$ with the sum of empirical scores $\psi$ within that cell. Since this generalisation would not change the main analysis except for complications in notations, we focus on the case with exactly one observation in each cell.

theorem[Consistent variance estimation] Under the same conditions as in Theorem (ref), as well as \begin{align} & \mathbb{E}\left[ \sup_{\eta\in \Gamma_n(\eta_0)} \left\Vert \psi(X,\theta_0,\eta)\right\Vert^2\right]< \infty, , \\ &\mathbb{E}\left[\sup_{\eta\in \Gamma_n(\eta_0)} \left\Vert \partial_\theta\psi(X,\theta_0,\eta) \right\Vert^2\right]< \infty , \\ &\mathbb{E}\left[\sup_{\eta\in\Gamma_n(\eta_0)}B^2(X,\eta)\right]<\infty , \end{align} { where $B(X,\eta)$ is defined in Assumption (ref)(iii)}, then $\widehat{V}\overset{p}{\to}V$ as $n\to\infty$.

A proof can be found in Section (ref) in the appendix.

Theorem (ref) establishes the consistency of the variance estimator using the full sample, under the same non-degeneracy condition as in the second statement of Theorem (ref). The extra moment conditions mildly strengthen the moment conditions in Theorem (ref). As in Theorem (ref), these local integrability conditions can be delivered by integrability conditions of the score, Jacobian, and the Hessian, given enough smoothness in $\eta$.

$\widehat{V}$ is positive semi-definite by construction because each $\widehat\Psi_{N,k}(\widehat\theta)$ is positive semi-definite mechanically. This can be seen easily in the case $K=2$, in which case each $\widehat\Psi_{N,k}(\widehat\theta)$ reduces to a one-way cluster variance estimator. A caveat is that when none of the cluster matters, e.g., i.i.d. data, the non-degeneracy condition can fail, and this variance estimator would overestimate the asymptotic variance $V$ and result in a conservative test, which is well-known in the cluster robust inference literature (e.g., see mackinnon2021wild). A potential fix for this issue is to remove double-counting terms in $\widehat\Psi_{N}(\widehat\theta)$ by defining

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

where $I_{\bm e}\{\bm i, \bm j\}= \mathbbm 1\{\bm{e}\odot \bm{i} = \bm{e}\odot \bm{j}\} $, which instead indicates whether the two observations share the same clusters over the whole support of $\bm e$. This is basically a DML version of the $K$-way generalisation of variance estimator proposed in Cameron2011, referred to as the CGM estimator\footnote{Other candidates for inference procedures include the modified multiway empirical likelihood and the modified multiway jackknife variance estimator proposed in chiang2024multiway. While these methods are computationally more demanding, they can potentially deliver improved higher-order asymptotic properties; see Theorem 3 therein.}. In a parametric setting, DDG2018 shows that these two types of variance estimators are both consistent for the asymptotic variance under non-degeneracy. However, the CGM estimator involves more terms to calculate and is not guaranteed to be positive semi-definite. In practice, it is rarely true that multi-dimensional data is i.i.d because of the common existence of unobserved heterogeneous effects.

For some degenerate yet cluster-dependent scenarios, such as those studied by menzel2021bootstrap, the estimator itself may fail to satisfy asymptotic normality. In such cases, standard inference procedures—including CGM-type variance estimators as well as the approach proposed in this paper—become invalid. In this context, menzel2021bootstrap proposes bootstrap-based inference methods that are uniformly valid, but they require the choice of tuning parameters and is typically conservative. In the two-way clustering setting, davezies2025analytic develop a simple analytical inference that remains valid under non-Gaussian degeneracy, while avoiding the need for tuning parameters. Complementarily, hounyo2025projection propose bootstrap procedures that are adaptive across a range of non-degenerate and (Gaussian) degenerate cases.

Extending these approaches to our setting is substantially more involved due to the presence of two-step estimators and machine learning-based first stages. In particular, degeneracy changes the effective stochastic order of the leading term, thereby tightening the rate requirements on the first-stage estimators to ensure that their estimation error is asymptotically negligible. Moreover, incorporating the strategy of davezies2025analytic is highly non-trivial in our framework, as it relies on conditioning arguments that, in the presence of multiway dependence and generated regressors, further complicate the first-stage convergence requirements. Addressing these challenges would require new techniques, and we leave them for future research.”

\vskip 0.15in

Maximum Inequalities under separate exchangeability

In this section, we establish inequalities that control the $q$-th moment of the supremum of the empirical process, $ \mathbb{E}\bigl[\|\mathbb{G}_{n}\|_\mathcal{F}^q \bigr] , $ for some $q \in [1,\infty)$ for SE arrays. Throughout this section, assume without loss of generality that $\mathbb{E}[f(X_{\bm{1}})] = 0$ for all $f \in \mathcal{F}$. Before presenting the main results, let us first introduce the Hoeffding-type decomposition from chiang2023inference. For any $\bm{i}\in [\bm{N}]$, define

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

We then define recursively for $k=1,2,\dots, K$ that

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

and for $\bm{e}\in \bigcup_{k=2}^K\mathcal{E}_k$ set

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

Note that by the AHK representation (ref), for a fixed $\bm{e}$ the distributions of \[ (P_{\bm{e}}f)\Bigl(\{U_{\bm{i}\odot \bm{e}'}\}_{\bm{e}'\le \bm{e}}\Bigr) \quad\text{and}\quad (\pi_{\bm{e}}f)\Bigl(\{U_{\bm{i}\odot \bm{e}'}\}_{\bm{e}'\le \bm{e}}\Bigr) \] do not depend on the index $\bm{i}$. Hence, we shall write $P_{\bm{e}}f$ and $\pi_{\bm{e}}f$ for a generic $\bm{i}$.

Now, fix any $1\le k\le K$ and let $\bm{e}\in \mathcal{E}_k$. Then, by Lemma 1 in chiang2023inference, for any $\ell\in \mathrm{supp}(\bm{e})$ the random variable $ (\pi_{\bm{e}}f)\Bigl(\{U_{\bm{i}\odot \bm{e}'}\}_{\bm{e}'\le \bm{e}}\Bigr) $ has mean zero conditionally on $\{U_{\bm{i}\odot \bm{e}'}\}_{\bm{e}'\le \bm{e}-\bm{e}_{\ell}}$. In addition, define $ I_{\bm{N},\bm{e}} = \{\bm{i}\odot \bm{e} : \bm{i}\in [\bm{N}]\}. $ Then, we have $ \bigl|I_{\bm{N},\bm{e}}\bigr| = \prod_{k'\in\mathrm{supp}(\bm{e})} N_{k'}. $ Accordingly, define

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

We now obtain the Hoeffding-type decomposition

align[align omitted — 117 chars of source]

To bound $\mathbb{E}\bigl[ \|\mathbb{G}_{n}(f)\|_{\mathcal{F}} \bigr] $, it thus suffices to control each individual term $\mathbb{E}\bigl[\|H_{\bm{N}}^{\bm{e}}(f)\|_{\mathcal{F}}\bigr]$ separately.

Finally, fix any $1\le k\le K$ and $\bm{e}\in \mathcal{E}_k$. Define the uniform entropy integral by

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

where $ P_{\bm{e}}\mathcal{F} := \{P_{\bm{e}}f : f\in \mathcal{F}\}, $ and the supremum is taken over all finite discrete distributions $Q$.

The following result is a general global maximal inequality for SE empirical processes with an arbitrary index order $K$ and for a general order of moment $q\in[1,\infty)$. Its proof follows the arguments in the proof of Corollary B.1 in chiang2023inference with some modifications to account for a more general class of functions.

theorem[Global maximal inequality for SE processes] Suppose $\mathcal{F}:\mathcal{S}\to \mathbb{R}$ is a pointwise measurable class of functions. Let $(X_{\bm{i}})_{{\bm{i}}\in [\bm{N}]}$ be a sample from $S$-valued separately exchangeable random vectors $(X_{\bm{i}})_{{\bm{i}}\in \mathbb{N}^K}$. Pick any $1 \le k \le K$ and $\bm{e} \in \mathcal{E}_{k}$. Then, for any $q \in [1,\infty)$, we have \[ |I_{{\bm{N}},{\bm{e}}}|^{1/2}\left (\mathbb{E} \left [ \left \| H_{{\bm{N}}}^{\bm{e}}(f) \right \|_{\mathcal{F}}^{q} \right ] \right)^{1/q} \lesssim J_{\bm{e}}(1) \| F \|_{P,q\vee 2}. \]

A proof can be found in Section (ref) in the appendix.

Although the global maximal inequality works for general $q$, in the case that the supremum of the first absolute moment is concerned, local maximal inequalities usually provides shaper bounds. The following is a novel local maximal inequality for SE empirical processes.

theorem[Local maximal inequality for SE processes] Suppose $\mathcal{F}:\mathcal{S}\to \mathbb{R}$ is a pointwise measurable class of functions. Let $(X_{\bm{i}})_{{\bm{i}}\in [\bm{N}]}$ be a sample from $S$-valued separately exchangeable random vectors $(X_{\bm{i}})_{{\bm{i}}\in \mathbb{N}^K}$. Set ${\bm{e}}\in \{0,1\}^K$ and let $\sigma_{\bm{e}}$ be a constant such that $\sup_{f\in\mathcal{F}}\|P_{\bm{e}} f \|_{P,2}\le \sigma_{\bm{e}} \le \|P_{\bm{e}} F\|_{P,2}$, \[ \delta_{\bm{e}}=\sigma_{\bm{e}}/\|P_{\bm{e}} F\|_{P,2}\quad \text{ and } \quad M_{\bm{e}}=\max_{t\in [n]}(P_{\bm{e}} F)\left(\{U_{(t,...,t)\odot {\bm{e}}'}\}_{{\bm{e}}'\le {\bm{e}}}\right),\] then \begin{align*} |I_{{\bm{N}},{\bm{e}}}|^{1/2}\mathbb{E} \left [ \left \| H_{{\bm{N}}}^{\bm{e}}(f) \right \|_{\mathcal{F}} \right ] \lesssim J_{\bm{e}}(\delta_{\bm{e}}) \|P_{\bm{e}} F\|_{P,2} +\frac{J_{\bm{e}}^2(\delta_{\bm{e}})\|M_{\bm{e}}\|_{P,2}}{\sqrt{n}\delta_{\bm{e}}^2}. \end{align*}

A proof can be found in Section (ref) in the appendix.

remarkAlthough our proof strategy broadly follows that of Theorem 5.1 in ChenKato2019b for U-processes with modifications accounting for SE structures, a key divergence arises. In Theorem 5.1, the Hoffmann–J\o rgensen inequality (which requires independence) is applied via the classical $U$-statistic technique of Hoeffding averaging (see, for example, Section 5.1.6 in serfling1980approximation). However, the more intricate dependence structure inherent to separately exchangeable arrays renders Hoeffding averaging inapplicable in our context. To address this challenge, we introduce an alternative approach by establishing Lemma (ref), which partitions the index set $I_{{\bm{N}},{\bm{e}}}$ into transversal groups of size $n$ (see Lemma (ref) for its definition). Together with AHK representation (ref), this yields i.i.d. elements within each group, thereby facilitating the application of the Hoffmann–J\o rgensen inequality. $\lozenge$

In practice, bounding the uniform entropy integrals appearing on the right-hand side of maximal inequalities can be involved. Fortunately, many function classes arising in econometrics, machine learning, and statistics can be shown to be of VC-type, in the sense of (ref). Under this assumption, the entropy terms entering the maximal inequalities admit substantially simpler bounds. We therefore derive a local maximal inequality under the VC-type condition, which is the version utilised in the proof of Theorem (ref).

corollaryUnder the same setting as in Theorem (ref). In addition, suppose $\mathcal{F}$ is of VC-type with characteristics $A\ge (e^{2(K-1)}/16)\vee e$ and $v\ge 1$, then for each ${\bm{e}}\in \mathcal{E}_k$, one has \begin{align*} |I_{{\bm{N}},{\bm{e}}}|^{1/2}\mathbb{E} \left [ \left \| H_{{\bm{N}}}^{\bm{e}}(f) \right \|_{\mathcal{F}} \right ] \lesssim& \sigma_{\bm{e}}\{v\log (A\|P_{\bm{e}} F\|_{P,2}/\sigma_{\bm{e}})\}^{k/2} +\frac{\|M_{\bm{e}}\|_{P,2}}{\sqrt{n}}\{v\log (A \|P_{\bm{e}} F\|_{P,2}/\sigma_{\bm{e}})\}^k \\ \lesssim& \sigma_{\bm{e}}\{v\log (A\vee \overline N)\}^{k/2} +\frac{\|M_{\bm{e}}\|_{P,2}}{\sqrt{n}}\{v\log (A\vee \overline N)\}^k. \end{align*}

A proof is provided in Section (ref) in the appendix.

Conclusion

This paper develops a cross-fitting-free asymptotic theory for two-step debiased GMM estimators under multiway clustered dependence and shows that valid inference in such settings hinges on new empirical process techniques. By combining orthogonal moment conditions with a localisation-based argument, we demonstrate that the impact of high-dimensional or nonparametric nuisance estimation can be controlled without sample splitting, even when the effective sample size is determined by the number of independent cluster units. The resulting estimators are shown to be asymptotically linear and normal under separately exchangeable sampling, providing a practical inference framework for empirically relevant clustered environments. A central contribution is the derivation of new global and local maximal inequalities for possibly uncountable, pointwise measurable function classes under multiway dependence, which fill a gap in the existing theory and may be useful beyond the DML context. Future work may further explore stability-type conditions under multiway clustering and extend these tools to other two-step and high-dimensional problems with complex dependence structures.