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.
77,132 characters
Limit Theory for U-Statistics under Clustered and Weakly Dependent Data
\hypersetup{pageanchor=false}
\pagenumbering{gobble}
\title{Limit Theory for U-Statistics under Clustered and Weakly Dependent Data}
\date{\today }
\author{Emmanuel Selorm Tsyawo\footnote{email: [email removed], Department of Economics, Finance and Legal Studies, Culverhouse College of Business, University of Alabama}}
\maketitle
\thispagestyle{empty}
\begin{abstract}
\noindent This paper develops asymptotic theory and feasible inference for unbounded-kernel order-$k$ $U$-statistics under clustered sampling and weakly dependent time-series sampling. The analysis first builds the complete order-$2$ pipeline, moving from clustered data to exact $m$-dependence and then to near-epoch dependence. The same logic is subsequently extended to general order $k\geq2$. Under clustered sampling, the theory allows arbitrary within-cluster dependence and growing, unbalanced cluster sizes. Under weak dependence, an i.i.d.-based approximating sequence carries the exact-$m$ theory to near-epoch-dependent processes. The common combinatorial device is a vertical rearrangement, which isolates sampling-generic tuples, where the first-order Hoeffding projection is analysed, from collision terms and higher-order degenerate remainders, which are controlled explicitly. Cluster-robust and HAC estimators of the covariance of the first-order projection, needed for feasible inference, are shown to be consistent.
\vspace{.5cm}
\noindent \textit{Keywords:} $U$-statistics, clustered sampling, weak dependence, near-epoch dependence, central limit theorem, cluster-robust inference, HAC covariance estimation
\vspace{.5cm}
\noindent \textit{JEL classification: C12, C14, C21, C22}
\end{abstract}
\begin{refsection}
\newpage
\hypersetup{pageanchor=true}
\pagenumbering{arabic}
\section{Introduction}\label{Sect:Introduction}
$U$-statistics arise naturally in economics and finance whenever the object of interest is built from pairwise or higher-order comparisons rather than from one-observation moments. Direct economic examples include inequality, distribution-shape, and rank-dependence measures such as the Gini mean difference, L-moments, Kendall's tau, and Spearman's rho, e.g., \citet{bhattacharya2007inference,darku2020gini,hosking1990lmoments,dehling2017testing}. $U$-statistics also appear more broadly in score and criterion functions associated with pairwise differencing estimators, e.g., \citet{honore1994pairwise,jochmans2013pairwise}, integrated conditional moment estimators and \(\chi^2\) specification tests, e.g., \citet{dominguez2004consistent,escanciano2018simple,jiang-tsyawo-2026consistent}, and gravity models, e.g. \citet{jochmans-2017-two_way_gravity,yang-zhang-2023-three_way_gravity}. Under dependence, these objects are often easy to motivate substantively but harder to analyse, because their summands overlap and the usual independent-tuple geometry breaks down. The difficulty is therefore not only point estimation, but also feasible covariance estimation for the first-order projection, whether the $U$-statistic is itself the target of inference or enters via an estimation or specification testing framework.
The results developed here treat these objects under two dependence structures usually handled separately: clustered sampling and weakly dependent time series. The unifying combinatorial device is a vertical rearrangement. In the clustered design, columns are formed so that no column contains two observations from the same cluster; in the fixed-$m$ time-series design, the analogue is the set of vertical partitions modulo $(m+1)$, whose observations are more than $m$ periods apart. The leading term of an order-$k$ $U$-statistic then comes from sampling-generic tuples: cluster-generic tuples in the clustered design and lag-generic tuples in the fixed-$m$ design, where ordinary Hoeffding projection logic applies. The remaining tuples are same-cluster collision terms in the clustered design and lag-collision terms in the temporal design. This makes the proof modular, i.e., the vertical rearrangement standardises the treatment of collision terms and higher-order degenerate terms, while the remaining first-order projection is handled by the probability theory appropriate to the sampling scheme.
\citet{hansen2019asymptotic} gives the clustered-sum central limit theorem that underlies the linear projection step, allowing arbitrary within-cluster dependence and heterogeneous, unbounded cluster growth. The clustered results below extend that input to $U$-statistics by isolating same-cluster collision terms, extracting the Hoeffding linear term, and controlling the higher-order degenerate remainders. This complements bounded-cluster two-sample $U$-statistic results such as \citet{Lee-Dehling-2005-generalized}, where the residual terms are lower order under uniformly bounded cluster sizes. It also sits alongside recent clustered rank-statistic work, such as the two-sample Mood statistic of \citet{suzuki2025mood}, which allows cluster-size heterogeneity in a two-sample setting.
For weakly dependent time series, the main nearby results cover related but distinct settings. \citet{sen1963properties} studies $U$-statistics for stationary $m$-dependent processes by separating non-serial tuples from tuples containing short-lag interactions and reducing the limit theory to the first-order projection. \citet{fischer2016multivariate} proves general-order $U$-statistic central limit theorems for strongly mixing sequences with bounded kernels satisfying an extended variation condition, and leaves the near-epoch-dependent multivariate extension as a conjectural direction. \citet{fischer2017robust} develops a general NED theory for multivariate kernels and Generalised Linear (GL) statistics, including an unbounded-kernel CLT and a HAC variance estimator in the bounded-kernel case. \citet{dehling2017testing}, working with P-near-epoch dependence on an absolutely regular process, treats the order-2 Kendall setting and uses empirical Hoeffding projections for HAC inference. The weak-dependence results below instead work within the present $L_2$-NED framework, with kernel smoothness or tail control used to transfer finite-window approximation to the $U$-statistic, while placing the time-series theory inside the same order-2/order-$k$ progression as the clustered results.
The paper makes three contributions. First, it identifies a common combinatorial structure behind clustered sampling and finite-range temporal dependence, where, after a suitable rearrangement, sampling-generic tuples are built from pairwise-independent units, while same-cluster and lag collisions are counted and controlled. Second, it uses that structure to separate the common remainder argument from the sampling-specific linear-projection step, yielding order-$k$ $U$-statistic laws of large numbers and central limit theorems under arbitrary within-cluster dependence with growing, heterogeneous cluster sizes, and under weakly dependent time series carried from exact $m$-dependence to near-epoch dependence through an i.i.d.-based approximating sequence. Third, it pairs that limit theory with feasible inference by providing a cluster-robust plug-in estimator in the clustered case and a HAC estimator based on the empirical first-order projection in the time-series case.
The remainder of the paper is organised as follows. \Cref{Sect:Order2Pipeline} develops the order-$2$ results, while \Cref{Sect:OrderKTheory} extends them to general order $k$. \Cref{Sect:FeasibleInference} develops the corresponding covariance estimators. \Cref{Sect:Empirical} presents the empirical applications to which simulations in \Cref{Sect:Simulations} are calibrated. Finally, \Cref{Sect:Conclusion} concludes. The proofs of all theoretical results are collected in the appendix.
\paragraph{Notation:} $\|X\|_p:=\big(\mathbb{E}[\|X\|^p]\big)^{1/p}$ denotes the $L^p$-norm of a random vector $X$ for $p\geq1$, where the inner $\|\cdot\|$ is the Euclidean norm; $[N]:=\{1,\ldots,N\}$ for $N\in\mathbb N$; and $\mathbb{V}[\cdot]$ denotes variance. For an ordered tuple \(\mathbf i=(i_1,\ldots,i_k)\in[N]^k\), \(\mathbf i_{\neq}\) means \(i_a\neq i_b\) for all \(a\neq b\), and such tuples are \emph{index-generic}. The set $\mathcal{G}_g$ collects units in cluster $g\in[G]$, $n_g$ is its cardinality, and $g(i):=\sum_{g=1}^G g\mathbbm{1}\{\mathcal{G}_g\ni i\}$ denotes the cluster of unit $i\in[n]$; the map \(g\) is applied componentwise to tuples, so \(g(\mathbf i):=(g(i_1),\ldots,g(i_k))\). Under clustered sampling, \(g(\mathbf i)_{\neq}\) means that the tuple has pairwise distinct cluster labels, and such tuples are \emph{cluster-generic}; under exact \(m\)-dependence, an index-generic tuple is \emph{lag-generic} when \(|i_a-i_b|>m\) for all \(a\neq b\). When the sampling scheme is clear, \emph{sampling-generic} means cluster-generic in the clustered design and lag-generic in the exact-\(m\) or \(m\)-dependent approximating design; the corresponding collision terms are index-generic tuples that fail the relevant sampling-generic condition. All random variables are defined on an underlying probability space \((\Omega,\mathcal A,\mathbb{P})\); the state space is standard Borel, kernels are measurable, and conditional expectations are taken as fixed measurable versions. The constants $r\geq2$, $C>0$, and $\lambda>0$ are used throughout, and the value of $C$ may change across occurrences.
\section{Order-\texorpdfstring{$2$}{2} Theory}\label{Sect:Order2Pipeline}
This section develops the full argument for order-$2$ $U$-statistics. It begins with clustered sampling, where the vertical rearrangement makes same-cluster collision terms visible, and then gives the exact-$m$ time-series analogue, where vertical partitions play the same role. The near-epoch-dependent extension is added next, by approximating the observed process with exact finite-range versions.
\subsection{Clustered sampling}\label{Sect:Cluster}
Given a sampled clustered array $\{W_i\}_{i\in \mathcal{G}_g}, \ g\in [G]$, the following imposes the sampling assumption.
\begin{assumption}[Sampling]\label{ass:sampling}
$\{W_i\}_{i\in \mathcal{G}_g}, \ g\in [G] $ are independently distributed across clusters $g \in [G] $.
\end{assumption}
\noindent \Cref{ass:sampling} leaves within-cluster dependence unrestricted and imposes independence only across clusters. Cluster sizes may be equal or unequal. The framework therefore covers short balanced panels, unbalanced panels, clustered data, cross-sectional data, and repeated cross-sections with arbitrary within-cluster dependence.
The order-$2$ clustered case is the simplest place to see the full mechanism. For intuition, sort the clusters in decreasing order of size so that the number of columns equals the largest cluster size. Suppose the ordered cluster sizes are $4,3,2,1$, and write $W_{g,a}$ for the $a$-th observation in cluster $g$. In the following display, box styles identify clusters and braces identify the vertical partitions:
\begingroup
\tikzset{
clusterbox/.style={draw, rounded corners=2pt, minimum width=1.35cm,
minimum height=0.48cm, inner sep=1pt, align=center},
clusterone/.style={clusterbox, fill=gray!8, draw=black, line width=0.45pt},
clustertwo/.style={clusterbox, fill=gray!22, draw=black, line width=0.45pt},
clusterthree/.style={clusterbox, fill=white, draw=black, dashed, line width=0.55pt},
clusterfour/.style={clusterbox, fill=gray!35, draw=black, dotted, line width=0.65pt}
}
\[
\begin{array}{cccc}
\overbrace{
\begin{array}{c}
\tikz[baseline=(x.base)] \node[clusterone] (x) {$W_{1,1}$};\\[2pt]
\tikz[baseline=(x.base)] \node[clustertwo] (x) {$W_{2,1}$};\\[2pt]
\tikz[baseline=(x.base)] \node[clusterthree] (x) {$W_{3,2}$};
\end{array}}^{\mathcal M_n(1)}
&
\overbrace{
\begin{array}{c}
\tikz[baseline=(x.base)] \node[clusterone] (x) {$W_{1,2}$};\\[2pt]
\tikz[baseline=(x.base)] \node[clustertwo] (x) {$W_{2,2}$};\\[2pt]
\tikz[baseline=(x.base)] \node[clusterfour] (x) {$W_{4,1}$};
\end{array}}^{\mathcal M_n(2)}
&
\overbrace{
\begin{array}{c}
\tikz[baseline=(x.base)] \node[clusterone] (x) {$W_{1,3}$};\\[2pt]
\tikz[baseline=(x.base)] \node[clustertwo] (x) {$W_{2,3}$};\\[2pt]
\varnothing
\end{array}}^{\mathcal M_n(3)}
&
\overbrace{
\begin{array}{c}
\tikz[baseline=(x.base)] \node[clusterone] (x) {$W_{1,4}$};\\[2pt]
\tikz[baseline=(x.base)] \node[clusterthree] (x) {$W_{3,1}$};\\[2pt]
\varnothing
\end{array}}^{\mathcal M_n(4)}
\end{array}.
\]
\endgroup
The useful property is columnwise: each brace covers at most one observation from any cluster, even with unequal cluster sizes. Thus observations read down a column are cross-cluster and independent under \Cref{ass:sampling}. Same-cluster pairs must occur across columns and are counted as same-cluster collision terms. The appendix gives the exact indexing for the general unbalanced array.
The next assumption, which is \citet[Assumption 2]{hansen2019asymptotic}, regulates the heterogeneity in cluster sizes.
\begin{assumption}[Cluster Size Heterogeneity]\label{ass:clus_het} For some $ 2 \leq r < \infty $,
\[
\frac{\Big( \sum_{g=1}^G n_g^r \Big)^{2/r}}{n} \leq C < \infty, \quad \text{ and } \quad \max_{g\in [G]}\frac{n_g^2}{n} \rightarrow 0 \quad \text{as} \quad n\rightarrow \infty.
\]
\end{assumption}
\noindent Let $ \displaystyle m_n:=\max_{g\in [G]}n_g$ denote the maximum cluster size; then \Cref{ass:clus_het} implies $m_n^2/n\to0$, so $m_n/\sqrt n\to0$ and, \emph{a fortiori}, $m_n/n\to0$.
Denote the U-statistic by
\[
U_{n,2}:= \frac{1}{n(n-1)}\sum_{i=1}^n \sum_{j\neq i}^n \varphi(W_i,W_j),
\]
where, without loss of generality, the scalar kernel $\varphi$ is symmetric and centred on cluster-generic pairs, i.e., \( \mathbb{E}[\varphi(W_i,W_j)]=0 \) whenever \(j\notin\mathcal{G}_{g(i)}\). No centring is imposed on within-cluster pairs. Define the first-order projection \( \varphi_j^{(1)}(W_i):= \mathbb{E}[\varphi(W_i,W_j) \mid W_i] \), its cross-cluster average
\( \displaystyle
\bar\varphi_n(W_i):= \frac{1}{n-n_{g(i)}} \sum_{j=1}^{n} \mathbbm{1}\{j\not\in \mathcal{G}_{g(i)} \} \varphi_j^{(1)}(W_i),
\)
and the scalar variance term
\[
\sigma_{n,2}^2:= \frac{4}{n} \sum_{g=1}^G \mathbb{E}\Big[ \Big( \sum_{i\in\mathcal{G}_g} \bar\varphi_n(W_i) \Big)^2 \Big].
\]
\noindent The following result provides a limit theory for $ U_{n,2} $ under clustered sampling.
\begin{theorem}\label{thm:CLT_Ustats_Clus}
Suppose
\begin{equation}\label{eq:clus_ui_2}
\lim_{M\rightarrow \infty} \sup_{i,j \leq n } \mathbb{E}\big[ | \varphi(W_i,W_j) |^r \mathbbm{1}\{ | \varphi(W_i,W_j) | > M \} \big] = 0
\end{equation}
and
\Cref{ass:sampling,ass:clus_het} hold. Then
(a) \( \displaystyle U_{n,2} = \frac{2}{n} \sum_{i=1}^n \bar\varphi_n(W_i) + \mathcal O_p\Big(\frac{m_n}{n}\Big) \);
(b) \( \displaystyle U_{n,2}=\mathcal O_p\Big(\sqrt{\frac{m_n}{n}}\Big)=o_p(1) \).
If, in addition, \( \sigma_{n,2} \geq \lambda > 0 \), then
(c) \( \displaystyle \sigma_{n,2}^{-1}\sqrt{n} U_{n,2} \xrightarrow{d} \mathcal{N}(0,1) \).
\end{theorem}
\noindent The theorem makes the order-$2$ mechanism explicit. Part (a) gives the asymptotically linear, Hájek-type representation, with the non-linear remainder reduced to order $m_n/n$ by counting same-cluster collisions and the degenerate term. Part (b) is then the corresponding weak law, while part (c) applies the clustered central limit theorem to the leading projection under a non-degeneracy condition. When $m_n=1$, i.e., when $n_g=1$ for every cluster, $\mathcal{G}_{g(i)}=\{i\}$ for every $i$, $\bar\varphi_n(W_i)$ reduces to the usual Hájek projection $\frac{1}{n-1}\sum_{j\neq i}\varphi_j^{(1)}(W_i)$ of an iid $U$-statistic, and all three parts recover the classical iid $U$-statistic asymptotics as a special case; see, e.g., \citet[Theorem 1, Sect. 3.2.1]{lee1990ustatistics}.
\subsection[Exact m-dependence]{Exact $m$-dependence}\label{Sect:FixedM}
The clustered order-$2$ argument has a direct temporal analogue. Same-cluster collision terms are replaced by lag-collision terms, and the vertical columns of the clustered rearrangement are replaced by vertical partitions modulo $m+1$. Exact finite-range dependence is the cleanest setting for seeing this analogy because those vertical partitions comprise independent samples. This lag-generic/lag-collision distinction is close in spirit to \citet{sen1963properties}'s distinction between non-serial tuples and tuples containing serial interactions.
The following definition, adapted from \citet[Definition 4.1]{henze2024asymptotic}, formalises $m$-dependence.
\begin{definition}[$m$-Dependence]\label{def:m_dep}
Let $m\in\mathbb{N}_0:=\{0,1,2,\ldots\}$. A sequence $\{W_i\}_{i\geq1}$ is called \emph{$m$-dependent} if the $\sigma$-fields $\sigma(W_1,\ldots,W_s)$ and $\sigma(W_{s+m+j}:j\geq1)$ are independent for every $s\geq1$.
\end{definition}
\noindent Equivalently, any two blocks of the sequence separated by more than $m$ time periods are independent, while dependence within an $m$-neighbourhood is unrestricted. The case $m=0$ recovers an independent sequence $\{W_i\}$. The resulting mutual independence of index sets pairwise separated by more than $m$ is recorded as a standard fact for $m$-dependent sequences by \citet[Sect.~2.2]{janson2023asymptotic}.
The following stationarity assumption is useful in simplifying the exposition and arguments.
\begin{assumption}[Stationarity]\label{ass:stationary}
$\{W_i\}_{i\geq1}$ is stationary: for every $j\geq1$ and $\iota\geq0$, the distribution of $(W_j,\ldots,W_{j+\iota})$ does not depend on $j$.
\end{assumption}
\noindent Extending the argument to heterogeneous, non-stationary $m$-dependence, the temporal analogue of \Cref{ass:clus_het}, would come at the cost of more notation. \Cref{ass:stationary} is therefore imposed throughout the temporal pipeline for ease of exposition.
Let $W_i$, $i\in[n]$, generically denote data satisfying \Cref{def:m_dep,ass:stationary}. Denote the U-statistic by
\[
U_{n,2}:= \frac{1}{n(n-1)}\sum_{i=1}^n \sum_{j\neq i} \varphi(W_i,W_j),
\]
where, without loss of generality, the scalar kernel $\varphi$ is symmetric and centred on lag-generic pairs, i.e., $\mathbb{E}[\varphi(W_i,W_j)]=0$ whenever $|i-j|>m$. No centring is imposed on lag-collision pairs. For any $i,j$ with $|i-j|>m$, \Cref{def:m_dep} gives $W_i\protect\mathpalette{\protect\independenT}{\perp} W_j$, so the first-order projection $\mathbb{E}[\varphi(W_i,W_j)\mid W_i]$ does not depend on which distant co-index $j$ is used. By \Cref{ass:stationary}, it therefore equals a single function
\(
\varphi^{(1)}(W_i):=\mathbb{E}[\varphi(W_i,W_j)\mid W_i]
\)
for every $j$ with $|i-j|>m$, dropping the co-index subscript carried by $\varphi_j^{(1)}(W_i)$ in \Cref{Sect:Cluster}.
Partition $[n]$ into \(\mathcal M_{n,m}(\ell):=\{i\in[n]:i\bmod(m+1)=\ell\}\), \(\ell=0,\ldots,m\). Consecutive elements of each partition are spaced $m+1$ apart and are therefore mutually independent. Call a pair $(i,j)$ a \emph{lag-collision} if $|i-j|\leq m$ and \emph{lag-generic} otherwise; no lag-collision pair lies within a single partition. The following extends \Cref{thm:CLT_Ustats_Clus} to the $ m $-dependence setting.
\begin{theorem}[Fixed-\texorpdfstring{$m$}{m} CLT]\label{thm:CLT_Ustats_mdep}
Suppose
\begin{equation}\label{eq:mdep_ui_2}
\lim_{M\rightarrow\infty}\sup_{i,j\leq n}\mathbb{E}\big[|\varphi(W_i,W_j)|^r\mathbbm{1}\{|\varphi(W_i,W_j)|>M\}\big]=0
\end{equation}
for some $r\geq2$, and $\{W_i\}_{i\geq1}$ satisfies \Cref{def:m_dep} for some fixed $m\in\mathbb N_0$ and \Cref{ass:stationary}. Then
(a) \( \displaystyle U_{n,2}=\frac{2}{n}\sum_{i=1}^n \varphi^{(1)}(W_i)+\mathcal O_p\Big(\frac{m+1}{n}\Big) \);
(b) \( \displaystyle U_{n,2}=\mathcal O_p\Big(\sqrt{\frac{m+1}{n}}\Big)=o_p(1) \).
If, in addition, \(\sigma_{m,2}^2\geq\lambda>0\), with
\[
\sigma_{m,2}^2
:=
4\Big(\mathbb{V}\big[\varphi^{(1)}(W_1)\big]
+2\sum_{h=1}^m\mathrm{Cov}\big(\varphi^{(1)}(W_1),\varphi^{(1)}(W_{1+h})\big)\Big),
\]
then (c) \( \displaystyle \sigma_{m,2}^{-1}\sqrt{n}\,U_{n,2}\xrightarrow{d}\mathcal N(0,1) \).
\end{theorem}
\noindent This is the temporal counterpart of the clustered order-$2$ result. Parts (a)--(b) follow once same-cluster collisions are replaced by lag-collisions: only $\mathcal O(nm)$ of the $\mathcal O(n^2)$ ordered pairs are non-generic, and the degenerate Hoeffding remainder again contributes only $\mathcal O_p(n^{-1})$. The leading term is therefore the first-order projection. Since $\{\varphi^{(1)}(W_i)\}_{i\geq1}$ is itself stationary and $m$-dependent, part (c) follows from the central limit theorem for fixed-$m$ sequences \citep[Eq.\ (4.2), Theorem 4.6]{henze2024asymptotic}. In the exact-$m$ case, that variance has a finite auto-covariance formula. The explicit \(m\)-dependence of these rates is retained because the NED extension below reuses the same bounds along a sequence \(m=m_n\to\infty\).
\subsection{Near-epoch dependence}\label{Sect:NED}
The exact-$m$ result completes the order-$2$ argument when dependence has a fixed range. Near-epoch dependence relaxes that restriction by replacing the observed process with finite-range approximations. In the clustered design, the vertical rearrangement is data-driven through the largest cluster size \(m_n\); under exact \(m\)-dependence, the corresponding partition width is structural and equals \(m+1\); under NED, the sequence \(m_n\) is instead an \emph{approximating-window device} used only to transfer the exact-\(m\) argument. Under exact \(m\)-dependence, the span \(m\) marks a region of unrestricted local dependence and exact independence beyond it. Under NED, the same span appears only in the approximating proxy \(W_i^{(m)}\), which depends on a block of \(m+1\) nearby i.i.d. innovations. Thus the local dependence inside the proxy is structured rather than arbitrary, and the original process is recovered by letting the approximation window widen. For $W_i$ generated by an i.i.d.-driven process, write
\[
\mathcal F_a^b:=\sigma(\mathcal E_t:a\leq t\leq b),\qquad a\leq b\in\mathbb Z\cup\{\pm\infty\},
\]
where $\{\mathcal E_t\}_{t\in\mathbb Z}$ is an i.i.d.\ base sequence and each $W_i$ is $\mathcal F_{-\infty}^{\infty}$-measurable. For $m\in\mathbb N_0$, define
\(
W_i^{(m)}:=\mathbb{E}\big[W_i\mid\mathcal F_{i-\floor{m/2}}^{i+\ceil{m/2}}\big],
\)
where $\mathcal F_{i-\floor{m/2}}^{i+\ceil{m/2}}$ is the sigma-field generated by a retained innovation block of width $m+1$ around date $i$. This convention indexes the approximation by its dependence span, thereby rendering the approximating sequence $\{W_i^{(m)}\}$ exactly $m$-dependent.
The exact-\(m\) theorem can then be used for the approximating statistic after centring the truncated kernel under its generic law,
\[
U_{n,2}^{[m]}:= \frac{1}{n(n-1)}\sum_{i=1}^n \sum_{j\neq i} \varphi(W_i^{(m)},W_j^{(m)}).
\]
The main remaining step is to make the difference between the original and approximating $U$-statistics negligible. Process approximation alone is not enough for this: one also needs a kernel-regularity bridge that transfers the input approximation rate \(\nu_m=\|W_i-W_i^{(m)}\|_2\) into a kernel approximation rate.
This finite-window approximation motivates the following $L_2$-NED condition, adapted to the sequence formulation in \citet[Def.~17.1]{davidson1994stochastic} and specialised to the symmetric retained blocks used here.
\begin{assumption}[Near-epoch dependence of $\{W_i\}$]\label{ass:ned}
There is an i.i.d.\ base sequence $\{\mathcal E_t\}_{t\in\mathbb Z}$ such that each $W_i$ is $\mathcal F_{-\infty}^{\infty}$-measurable, where \(\mathcal F_a^b:=\sigma(\mathcal E_t:a\leq t\leq b)\). Moreover, $\{W_i\}$ is $L_2$-near-epoch dependent on this base sequence, i.e., there is a sequence of constants $\nu_m\to0$ such that
\[
\|W_i-W_i^{(m)}\|_2=\big\|W_i-\mathbb{E}[W_i\mid\mathcal F_{i-\floor{m/2}}^{i+\ceil{m/2}}]\big\|_2\leq\nu_m\quad\text{for every }i\in\mathbb N\text{ and }m\in\mathbb N_0.
\]
\end{assumption}
\noindent Davidson's scaling constants $d_i$ are absorbed into $\nu_m$: since $\big(W_i,W_i^{(m)}\big)$ has, by construction, the same joint law for every $i$ under \Cref{ass:stationary}, the left-hand side does not depend on $i$, so $d_i\equiv1$ entails no loss of generality.
Following \citet[p.~262]{davidson1994stochastic}, $\{W_i\}$ is $L_2$-NED \emph{of size $-\zeta_0$} if $\nu_m=\mathcal O(m^{-\zeta})$ for some $\zeta>\zeta_0$. \Cref{ass:ned} controls the input approximation error \(W_i-W_i^{(m)}\), but the $U$-statistic also requires a kernel-level bridge from input approximation to kernel approximation. Standard kernels exhibit two useful perturbation behaviours. The Gini mean difference kernel \(\varphi(x,y)=|x-y|\) is globally Lipschitz, while rank kernels such as Kendall's tau and Spearman-type sign kernels are not pointwise Lipschitz but can be controlled after taking expectations. These kernels are naturally handled by an expected-H\"older condition. By contrast, the squared-difference kernel satisfies \( |(w-y)^2-(\widetilde w-y)^2|=|w-\widetilde w|\,|w+\widetilde w-2y| \), so its local slope is available but random and unbounded; this is better captured by a data-dependent modulus with tail and moment control. Because both behaviours arise for common $U$-statistic kernels, the next assumption supplies either of two non-nested one-coordinate perturbation bridges. For the data-dependent modulus branch below, given a candidate modulus \(L(w,\widetilde w,w^\star)\), write \(L_{ij}^{(m,1)}:=L(W_i,W_i^{(m)},W_j)\) and \(L_{ij}^{(m,2)}:=L(W_j,W_j^{(m)},W_i^{(m)})\).
\begin{assumption}[Kernel perturbation for the order-$2$ kernel]\label{ass:kernel_perturb_2}
One of the following two branches holds.
\begin{enumerate}[label=(\alph*)]
\item \emph{Expected-H\"older branch.} There exist $C<\infty$ and \(\alpha\in(0,1]\) such that, uniformly over \(i\neq j\) and \(m\in\mathbb N_0\),
\( \displaystyle
\|\varphi(W_i,W_j)-\varphi(W_i^{(m)},W_j)\|_2
+
\|\varphi(W_i^{(m)},W_j)-\varphi(W_i^{(m)},W_j^{(m)})\|_2
\leq C\nu_m^\alpha,
\)
with the same bound understood for the corresponding generic independent copies used in projection arguments. In this branch, set \(\rho_m:=\nu_m^\alpha\).
\item \emph{Data-dependent modulus branch.} There is a measurable function \(L(w,\widetilde w,w^\star)\), taking values in \([0,\infty)\), such that, for all random variables \(W,\widetilde W,W^\star\) defined on the same probability space, \(|\varphi(W,W^\star)-\varphi(\widetilde W,W^\star)|\leq L(W,\widetilde W,W^\star)\|W-\widetilde W\|\). There exists a \(\delta>0\) such that \(\displaystyle \sup_{i\neq j,\,m\in\mathbb N_0,\,a=1,2}\mathbb{E}\big[L_{ij}^{(m,a)}\mathbbm{1}\{L_{ij}^{(m,a)}>M\}\big]=\mathcal O(M^{-\delta})\) as \(M\to\infty\). For some \(r>2\), \(\displaystyle \sup_{i\neq j,\,m\in\mathbb N_0}\big\|\max\{|\varphi(W_i,W_j)|,|\varphi(W_i^{(m)},W_j)|,|\varphi(W_i^{(m)},W_j^{(m)})|\}\big\|_r<\infty\). The same tail and moment bounds are understood for the corresponding generic independent copies used in projection arguments. Set \(\beta:=\tau/(1+\tau)\), where \(\tau:=(1+\delta)(r-2)/(2r)\), and set \(\rho_m:=\nu_m^\beta\).
\end{enumerate}
\end{assumption}
\noindent In either branch, \(\rho_m\) is the effective kernel-approximation rate induced by the NED approximation. Both branches keep the same one-coordinate perturbation format -- only one argument is replaced, while the co-index is read with the relevant independent or mixed original/truncated law. The case \(\alpha=1\) in the expected-H\"older branch is the expected Lipschitz condition along the NED couplings, and a deterministic Lipschitz bound verifies it immediately. Some kernels with unbounded local moduli may still fall into this branch when the relevant NED couplings have enough integrability to absorb the local slope into the finite constant \(C\).\footnote{For example, under common-mean, common-scale Gaussian NED couplings, the non-globally-Lipschitz kernel \(\varphi(x,y)=(x-y)^2\) can satisfy the expected branch along the relevant couplings: with \(D:=W-\widetilde W\), \(A:=W+\widetilde W-2W^\star\), and \(W^\star\protect\mathpalette{\protect\independenT}{\perp}(W,\widetilde W)\), joint normality gives \(D\protect\mathpalette{\protect\independenT}{\perp} A\), so \(\|\varphi(W,W^\star)-\varphi(\widetilde W,W^\star)\|_2=\|D\|_2\|A\|_2\), while \(\|A\|_2\leq2\sqrt2\,\sigma\). Thus \(C=2\sqrt2\,\sigma\) is admissible for those couplings.}
The two branches are complementary. The Gini mean difference sits in the overlap: it satisfies the expected-H\"older branch with \(\alpha=1\), and it also satisfies the data-dependent modulus branch with a constant modulus. For globally Lipschitz kernels, the expected-H\"older route is simpler and sharper. For kernels with random unbounded slopes, including sample variances, squared deviations, and quadratic loss functions, the data-dependent modulus branch can be the more natural route. For sign kernels, the expected branch may be verified by a conditional crossing-probability argument rather than by pointwise continuity.\footnote{A sufficient condition is that the signed distance \(\Delta\) to the relevant threshold has a conditional density, given the perturbation \(D_m\) or the relevant finite-window information, uniformly bounded near zero. Then a sign changes only when \(\Delta\) lies between \(0\) and \(D_m\), so \(\mathbb{P}(\mathrm{crossing}\mid D_m)\leq C|D_m|\), and \(\|\mathrm{sign}(\Delta)-\mathrm{sign}(\Delta-D_m)\|_2=\mathcal O(\nu_m^{1/2})\) whenever \(\|D_m\|_2\leq\nu_m\). Kendall's tau follows directly, while Spearman-type products follow from \(ab-\widetilde a\widetilde b=a(b-\widetilde b)+(a-\widetilde a)\widetilde b\).}
The effective perturbation bound in \Cref{lem:effective_perturbation}, evaluated at \(k=2\), then gives
\(
\big\|U_{n,2}-U_{n,2}^{[m]}\big\|_2\leq C\rho_m.
\)
The proof applies the exact-\(m\) result to this centred approximating kernel and then transfers the conclusion back to the original statistic. To pass from the fixed-\(m\) approximation to the original NED statistic, choose \(m=m_n\) with \(m_n\to\infty\), \(m_n=o(\sqrt n)\), and \(\sqrt n\,\rho_{m_n}=o(1)\). The first condition lets the approximating sequence approach the original one; the second makes the fixed-\(m\) projection remainder negligible after \(\sqrt n\) scaling; the third makes the NED approximation error negligible. The following provides a polynomial rate illustration.
\begin{example}[Kernel--dependence rate balance]\label{ex:rate_balance}
If \(\rho_m=\mathcal O(m^{-\eta})\) and \(m_n\sim n^\gamma\), then \(m_n=o(\sqrt n)\) requires \(\gamma<1/2\), while \(\sqrt n\,\rho_{m_n}=o(1)\) requires \(\gamma\eta>1/2\). Hence the admissible region is \(\gamma\in(1/(2\eta),1/2)\), which is non-empty exactly when \(\eta>1\). In the expected-H\"older branch, \(\eta=\alpha\zeta\) when \(\nu_m=\mathcal O(m^{-\zeta})\); in the data-dependent modulus branch, \(\eta=\beta\zeta\). Thus feasibility is governed jointly by kernel regularity and the finite-window approximation rate of the data. Smoother kernels permit slower NED approximation, whereas rougher kernels require faster NED approximation to compensate. For instance, in the expected-H\"older branch with \(\alpha=1/2\) and \(\zeta=3\), any \(\gamma\in(1/3,1/2)\) is admissible.
\end{example}
All ingredients are now in place for the order-$2$ limiting result under NED data.
In the NED subsections, centring is understood with respect to independent copies of the generic marginal law -- cf. \citet{fischer2016multivariate}.
\begin{theorem}[Order-\texorpdfstring{$2$}{2} NED CLT]\label{thm:CLT_NED_2}
Let $\varphi(w_1,w_2)$ be symmetric. Suppose \Cref{ass:stationary,ass:ned,ass:kernel_perturb_2} hold, the uniform integrability condition \eqref{eq:mdep_ui_2} holds for $\{\varphi(W_i,W_j)\}$ and uniformly for $\{\varphi(W_i^{(m)},W_j^{(m)})\}$, \(\rho_m=\mathcal O(m^{-\eta})\) for some \(\eta>1\), and
\( \displaystyle
\mathbb{V}\Big[ \frac{1}{\sqrt{n}} \sum_{i=1}^n\varphi^{(1)}(W_i)\Big]\to \sigma^2>0.
\)
If $m_n\to\infty$, $m_n=o(\sqrt n)$, and $\sqrt n\,\rho_{m_n}=o(1)$, then
(a) $\displaystyle U_{n,2}=\frac2n\sum_{i=1}^n\varphi^{(1)}(W_i)+o_p(n^{-1/2})$;
(b) $\displaystyle U_{n,2}=\mathcal O_p(n^{-1/2})=o_p(1)$; and
(c) $\displaystyle (2\sigma)^{-1}\sqrt n\,U_{n,2}\xrightarrow{d}\mathcal N(0,1)$.
\end{theorem}
\noindent The same projection logic is still at work. The statistic is reduced to the linear term $2n^{-1}\sum_i\varphi^{(1)}(W_i)$, but now near-epoch dependence supplies the approximation step that replaces the fixed-$m$ argument. The variance is normalised through the asymptotic variance of the linear term.
\begin{remark}
The i.i.d.\ base formulation is used to make the truncation step exact. Once $W_i$ is replaced by $W_i^{(m)}=\mathbb{E}[W_i\mid\mathcal F_{i-\floor{m/2}}^{i+\ceil{m/2}}]$, retained innovation blocks attached to observations more than $m$ dates apart are non-overlapping and therefore independent, so the approximating sequence is exactly $m$-dependent; the formal verification is given in \Cref{app:auxiliary}. The observed process $\{W_i\}$ may still exhibit persistent and non-linear dependence through its underlying i.i.d.-innovation representation.\footnote{For treatments of shift transformations and their role in representing stochastic sequences, see, for example, \citet[Sects.~13.1--13.2]{davidson2021stochastic}.} This formulation covers many standard econometric data-generating processes, including linear processes such as MA($\infty$) and stationary ARMA models, volatility recursions such as ARCH/GARCH-type specifications, and smooth non-linear autoregressions, whenever the usual stability, contraction, and moment conditions deliver finite-window approximation errors, e.g., \citet[Examples 17.3--17.4]{davidson2021stochastic}. A broader NED-on-mixing formulation would require an additional coupling layer, which is not pursued here.
\end{remark}
\section{Order-\texorpdfstring{$k$}{k} Theory}\label{Sect:OrderKTheory}
The order-$2$ pipeline in the preceding section isolates the main ideas: arrange the data so that independent pieces are visible, control the relevant same-cluster or lag-collision terms, reduce the statistic to its first-order projection, and then apply the relevant central limit theorem to that projection. The order-$k$ theory keeps the same structure.
\subsection[Clustered order-k U-statistic]{Clustered order-$k$ U-statistic}\label{Sect:Orderk}
Let $\varphi(w_1,\ldots,w_k)$, $k\geq2$, be symmetric and centred on cluster-generic $k$-tuples, i.e. $\mathbb{E}[\varphi(W_{i_1},\ldots,W_{i_k})]=0$ whenever the indices are distinct and belong to distinct clusters. The order-$k$ $U$-statistic is
\[
U_{n,k}:= \frac{1}{\mathcal{P}(n,k)}\sum_{(i_1,\ldots,i_k)_{\neq}}\varphi(W_{i_1},\ldots,W_{i_k}),
\]
where $\mathcal{P}(n,k):=n(n-1)\cdots(n-k+1)$ denotes the number of permutations of $k$ distinct indices from $[n]$. Define the first-order projection $\varphi^{(1)}_{(i_2,\ldots,i_k)}(W_i):=\mathbb{E}[\varphi(W_i,W_{i_2},\ldots,W_{i_k})\mid W_i]$. Let
\(
Q_{n,i}:=\#\{(i_2,\ldots,i_k)_{\neq}: g(i,i_2,\ldots,i_k)_{\neq}\},
\)
and define the cluster-generic co-index average
\[
\bar\varphi_{n,k}^{(1)}(W_i):=\frac{1}{Q_{n,i}}\sum_{\substack{(i_2,\ldots,i_k)_{\neq}\\ g(i,i_2,\ldots,i_k)_{\neq}}}\varphi^{(1)}_{(i_2,\ldots,i_k)}(W_i),
\]
where the sum ranges over ordered co-index tuples for which \((i,i_2,\ldots,i_k)\) is cluster-generic. Hence \(\mathbb{E}[\bar\varphi_{n,k}^{(1)}(W_i)]=0\). \(Q_{n,i}=n-n_{g(i)}\) and \(\bar\varphi_{n,2}^{(1)}(W_i)=\bar\varphi_n(W_i)\) are recovered at $k=2$. Let $\displaystyle \sigma_{n,k}^2:= \frac{k^2}{n} \sum_{g=1}^G \mathbb{E}\Big[ \Big( \sum_{i\in\mathcal{G}_g} \bar\varphi_{n,k}^{(1)}(W_i) \Big)^2 \Big]$, which reduces to $\sigma_{n,2}^2$ at $k=2$.
\begin{namedtheorem}[Order-\texorpdfstring{$k$}{k} CLT]{$1'$}\label{thm:CLT_Ustats_Clus_k}
Suppose
\begin{equation}\label{eq:clus_ui_k}
\lim_{M\rightarrow \infty} \sup_{\mathbf i}
\mathbb{E}\big[ | \varphi(W_{i_1},\ldots,W_{i_k}) |^r
\mathbbm{1}\{ | \varphi(W_{i_1},\ldots,W_{i_k}) | > M \} \big] = 0
\end{equation}
and \Cref{ass:sampling,ass:clus_het} hold. Then (a)
\( \displaystyle
U_{n,k} = \frac{k}{n} \sum_{i=1}^n \bar\varphi_{n,k}^{(1)}(W_i)
+ \mathcal O_p\Big(\frac{m_n}{n}\Big);
\)
(b) \( \displaystyle U_{n,k}=\mathcal O_p\Big(\sqrt{\frac{m_n}{n}}\Big)=o_p(1) \).
If, in addition, \( \sigma_{n,k} \geq \lambda > 0 \), then
(c) \( \displaystyle \sigma_{n,k}^{-1}\sqrt{n} \, U_{n,k} \xrightarrow{d} \mathcal{N}(0,1) \).
\end{namedtheorem}
\noindent The structure is unchanged from the order-$2$ case: the statistic is linearised by its first-order projection, and the extra work lies in showing that all higher-order collision and degenerate terms remain of smaller order after counting generic and non-generic $k$-tuples.
\subsection[Order-k fixed-m CLT]{Order-$k$ fixed-$m$ CLT}\label{Sect:FixedMk}
The fixed-$m$ order-$k$ theorem is the temporal counterpart of \Cref{thm:CLT_Ustats_Clus_k}. Let $\varphi(w_1,\ldots,w_k)$, $k\geq2$, be symmetric and centred on lag-generic $k$-tuples, i.e., $\mathbb{E}[\varphi(W_{i_1},\ldots,W_{i_k})]=0$ whenever every pair of coordinates is more than $m$ apart, and denote the order-$k$ $U$-statistic by
\[
U_{n,k}:=\frac1{\Perm nk}\sum_{(i_1,\ldots,i_k)_{\neq}}\varphi(W_{i_1},\ldots,W_{i_k}).
\]
Call a $k$-tuple of distinct indices \emph{lag-generic} if every pair of its coordinates is more than $m$ apart, matching the lag-generic/lag-collision distinction in \Cref{Sect:FixedM}. For such tuples, $W_{i_1},\ldots,W_{i_k}$ are mutually independent by \Cref{def:m_dep} (sorting the coordinates, every consecutive gap exceeds $m$). By stationarity, the first-order projection is the same for every lag-generic co-index tuple, i.e., $\mathbb{E}[\varphi(W_i,W_{i_2},\ldots,W_{i_k})\mid W_i]=\varphi^{(1)}(W_i)$ whenever $(i,i_2,\ldots,i_k)$ is lag-generic. Thus the fixed-$m$ linear term uses the single projection $\varphi^{(1)}$, unlike the clustered case that uses the cluster-generic co-index average $\bar\varphi_{n,k}^{(1)}$.
Let $\sigma_{m,k}^2:=k^2\Big(\mathbb{V}\big[\varphi^{(1)}(W_1)\big]+2\sum_{h=1}^m\mathrm{Cov}\big(\varphi^{(1)}(W_1),\varphi^{(1)}(W_{1+h})\big)\Big)$, which reduces to $\sigma_{m,2}^2$ at $k=2$.
\begin{namedtheorem}[Order-\texorpdfstring{$k$}{k} Fixed-\texorpdfstring{$m$}{m} CLT]{$2'$}\label{thm:CLT_Ustats_mdep_k}
Suppose
\begin{equation}\label{eq:mdep_ui_k}
\lim_{M\rightarrow\infty}\sup_{\mathbf i}
\mathbb{E}\big[|\varphi(W_{i_1},\ldots,W_{i_k})|^r
\mathbbm{1}\{|\varphi(W_{i_1},\ldots,W_{i_k})|>M\}\big]=0
\end{equation}
for some $r\geq2$, and $\{W_i\}_{i\geq1}$ satisfies \Cref{def:m_dep} for some fixed $m\in\mathbb N_0$ and \Cref{ass:stationary}. Then (a)
\( \displaystyle
U_{n,k}=\frac kn\sum_{i=1}^n\varphi^{(1)}(W_i)+\mathcal O_p\Big(\frac{m+1}{n}\Big);
\)
(b) \( \displaystyle U_{n,k}=\mathcal O_p\Big(\sqrt{\frac{m+1}{n}}\Big)=o_p(1) \).
If, in addition, $\sigma_{m,k}\geq\lambda>0$, then
(c) \( \displaystyle
\sigma_{m,k}^{-1}\sqrt n\,U_{n,k}
\xrightarrow{d}\mathcal N(0,1)
\).
\end{namedtheorem}
\noindent At $k=2$, \Cref{thm:CLT_Ustats_mdep_k} recovers \Cref{thm:CLT_Ustats_mdep}. The lag-generic/lag-collision split has a classical antecedent in \citet{sen1963properties}, while \citet{malevich_abdalimov_1982} studies stationary $U$-statistics under $m$-dependence in both the fixed-$m$ case and regimes with $m=m(n)\to\infty$. For symmetric $\varphi$, \Cref{thm:CLT_Ustats_mdep_k} is also consistent with the general fixed-$m$ $U$-statistic CLT of \citet[Thm.~3.3]{janson2023asymptotic}.
\subsection[Order-k NED CLT]{Order-$k$ NED CLT}\label{Sect:DoubleLimit}
The NED extension uses the same approximating sequence and rate balance introduced in \Cref{Sect:NED}. The order-$k$ additions are the order-$k$ version of the same two-branch kernel perturbation condition. For the data-dependent modulus branch, given a candidate modulus \(L(w,\widetilde w,w_2^\star,\ldots,w_k^\star)\), write \(L_{\mathbf i}^{(m,a)}\) for \(L\) evaluated at \((W_{i_a},W_{i_a}^{(m)})\), with earlier co-ordinates truncated and later co-ordinates left original, \(a=1,\ldots,k\). Also, for \(a=0,\ldots,k\), write \(Z_{\mathbf i}^{(m,a)}\) for \((W_{i_1},\ldots,W_{i_k})\) with the first \(a\) coordinates replaced by their \(m\)-approximants.
\begin{namedassumption}[Kernel perturbation for the order-$k$ kernel]{$5'$}\label{ass:kernel_perturb_k}
One of the following two branches holds.
\begin{enumerate}[label=(\alph*),itemsep=0pt,topsep=0pt]
\item \emph{Expected-H\"older branch.} There exist $C<\infty$ and \(\alpha\in(0,1]\) such that, uniformly over \(\mathbf i_{\neq}\), \(m\in\mathbb N_0\), and \(1\leq a\leq k\),
\[
\|\varphi(Z_{\mathbf i}^{(m,a-1)})-\varphi(Z_{\mathbf i}^{(m,a)})\|_2
\leq C\nu_m^\alpha,
\]
with the same bound understood for the corresponding generic independent copies used in projection arguments.
In this branch, set \(\rho_m:=\nu_m^\alpha\).
\item \emph{Data-dependent modulus branch.} There is a measurable function
\(L(w,\widetilde w,w_2^\star,\ldots,w_k^\star)\),
taking values in \([0,\infty)\), such that, for all random variables
\(W,\widetilde W,W_2^\star,\ldots,W_k^\star\)
defined on the same probability space,
\[
\big|\varphi(W,W_2^\star,\ldots,W_k^\star)
-\varphi(\widetilde W,W_2^\star,\ldots,W_k^\star)\big|
\leq
L(W,\widetilde W,W_2^\star,\ldots,W_k^\star)\,\|W-\widetilde W\|.
\]
There is $\delta>0$ such that
\( \displaystyle
\sup_{\mathbf i_{\neq},\,m\in\mathbb N_0,\,1\le a\le k}
\mathbb{E}\big[L_{\mathbf i}^{(m,a)}\mathbbm{1}\{L_{\mathbf i}^{(m,a)}>M\}\big]
= \mathcal O(M^{-\delta})\quad\text{as }M\to\infty.
\)
For some \(r>2\),
\begin{equation}\label{eq:random_modulus_kernel_moment_k}
\sup_{\mathbf i_{\neq},\,m\in\mathbb N_0,\,0\le a\le k}
\|\varphi(Z_{\mathbf i}^{(m,a)})\|_r
<\infty,
\end{equation}
The same tail and moment bounds are understood for the corresponding generic independent copies used in projection arguments. Set \(\beta:=\tau/(1+\tau)\), where \(\tau:=(1+\delta)(r-2)/(2r)\), and set \(\rho_m:=\nu_m^\beta\).
\end{enumerate}
\end{namedassumption}
\noindent At \(k=2\), this is the order-$k$ analogue of \Cref{ass:kernel_perturb_2}. It retains the one-coordinate perturbation format in both branches.
The theorem below provides the order-k extension of the NED limit theory.
\begin{namedtheorem}[Order-\texorpdfstring{$k$}{k} NED CLT]{$3'$}\label{thm:CLT_double_limit}
Let $\varphi(w_1,\ldots,w_k)$, $k\geq2$, be symmetric. Suppose that \Cref{ass:stationary,ass:ned,ass:kernel_perturb_k} hold, the uniform integrability condition \eqref{eq:mdep_ui_k} holds for $\{\varphi(W_{i_1},\ldots,W_{i_k})\}$ and uniformly in $m$ for $\{\varphi(W_{i_1}^{(m)},\ldots,W_{i_k}^{(m)})\}$, \(\rho_m=\mathcal O(m^{-\eta})\) for some \(\eta>1\), and
\( \displaystyle
\mathbb{V}\Big[ \frac{1}{\sqrt{n}} \sum_{i=1}^n\varphi^{(1)}(W_i)\Big]\to \sigma^2>0.
\)
If $m_n\to\infty$, $m_n=o(\sqrt n)$, and $\sqrt n\,\rho_{m_n}=o(1)$, then
(a) $\displaystyle U_{n,k}=\frac kn\sum_{i=1}^n\varphi^{(1)}(W_i)+o_p(n^{-1/2})$;
(b) $\displaystyle U_{n,k}=\mathcal O_p(n^{-1/2})=o_p(1)$; and
(c) $\displaystyle (k\sigma)^{-1}\sqrt n\,U_{n,k}\xrightarrow{d}\mathcal N(0,1)$.
\end{namedtheorem}
\noindent \Cref{thm:CLT_double_limit} closes the order-$k$ NED extension: it delivers a CLT for $U_{n,k}$ under near-epoch dependence, with no exact-$m$-dependence requirement on $\{W_i\}$ itself.
\section{Feasible Inference}\label{Sect:FeasibleInference}
The limit theorems above reduce both sampling designs to the covariance of the first-order projection. In the clustered case that covariance is estimated by cluster aggregation, while under near-epoch dependence it is estimated by a long-run covariance procedure.
\subsection{Cluster-robust projection covariance}\label{Sect:ClusterEstimation}
Inference via \Cref{thm:CLT_Ustats_Clus_k} requires an estimator of the covariance of the order-$k$ first-order projection. The scalar target is $\sigma_{n,k}^2$ from \Cref{Sect:Orderk}. For each \(i\in[n]\), define the all-coindex leave-one-out average
\[
\tilde\varphi_{n,k,i}
:=
\frac1{\mathcal{P}(n-1,k-1)}
\sum_{\substack{(i_2,\ldots,i_k)_{\neq}\\ i_a\in[n]\setminus\{i\},\ a=2,\ldots,k}}
\varphi(W_i,W_{i_2},\ldots,W_{i_k}).
\]
Then \(U_{n,k}=n^{-1}\sum_{i=1}^n\tilde\varphi_{n,k,i}\). Centre this observation-level vector by setting
\(
\hat\varphi^{(1)}_{n,k,i}:=\tilde\varphi_{n,k,i}-U_{n,k}.
\)
The natural cluster-robust plug-in estimator is
\[
\hat\sigma_{n,k}^2:=\frac{k^2}{n}\sum_{g=1}^G\Big(\sum_{i\in\mathcal{G}_g}\hat\varphi^{(1)}_{n,k,i}\Big)^2.
\]
\noindent Thus the statistic and covariance estimator are built from the same \(n\)-vector, i.e., the average of \(\tilde\varphi_{n,k,i}\) gives the $U$-statistic, while its centred version gives the cluster summands for feasible inference. The proof compares this all-coindex empirical projection with the infeasible cluster-generic projection after cluster aggregation.
\begin{theorem}[Consistent estimation of $\sigma_{n,k}^2$]\label{thm:cluster_var_consistency}
Under the assumptions of \Cref{thm:CLT_Ustats_Clus_k}, including $\sigma_{n,k}\geq\lambda>0$, $\hat\sigma_{n,k}^2/\sigma_{n,k}^2\xrightarrow{p}1$.
\end{theorem}
\subsection{HAC projection covariance}\label{Sect:Estimation}
\citet{dehling2017testing} considers this problem for Kendall's tau, an order-$2$ $U$-statistic, under near-epoch dependence. The construction combines a plug-in estimator of the first-order projection with a HAC (heteroscedasticity- and autocorrelation-consistent) kernel estimator of the long-run variance, building on \citet{dejong2000consistency}. \citet{fischer2017robust} develops the same broad oracle-plus-plug-in HAC strategy for multivariate kernels in a bounded-kernel NED setting. This section extends that plug-in-plus-HAC architecture to the order-$k$ projection in \Cref{thm:CLT_double_limit}, under the present $L_2$-NED framework and either of the complementary kernel-regularity branches.
For $i=1,\ldots,n$, use the same leave-one-out empirical projection,
\[
\hat\varphi^{(1)}_{n,k,i}:=\tilde\varphi_{n,k,i}-U_{n,k},
\]
as defined in \Cref{Sect:ClusterEstimation}.
Let
\(
\hat\rho(r) := \frac1n\sum_{i=1}^{n-r}\hat\varphi^{(1)}_{n,k,i}\,\hat\varphi^{(1)}_{n,k,i+r},\quad r=0,1,\ldots,n-1
\) be the empirical covariance for lag $r$, and define the HAC-type long-run variance estimator
\[
\hat\sigma_n^2 := \hat\rho(0)+2\sum_{r=1}^{n-1}\kappa\Big(\frac r{b_n}\Big)\hat\rho(r)
\]
where $\kappa(\cdot)$ is a HAC weight function and $b_n$ is a bandwidth sequence.
The following regularity condition follows \citet[Assumption 2.6]{dehling2017testing}.
\begin{assumption}[HAC weights and bandwidth]\label{ass:hac}
$\kappa:\mathbb R\to[-1,1]$ satisfies $\kappa(0)=1$, $\kappa(x)=\kappa(-x)$ for all $x$, $\int|\kappa(x)|\,dx<\infty$, its Fourier transform $f(\xi):=(2\pi)^{-1}\int\kappa(x)e^{i\xi x}\,dx$ satisfies $\int|f(\xi)|\,d\xi<\infty$, and $\kappa$ is continuous at $0$ and at all but finitely many points. The bandwidth satisfies $b_n\to\infty$ and $b_n=o(\sqrt n)$.
\end{assumption}
\noindent \Cref{ass:hac} is satisfied by the Bartlett, Parzen, quadratic-spectral, and Tukey--Hanning weight functions; see \citet{dehling2017testing} and \citet{dejong2000consistency}.
\begin{theorem}[Consistent estimation of $\sigma^2$]\label{thm:hac_consistency}
Under the assumptions of \Cref{thm:CLT_double_limit}, with the uniform integrability condition \eqref{eq:mdep_ui_k} imposed for some \(r>2\), and if \Cref{ass:hac} also holds, $\hat\sigma_n^2\xrightarrow{p}\sigma^2$.
\end{theorem}
\section{Empirical Applications}\label{Sect:Empirical}
The empirical applications mirror the paper's theoretical progression. Each application pairs a familiar order-2 statistic with a genuine order-3 statistic under a common sampling scheme, so the empirical section moves from order-2 intuition to order-$k$ feasible inference in the same spirit as the theory. The county-income application uses the Gini coefficient and L-skewness under state-level clustering, while the financial-return application uses Kendall's tau and Spearman's rho under weak time-series dependence. In both cases, the common issue is feasible inference rather than point estimation: the existing $U$-statistic literature provides only part of what is needed for these applications, and the procedures developed here supply the missing cluster-robust and weak-dependence inference steps used below.
\subsection{County-income inequality and asymmetry}\label{Sect:Empirical_Gini}
\citet{park_shin_2023_county_incomes} studies the distribution dynamics of U.S. county per-capita incomes from 1970 to 2017, including the role of government transfer programmes in compressing cross-county income differences. The application follows that question through two index-based summaries built from the same county panel: one for overall inequality and one for distributional asymmetry. The Gini coefficient is the mean absolute pairwise income gap divided by twice the mean income, so its numerator is a natural order-2 $U$-statistic. Counties are observed repeatedly over time and are naturally grouped within states, making state-clustered inference for annual cross-county inequality a direct use case for \Cref{Sect:ClusterEstimation}. Related work by \citet{bhattacharya2007inference} and \citet{darku2020gini} treats Gini inference under survey designs; the application instead uses administrative county panels and compares pre-transfer and post-transfer county-income measures through both the overall spread and the shape of the cross-county distribution.
Let $Y_{it}$ denote per-capita income in county $i$ in year $t$. For a fixed year $t$, let \(\lambda_{\iota t}\), \(\iota\in[3]\), denote the first three L-moments of the cross-county income distribution, written as
\[
\lambda_{1t}:=\mathbb{E}[\varphi_1(Y_{1t})],\quad
\lambda_{2t}:=\mathbb{E}[\varphi_2(Y_{1t},Y_{2t})],\quad \text{and} \quad
\lambda_{3t}:=\mathbb{E}[\varphi_3(Y_{1t},Y_{2t},Y_{3t})],
\]
with kernels
\(
\varphi_1(y_1):=y_1,\qquad
\varphi_2(y_1,y_2):=\frac12|y_1-y_2|,
\)
and
\[
\varphi_3(y_1,y_2,y_3):=\frac13\{\max(y_1,y_2,y_3)-2\operatorname{med}(y_1,y_2,y_3)+\min(y_1,y_2,y_3)\};
\]
see \citet[eqn. 2.4]{hosking1990lmoments}. The application reports two ratio transforms of adjacent L-moments, namely the Gini coefficient $ \displaystyle G_t:=\frac{\lambda_{2t}}{\lambda_{1t}} $ and the L-skewness $ \displaystyle \tau_{3t}:=\frac{\lambda_{3t}}{\lambda_{2t}} $. The Gini coefficient is the ratio of the second L-moment to the mean, while L-skewness is the ratio of the third to the second L-moment. This gives the inequality and asymmetry summaries a common order-statistic backbone while keeping the second object dimension free. See also \citet{sreelakshmi_asha_nair_2015_lmoments} for related discussion of L-moments and inequality.
The sample analogues are built from an order-1 mean, an order-2 $U$-statistic, and an order-3 $U$-statistic. In particular,
\[
\widehat\lambda_{1t}:=\frac1{n_t}\sum_{i=1}^{n_t}Y_{it},\qquad
\widehat\lambda_{2t}:=\frac1{2n_t(n_t-1)}\sum_{i=1}^{n_t}\sum_{j\neq i}|Y_{it}-Y_{jt}|,
\]
so that \(\widehat G_t=\widehat\lambda_{2t}/\widehat\lambda_{1t}\). Here \(i\in[n_t]\) indexes counties observed in year \(t\), and \(n_t\) is the number of such counties. Likewise, \(\widehat\tau_{3t}=\widehat\lambda_{3t}/\widehat\lambda_{2t}\), where \(\widehat\lambda_{3t}\) is the order-3 $U$-statistic sample analogue of \(\lambda_{3t}\). The clustered first-order projection for \(\widehat G_t\) therefore combines the empirical projections for \((\widehat\lambda_{1t},\widehat\lambda_{2t})\) through the delta method, while that for \(\widehat\tau_{3t}\) combines those for \((\widehat\lambda_{2t},\widehat\lambda_{3t})\). In both cases, the cluster-robust variance sums the resulting first-order projections within states, exactly as in \Cref{Sect:ClusterEstimation}. The reported standard errors therefore treat within-state county outcomes as potentially dependent while using the cross-state variation for inference. This order-3 addition asks a distinct empirical question: not only whether transfers are associated with a compressed overall inequality, but also whether they are associated with altered asymmetry of the county-income distribution.
The data are the replication files from \citet{park_shin_2023_county_incomes}. Following that construction, the pre-transfer series is pre-tax county per-capita income: government transfer income and military income are excluded from Bureau of Economic Analysis (BEA) personal income before dividing by county population. Annual Gini and L-skewness measures are computed from this pre-transfer series and from the post-transfer series that adds all government transfer receipts available in the replication package. The aim is to isolate an index-based implication of that transfer exercise and equip it with inference designed for clustered $U$-statistics. For the Gini, the estimand of interest is
\[
\Delta_t:=G_t^{\mathrm{post\ transfer}}-G_t^{\mathrm{pre\ transfer}},
\]
so negative values indicate that the post-transfer income measure is less unequal across counties than the pre-transfer measure. The corresponding L-skewness gap is defined analogously. \Cref{fig:park_shin_gaps} plots both annual paths. Each panel reports pointwise intervals and simultaneous sup-$t$ bands over calendar years; following the sup-$t$ principle of \citet{li_liao_2020_uniform_nonparametric}.
\begin{figure}[h!]
\centering
\caption{County-income order-2 and order-3 transfer effects, 1970--2017.}
\begin{subfigure}[t]{0.49\textwidth}
\centering
\includegraphics[width=\textwidth]{figures/park_shin_gini_transfer_gap_figure.pdf}
\caption{Gini coefficients and transfer-induced Gini gaps.}
\label{fig:park_shin_gini_gap}
\end{subfigure}\hfill
\begin{subfigure}[t]{0.49\textwidth}
\centering
\includegraphics[width=\textwidth]{figures/park_shin_l3_transfer_gap_figure.pdf}
\caption{L-skewness and transfer-induced gaps.}
\label{fig:park_shin_l3_gap}
\end{subfigure}
\label{fig:park_shin_gaps}
\begin{justify}
{\footnotesize \textit{Notes:} Panel (a) plots the pre-transfer and post-transfer cross-county Gini series in the upper sub-panel and post-transfer minus pre-transfer Gini in the lower sub-panel. Panel (b) plots the analogous L-skewness series and gap. Shaded bands are pointwise 95\% intervals using state-clustered empirical first-order projections; dashed bands are simultaneous 95\% sup-$t$ bands over calendar years, following the same sup-$t$ principle as in \citet{li_liao_2020_uniform_nonparametric}.}
\end{justify}
\end{figure}
The pre-transfer cross-county Gini rises from $0.1495$ in 1970 to $0.1841$ in 2017. By contrast, the post-transfer Gini remains nearly flat, moving from $0.1310$ to $0.1308$ over the same period. The post-transfer gap is already negative in 1970, $-0.0185$, and widens to $-0.0539$ by 2017, with a state-clustered standard error of $0.0029$. L-skewness gives the order-3 counterpart on a scale-free metric, but the signal is much weaker: the post-transfer gap is essentially zero in 1970, at about $-0.0003$, equals about $0.0043$ in 2017, and has a 2017 state-clustered standard error of about $0.0083$. Transfers therefore clearly compress cross-county inequality, whereas the evidence for a systematic effect on right-tail asymmetry is weak once the third L-moment is standardised. This agrees with \citet{park_shin_2023_county_incomes}'s substantive conclusion that transfers mitigate rising county-income inequality, while adding state-clustered uncertainty for annual $U$-statistic summaries. The simultaneous bands sharpen two features of the Gini figure. First, the zero-effect path is rejected uniformly: the upper dashed band for $\Delta_t$ is below zero in every year, with its maximum equal to $-0.0146$. Second, the transfer effect strengthens over time rather than remaining merely negative throughout the sample.
The L-skewness figure is more cautious at the order-3 level: the pre-transfer and post-transfer paths remain close, the gap changes sign over the sample, and the simultaneous bands do not support a uniform asymmetry effect analogous to the Gini result. Taken together, the two figures show that transfers compress the cross-county distribution, while the evidence for a systematic effect on distributional asymmetry is weaker once uncertainty is evaluated under arbitrary within-state dependence.
\subsection{Rank correlations in standardised financial returns}\label{Sect:Empirical_Tong_Hansen}
\citet{tong_hansen_2026_dynamic_factor_correlations} studies dynamic factor correlations in daily standardised returns for U.S. equities and factor portfolios. That empirical design is well suited to the weak-dependence side of the paper. The observations form a long daily time series, while the substantive object is dependence among factors, sectors, and stock returns. The application turns from linear dependence to rank dependence and again pairs an order-2 statistic with an order-3 statistic under the same sampling scheme. Kendall's tau gives the order-2 benchmark related to \citet{dehling2017testing}, while Spearman's rho provides the higher-order counterpart through its order-3 $U$-statistic representation.
The data in \citet{tong_hansen_2026_dynamic_factor_correlations} contain $4,278$ daily observations from January 3, 2007 to December 29, 2023. The variables are AR(1)-EGARCH standardised returns for the Fama--French five factors, the momentum factor, nine sector SPDR ETF factors, and $323$ individual stocks. The cited paper models dynamic linear correlation and sparse idiosyncratic correlation matrices. The exercise instead tracks rank dependence and equips it with HAC inference under weak serial dependence.
For the rank-correlation extension, let $Z_t=(X_t,Y_t)$ denote a bivariate standardised-return series, where $X_t$ is the market factor and $Y_t$ is a sector factor. Write
\[
\tau_K:=\mathbb{E}[\psi_2(Z_1,Z_2)]
\quad\text{with}\quad
\psi_2(z_1,z_2):=\mathrm{sign}(x_1-x_2)\mathrm{sign}(y_1-y_2),
\]
for Kendall's tau, and
\[
\rho_S:=\mathbb{E}[\psi_3(Z_1,Z_2,Z_3)]
\quad\text{with}\quad
\psi_3(z_1,z_2,z_3):=\frac12\sum_{\pi\in S_3}
\mathrm{sign}(x_{\pi(1)}-x_{\pi(2)})
\mathrm{sign}(y_{\pi(1)}-y_{\pi(3)}),
\]
for Spearman's rho, where \(S_3\) denotes the six permutations of \(\{1,2,3\}\). This is the symmetrised order-3 kernel associated with the basic non-symmetric term \(\mathrm{sign}(x_1-x_2)\mathrm{sign}(y_1-y_3)\). Under bounded conditional crossing densities for the relevant pairwise differences, the crossing-probability argument following \Cref{ass:kernel_perturb_2} verifies the expected-H\"older branch with \(\alpha=1/2\) for these rank kernels. For both statistics, standard errors come from the empirical first-order projection and a Bartlett HAC long-run variance estimator. For an order-$k$ statistic, \(\widehat{\mathrm{Var}}(U_{n,k})=k^2\hat\sigma_n^2/n\), so \(\mathrm{SE}(U_{n,k})=k\hat\sigma_n/\sqrt n\).\footnote{Throughout the paper, HAC implementations use the Bartlett weight function \(\kappa(x)=(1-|x|)\mathbbm{1}\{|x|\leq1\}\) with bandwidth \(b_n=\lfloor n^{1/4}\rfloor+1\).}
\begin{table}[h!]
\centering
\caption{Market-sector rank dependence in Tong--Hansen standardised returns, 2007--2023.}
\label{tab:tong_hansen_rank_hac}
\small
\begin{tabular}{lrrrrr}
\toprule
Sector & Pearson & Kendall $\hat\tau$ & SE & Spearman $\hat\rho_3$ & SE \\
\midrule
XLY & 0.8896 & 0.7032 & 0.0064 & 0.8781 & 0.0053 \\
XLK & 0.8906 & 0.7034 & 0.0067 & 0.8754 & 0.0057 \\
XLI & 0.8779 & 0.6831 & 0.0078 & 0.8589 & 0.0068 \\
XLF & 0.8383 & 0.6421 & 0.0086 & 0.8225 & 0.0083 \\
XLB & 0.8242 & 0.6206 & 0.0078 & 0.8069 & 0.0076 \\
XLV & 0.7661 & 0.5548 & 0.0097 & 0.7379 & 0.0106 \\
XLE & 0.6664 & 0.4787 & 0.0116 & 0.6560 & 0.0138 \\
XLP & 0.6769 & 0.4705 & 0.0108 & 0.6457 & 0.0129 \\
XLU & 0.4861 & 0.3179 & 0.0128 & 0.4528 & 0.0172 \\
\bottomrule
\end{tabular}
\begin{justify}
{\footnotesize \textit{Notes:} Each sector series is paired with the market factor \( \mathrm{MKT}\text{-}\mathrm{RF}_t \). Standard errors use empirical first-order projections and Bartlett HAC long-run variances with bandwidth \(b_n=\lfloor n^{1/4}\rfloor+1\).}
\end{justify}
\end{table}
\Cref{tab:tong_hansen_rank_hac} shows that rank dependence is strongest between the market factor and Consumer Discretionary (XLY), Information Technology (XLK), Industrials (XLI), Financials (XLF), and Materials (XLB). The weakest market-sector link is Utilities (XLU), but it remains strongly positive. The Spearman order-3 estimates are numerically almost identical to conventional rank Spearman correlations, with a maximum full-sample gap of about $0.00013$ across the market-factor pairs, confirming that the higher-order representation recovers the familiar rank-correlation object. The empirical payoff is feasible HAC inference for rank-based $U$-statistic dependence measures in a weakly dependent daily time series, now for an order-3 statistic alongside the order-2 Kendall benchmark treated by \citet{dehling2017testing}.
\section{Monte Carlo Simulations}\label{Sect:Simulations}
This section reports two Monte Carlo studies calibrated to the empirical applications in \Cref{Sect:Empirical_Gini,Sect:Empirical_Tong_Hansen}. The clustered arm uses the order-2 Gini mean difference and the order-3 third L-moment, while the weak-dependence arm uses the order-2 Kendall statistic and the order-3 Spearman statistic. In each case, the goal is the same: to assess the corresponding feasible variance estimator through finite-sample dispersion, standard-error calibration, and test size. Both simulation arms use \(2000\) replications. For each design, the reported diagnostics are the normalised median absolute deviation (MAD, or MAD/IQR in the clustered arm), the standard-error to Monte Carlo standard-deviation ratio (SE/SD), and the empirical \(5\%\) rejection frequency.
\subsection{Clustered sampling: Gini mean difference and the third L-moment}\label{Sect:Sim_Cluster}
The clustered simulation uses the same order-2/order-3 kernel pair as the county-income application. The order-2 Gini mean-difference kernel is centred as
\[
\varphi_2(w_1,w_2):=|w_1-w_2|-\theta_{\mathrm{GMD}},\qquad \theta_{\mathrm{GMD}}:=\mathbb{E}|W_1-W_2|,
\]
and the order-3 third L-moment kernel is centred as
\[
\varphi_3(w_1,w_2,w_3):=\frac13\{\max(w_1,w_2,w_3)-2\operatorname{med}(w_1,w_2,w_3)+\min(w_1,w_2,w_3)\}-\theta_{L3},
\]
where \(\theta_{L3}:=\mathbb{E}[\{\max(W_1,W_2,W_3)-2\operatorname{med}(W_1,W_2,W_3)+\min(W_1,W_2,W_3)\}/3]\), with the \(W_a\)'s independent draws from the calibrated marginal law. The DGP is calibrated to the county-income application. Cluster sizes are chosen deterministically from the empirical cluster-size distribution subject to the growth restriction on \(m_n\) in \Cref{ass:clus_het}, and the outcome is generated by applying the empirical quantile function to
\(
\Phi\{\sqrt\rho\,Z_{g(i)}+\sqrt{1-\rho}\,\epsilon_i\},\quad Z_g,\epsilon_i\overset{\mathrm{iid}}\sim\mathcal N(0,1),
\)
where \(\Phi\) is the standard-normal distribution function. Thus \(\rho\) tunes within-cluster dependence while preserving the calibrated marginal distribution. The resulting clusters are heterogeneous and unbalanced, with \(m_n=o(\sqrt n)\), and hence satisfy \Cref{ass:clus_het}.
The design parameters are \(n\in\{400,800,1200\}\) and \(\rho\in\{0,0.25,0.5\}\). The population targets are \(\theta_{\mathrm{GMD}}\) and \(\theta_{L3}\), computed exactly under the calibrated marginal distribution. Because each kernel is centred at its target, the object under study is the centred \(U\)-statistic, its cluster-robust standard error, and the associated cluster-robust \(t\)-statistic for testing the known null value \(0\). To show what is lost by ignoring the sampling scheme, \Cref{tab:sim_cluster_coverage} also reports the corresponding iid standard-error ratio and iid rejection rate alongside the cluster-robust diagnostics.
\begin{table}[h!]
\centering
\caption{Clustered arm: finite-sample behaviour of the cluster-robust \(t\)-statistic.}
\label{tab:sim_cluster_coverage}
\scriptsize
\resizebox{\textwidth}{!}{
\begin{tabular}{lrrrrrrrrrrr}
\toprule
\multicolumn{2}{c}{Design} & \multicolumn{5}{c}{GMD} & \multicolumn{5}{c}{Third L-moment} \\
\cmidrule(lr){1-2}\cmidrule(lr){3-7}\cmidrule(lr){8-12}
$(n,G_n,m_n)$ & $\rho$ & $\frac{\mathrm{MAD}}{\mathrm{IQR}}$ & $\frac{\mathrm{SE}}{\mathrm{SD}}_{\mathrm{rob}}$ & 5\%$_{\mathrm{rob}}$ & $\frac{\mathrm{SE}}{\mathrm{SD}}_{\mathrm{iid}}$ & 5\%$_{\mathrm{iid}}$ & $\frac{\mathrm{MAD}}{\mathrm{IQR}}$ & $\frac{\mathrm{SE}}{\mathrm{SD}}_{\mathrm{rob}}$ & 5\%$_{\mathrm{rob}}$ & $\frac{\mathrm{SE}}{\mathrm{SD}}_{\mathrm{iid}}$ & 5\%$_{\mathrm{iid}}$ \\
\midrule
$(400,163,7)$ & 0.00 & 0.039 & 0.959 & 0.080 & 0.962 & 0.077 & 0.016 & 0.930 & 0.102 & 0.932 & 0.102 \\
& 0.25 & 0.045 & 0.939 & 0.086 & 0.881 & 0.099 & 0.017 & 0.897 & 0.111 & 0.875 & 0.108 \\
& 0.50 & 0.052 & 0.945 & 0.090 & 0.784 & 0.140 & 0.018 & 0.909 & 0.119 & 0.815 & 0.141 \\
\midrule
$(800,294,8)$ & 0.00 & 0.029 & 0.988 & 0.064 & 0.990 & 0.066 & 0.012 & 0.959 & 0.085 & 0.960 & 0.083 \\
& 0.25 & 0.032 & 0.948 & 0.081 & 0.875 & 0.094 & 0.012 & 0.925 & 0.093 & 0.895 & 0.095 \\
& 0.50 & 0.036 & 0.955 & 0.083 & 0.761 & 0.142 & 0.013 & 0.921 & 0.102 & 0.796 & 0.130 \\
\midrule
$(1200,415,9)$ & 0.00 & 0.024 & 0.987 & 0.059 & 0.990 & 0.058 & 0.009 & 0.990 & 0.074 & 0.992 & 0.073 \\
& 0.25 & 0.026 & 0.978 & 0.065 & 0.896 & 0.089 & 0.010 & 0.975 & 0.075 & 0.937 & 0.086 \\
& 0.50 & 0.030 & 0.975 & 0.068 & 0.759 & 0.145 & 0.011 & 0.940 & 0.088 & 0.796 & 0.129 \\
\bottomrule
\end{tabular}
}
\begin{justify}
{\footnotesize \textit{Notes:} The design tuple is \((n,G_n,m_n)\), where \(n\) is the sample size, \(G_n\) is the number of clusters, and \(m_n=\max_{g\leq G_n}n_g\) is the maximum cluster size. GMD is the centred Gini mean difference (order 2) and the third L-moment is the centred order-3 kernel described in \Cref{Sect:Empirical}. \(\rho\) is the one-factor within-cluster dependence parameter. ``MAD/IQR'' is the Monte Carlo median absolute deviation of the centred statistic across replications, divided by the population interquartile range of the empirical county-income calibration, \(\mathrm{IQR}(W)=11{,}587.154\); original-unit MAD values are recovered by multiplying by this IQR. ``SE/SD'' is the ratio of the mean feasible standard error to the Monte Carlo standard deviation, with subscripts ``rob'' and ``iid'' denoting the cluster-robust and iid formulas. Each row uses \(2000\) replications; ``5\%'' is the empirical rejection rate of the corresponding two-sided \(5\%\) \(t\)-test.}
\end{justify}
\end{table}
\Cref{tab:sim_cluster_coverage} shows a clear separation between the robust and iid procedures once within-cluster dependence is present. Under the cluster-robust formula, the SE/SD ratios remain close to one, ranging from \(0.939\) to \(0.988\) for the Gini mean difference and from \(0.897\) to \(0.990\) for the third L-moment, while the corresponding \(5\%\) rejection rates improve steadily with \(n\). Even under \(\rho=0\), where the robust and iid procedures target the same population variance, finite-sample differences remain because the robust estimator is cluster-aggregated. The distortion is most visible for the order-3 L-moment: its cluster-robust rejection rate is \(10.2\%\) at \(n=400\), then falls to \(8.5\%\) and \(7.4\%\) as \(n\) increases, consistent with a more demanding third-order statistic and a moderate effective number of clusters. The iid comparator looks similar only when \(\rho=0\). Once \(\rho>0\), it understates dispersion, with SE/SD ratios falling as low as \(0.759\) for the Gini mean difference and \(0.796\) for the third L-moment, and its rejection rates rise to about \(14\%\). As expected, the order-3 statistic is the more demanding object, but the cluster-robust procedure remains well behaved across the calibrated designs.
\subsection{Near-epoch dependence: Kendall and Spearman rank correlations}\label{Sect:Sim_NED}
The weak-dependence simulation uses the same order-2/order-3 rank-correlation pair as the Tong--Hansen application. Kendall's tau is the order-2 statistic, while Spearman's rho is computed through its order-3 \(U\)-statistic representation.
For each replication, a bivariate stationary Gaussian AR(1) process \(W_t:=(X_t,Y_t)\) is generated through
\(
X_t=\rho X_{t-1}+\sqrt{1-\rho^2}\,u_t,\quad
Y_t=\rho Y_{t-1}+\sqrt{1-\rho^2}\,v_t,
\)
where \((u_t,v_t)\) are i.i.d.\ bivariate normal innovations with correlation \(\alpha\). Thus \(\rho\) tunes serial dependence here, in direct parallel with the clustered design. The dependence levels are calibrated to the full-sample Tong--Hansen market-sector estimates from \Cref{Sect:Empirical_Tong_Hansen}: high dependence uses MKT--XLK, middle dependence uses MKT--XLE, and low dependence uses MKT--XLU. The design parameters are \(n\in\{400,800,1200\}\) and \(\rho\in\{0,0.25,0.50\}\). The null value being tested is the known Gaussian-copula population rank correlation,
\[
\tau(\alpha)=\frac{2}{\pi}\arcsin(\alpha),\qquad
\rho_S(\alpha)=\frac{6}{\pi}\arcsin(\alpha/2),
\]
so the reported rejection rates are size checks for the HAC-based Wald statistic rather than power calculations against independence. As in the clustered arm, \Cref{tab:sim_ned_rank} also includes an iid comparator to show the effect of ignoring serial dependence.
\begin{table}[H]
\centering
\caption{NED rank-correlation arm: finite-sample behaviour of the HAC Wald statistic.}
\label{tab:sim_ned_rank}
\scriptsize
\begin{tabular}{lrrrrrrrrrrrr}
\toprule
& & & & \multicolumn{4}{c}{Kendall} & \multicolumn{4}{c}{Spearman} \\
\cmidrule(lr){5-8}\cmidrule(lr){9-12}
$n$ & Scenario & $\rho$ & $\alpha$ & $\mathrm{MAD}$ & $\frac{\mathrm{SE}}{\mathrm{SD}}_{\mathrm{rob}}$ & 5\%$_{\mathrm{rob}}$ & 5\%$_{\mathrm{iid}}$ & $\mathrm{MAD}$ & $\frac{\mathrm{SE}}{\mathrm{SD}}_{\mathrm{rob}}$ & 5\%$_{\mathrm{rob}}$ & 5\%$_{\mathrm{iid}}$ \\
\midrule
400 & High & 0.00 & 0.891 & 0.0101 & 1.0148 & 0.060 & 0.053 & 0.0083 & 1.0084 & 0.065 & 0.059 \\
& High & 0.25 & 0.891 & 0.0113 & 1.0094 & 0.049 & 0.059 & 0.0089 & 0.9956 & 0.062 & 0.068 \\
& High & 0.50 & 0.891 & 0.0134 & 0.9527 & 0.069 & 0.111 & 0.0107 & 0.9407 & 0.073 & 0.114 \\
& Middle & 0.00 & 0.666 & 0.0173 & 0.9925 & 0.051 & 0.047 & 0.0209 & 0.9897 & 0.053 & 0.050 \\
& Middle & 0.25 & 0.666 & 0.0185 & 0.9631 & 0.067 & 0.070 & 0.0225 & 0.9576 & 0.070 & 0.070 \\
& Middle & 0.50 & 0.666 & 0.0213 & 0.9435 & 0.071 & 0.130 & 0.0259 & 0.9451 & 0.072 & 0.128 \\
& Low & 0.00 & 0.486 & 0.0196 & 0.9767 & 0.062 & 0.059 & 0.0270 & 0.9715 & 0.068 & 0.065 \\
& Low & 0.25 & 0.486 & 0.0214 & 0.9720 & 0.059 & 0.065 & 0.0290 & 0.9672 & 0.062 & 0.073 \\
& Low & 0.50 & 0.486 & 0.0253 & 0.9366 & 0.075 & 0.132 & 0.0349 & 0.9309 & 0.078 & 0.127 \\
\midrule
800 & High & 0.00 & 0.891 & 0.0080 & 0.9782 & 0.058 & 0.055 & 0.0064 & 0.9737 & 0.062 & 0.059 \\
& High & 0.25 & 0.891 & 0.0077 & 1.0106 & 0.050 & 0.055 & 0.0063 & 0.9992 & 0.048 & 0.051 \\
& High & 0.50 & 0.891 & 0.0093 & 0.9545 & 0.069 & 0.114 & 0.0073 & 0.9471 & 0.070 & 0.116 \\
& Middle & 0.00 & 0.666 & 0.0124 & 1.0039 & 0.045 & 0.047 & 0.0152 & 0.9965 & 0.051 & 0.049 \\
& Middle & 0.25 & 0.666 & 0.0128 & 1.0045 & 0.048 & 0.056 & 0.0156 & 1.0035 & 0.048 & 0.057 \\
& Middle & 0.50 & 0.666 & 0.0151 & 0.9586 & 0.066 & 0.121 & 0.0183 & 0.9588 & 0.071 & 0.116 \\
& Low & 0.00 & 0.486 & 0.0142 & 0.9974 & 0.054 & 0.055 & 0.0193 & 0.9960 & 0.058 & 0.058 \\
& Low & 0.25 & 0.486 & 0.0152 & 0.9814 & 0.060 & 0.067 & 0.0206 & 0.9835 & 0.059 & 0.068 \\
& Low & 0.50 & 0.486 & 0.0184 & 0.9451 & 0.061 & 0.119 & 0.0253 & 0.9423 & 0.060 & 0.123 \\
\midrule
1200 & High & 0.00 & 0.891 & 0.0062 & 0.9812 & 0.057 & 0.052 & 0.0050 & 0.9797 & 0.058 & 0.057 \\
& High & 0.25 & 0.891 & 0.0063 & 0.9955 & 0.056 & 0.066 & 0.0050 & 0.9916 & 0.061 & 0.072 \\
& High & 0.50 & 0.891 & 0.0076 & 0.9370 & 0.060 & 0.130 & 0.0060 & 0.9383 & 0.061 & 0.120 \\
& Middle & 0.00 & 0.666 & 0.0103 & 0.9901 & 0.049 & 0.047 & 0.0125 & 0.9914 & 0.049 & 0.049 \\
& Middle & 0.25 & 0.666 & 0.0102 & 0.9984 & 0.054 & 0.061 & 0.0122 & 0.9933 & 0.054 & 0.064 \\
& Middle & 0.50 & 0.666 & 0.0129 & 0.9320 & 0.067 & 0.124 & 0.0155 & 0.9309 & 0.070 & 0.126 \\
& Low & 0.00 & 0.486 & 0.0118 & 0.9839 & 0.057 & 0.056 & 0.0161 & 0.9822 & 0.056 & 0.053 \\
& Low & 0.25 & 0.486 & 0.0120 & 0.9986 & 0.049 & 0.059 & 0.0164 & 0.9984 & 0.050 & 0.060 \\
& Low & 0.50 & 0.486 & 0.0147 & 0.9470 & 0.065 & 0.121 & 0.0201 & 0.9476 & 0.065 & 0.117 \\
\bottomrule
\end{tabular}
\begin{justify}
{\footnotesize \textit{Notes:} Kendall is Kendall's tau (order 2) and Spearman is Spearman's rho (order 3). High, middle, and low dependence use MKT--XLK, MKT--XLE, and MKT--XLU, respectively; $\rho$ is the AR(1) persistence parameter and $\alpha$ is the calibrated contemporaneous correlation. Each row uses $2000$ replications; ``MAD'' is the Monte Carlo median absolute deviation of the centred statistic across replications, ``SE/SD'' is the ratio of the mean feasible standard error to the Monte Carlo standard deviation, and ``5\%'' is the empirical rejection rate of the corresponding two-sided $5\%$ Wald test. Subscripts ``rob'' and ``iid'' denote the HAC and iid formulas.}
\end{justify}
\end{table}
\Cref{tab:sim_ned_rank} tells the same story for weak dependence. The HAC estimator remains well calibrated across dependence levels, with SE/SD$_{\mathrm{rob}}$ ratios between \(0.932\) and \(1.015\) for Kendall and between \(0.931\) and \(1.008\) for Spearman. The robust \(5\%\) rejection rates mostly stay in the \(0.045\) to \(0.078\) range, whereas the iid comparator again breaks down in the persistent designs: once \(\rho=0.5\), iid rejection rates move into the \(0.111\) to \(0.132\) range at \(n=400\) and remain clearly oversized even at \(n=1200\). Dispersion rises, as expected, in the smaller and more persistent designs, with the order-3 Spearman statistic again the more demanding object, but the HAC procedure remains stable across the calibrated designs.
\section{Conclusion}\label{Sect:Conclusion}
This paper develops a unified asymptotic theory and feasible inference framework for order-$k$ $U$-statistics under clustered and weakly dependent sampling. For each sampling scheme, the results include asymptotically linear representations of the U-statistic, a weak law of large numbers, a central limit theorem, and consistent variance estimation. This makes inference feasible for parameters of economic interest often cast as U-statistics, such as Gini coefficients and skewness parameters, under clustered and weakly dependent sampling schemes. The central organising idea is simple: after a vertical rearrangement, the difficult non-linear remainder can be analysed through sampling-generic tuples of pairwise-independent units, while same-cluster and lag collisions are counted and controlled explicitly. That step is common across sampling schemes; what remains sampling scheme-specific is the probability theory for the first-order projection and the way its covariance is estimated. The empirical and simulation results illustrate that this separation leads to usable cluster-robust and HAC procedures for both order-2 and higher-order statistics.
\vspace{0.35cm}
\noindent \texttt{Replication files:}
The replication package, including code, data, and output files, is available on the author's website.
\noindent \texttt{Declaration of AI use:}
During the preparation of this manuscript, the author used OpenAI's ChatGPT and Codex for research assistance. The author reviewed and edited all AI-assisted output and takes full responsibility for the content of the manuscript.
\printbibliography
\end{refsection}
\newpage
\setcounter{page}{1}
\begin{refsection}