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.
77,949 characters · 18 sections · 88 citation commands
Normal Approximation for U-Statistics with Cross-Sectional Dependence
\onehalfspacing
Cross-sectional dependence is pervasive in economic data. Unlike time series data, where observations follow a natural temporal ordering, cross-sectional data consist of dependent observations without an inherent ordering, as individuals may be interconnected in complex ways. In this paper, we study the asymptotic normal approximation for the distribution of an important class of statistics, known as second-order U-statistics, based on a cross-sectionally dependent sample \(\left\{X_{i}: i \in \mathcal{I}_{n}\right\}\). A second-order U-statistic can be written as
In the definition, \(H_{n}:\mathbb{R}^{d}\times \mathbb{R}^{d} \to \mathbb{R}\) is a symmetric measurable function known as the kernel of the U-statistic \(S_{n}\). With an independent and identically distributed (i.i.d.) sample and a fixed permutation symmetric kernel \(H\), this leads to an unbiased estimator of \(\theta = \operatorname{\mathrm{E}} H(X_{1},X_{2})\) after normalisation by the inverse of \(n(n-1)\). Therefore U-statistics can be thought of as a generalisation of sample averages and many important statistics are U-statistics or can be approximated by U-statistics serfling1980ApproximationTheorems. Examples include the sample variance estimators, Kendall's \(\tau\), and Cramer-von-Mises statistics vandervaart2000AsymptoticStatistics. Another motivating example is the nonparametric specification tests as in fan1996ConsistentModel and li2007NonparametricEconometrics for regression models, which will be studied in more details in Section (ref).
Although central limit theorems have been developed for U-statistics with independent or time series samples, existing theories do not apply under cross-sectional dependence and previous methods of proof are no longer applicable. This paper seeks to fill the gap and show conditions under which a normalized U-statistic \(S_n\) based on cross-sectionally dependent data converges in distribution to a normal random variable.
The difficulty of extending existing results on U-statistics to allow for cross-sectional dependence lies in the need to handle both the dependence induced by the nonlinear structure of U-statistics and the cross-sectional dependence in the underlying sample, where the lack of an ordering prohibits the use of classical methods such as the martingale method. In the following, we discuss these challenges and how we approach them.
The structure of U-statistics induces dependence among the summands due to the nonlinearity in the kernel \(H_n\). The classical approach to handle this challenge is to apply the famous Hoeffding decomposition. Let \(\hat{H}_{k}(x) = \operatorname{\mathrm{E}} H(X_{k}, x)\) and \(\theta_{ik} = \operatorname{\mathrm{E}} \hat{H}_{k}(X_{i}) = \operatorname{\mathrm{E}} \hat{H}_{i}(X_{k})\), we can write \(S_{n}\) as
where \(\hat{S}_{n} = 2 \sum_{i}\sum_{k\neq i} (\hat{H}_{k}(X_{i}) - \theta_{ik})\), and \({S}^{*}_{n} = \sum_i \sum_{k\neq i} \left(H(X_{i}, X_{k}) - \hat{H}_{i}(X_{k}) - \hat{H}_{k}(X_{i}) + \theta_{ik}\right).\) In the decomposition, \(\hat{S}_{n}\) is referred to as the projection part and \(S^{*}_{n}\) the degenerate part of the U-statistic \(S_{n}\). We say that a U-statistic is itself degenerate if the first-order projection terms vanish, that is, \(\hat{H}_{k}(X_{i}) - \theta_{ik} = 0\) almost surely for all distinct \(i,k \in \mathcal{I}_{n}\).
When the U-statistic of interest is non-degenerate, the Hoeffding decomposition reduces it to a sum of random variables with an error term \(S^{*}_{n}\). It is then straightforward to obtain central limit theorems with Slutsky's Lemma by demonstrating that, under suitable normalisation, the projection part \(\hat{S}_{n}\) converges in law to a normal random variable, while the degenerate part \(S^{*}_{n}\) is asymptotically negligible.
However, many useful statistics are degenerate in nature, such as those appearing in nonparametric specification tests (Section (ref)). For degenerate U-statistics, the analysis is more involved, since the linear approximation based on projection is now unavailable. In general the limiting distributions for degenerate U-statistics are Gaussian Chaos vandervaart2000AsymptoticStatistics and, in the case of U-statistics of order \(2\), it is a mixture of \(\chi^{2}\) random variables with weights depending on the spectrum of the kernel function.
It is of interest to develop conditions under which asymptotic normality can be recovered for a degenerate U-statistics. One of the classic results due to hall1984CentralLimit finds the following sufficient conditions for CLT of degenerate U-statistics with a random sample of independent and identically distributed random variables \(\left\{X_{i}: 1 \leq i \leq n\right\}\). If the kernel function \(H_n\) has bounded second moment for each $n$, and
where \(\Gamma_n(x, y) = \operatorname{\mathrm{E}}\left[H_n(X_{1}, x) H_n(X_{1}, y)\right]\), then a degenerate U-statistics \(S_{n}\) with symmetric and centered kernel \(H_n\) is asymptotically normally distributed, with mean $0$ and variance $2n^2 \operatorname{\mathrm{E}} H_n^2(X_1,X_2)$. These conditions will be a benchmark for our results on the degenerate U-statistics.
It is worth noting that hall1984CentralLimit's proof relies on martingale representation of \(S_{n}\). Generalisations to time series settings, such as fan1999CentralLimit, also rely on the ordering of the observations which are not available in the cross-sectional dependence settings. We now discuss how we model the cross-sectional dependence and how to handle this additional complication.
Cross-sectional dependence in the underlying sample \(\left\{X_{i}: i\in \mathcal{I}_{n}\right\}\) is a common phenomenon in economic data. There has been a growing interest in modelling cross-sectional dependence as well as in the estimation and inference under cross-sectional dependence. bolthausen1982CentralLimit studies central limit theorems for stationary and mixing random fields, where the random variables are indexed by their location on the lattices \(\mathbb{Z}^{d}\). anselin1988SpatialEconometrics provides review on the spatial statistics. conley1999GMMEstimation further studies GMM estimation and inference when there exists an economic distance among the observations, which is allowed to be imperfectly measured. jenish2009CentralLimit and jenish2012SpatialProcesses extend the previous results to allow for more general dependence structure, nonstationarity and potentially unbounded moments in the random fields. There is also a long tradition of studying dependency graph in the literature of Stein's methods, see for example, chen1986RateConvergence, baldi1989NormalApproximations, rinott1996MultivariateCLT and further extensions in chen2004NormalApproximation. Recent extensions allow for observations located on a stochastic network, for example kojevnikov2020LimitTheorems considers estimation and inference under network dependence where they assume there is a random network and the economic distance is determined by the length of the shortest path connecting two nodes. vainora2020NetworkDependence defines network stationarity and studies estimation and inference.
We adopt the approach to model the cross-sectional dependence by assuming that there is a measure of economic distance \(\mathbf{d}\) between the observations and generalise current results on sums of cross-sectionally dependent samples to U-statistics. We will consider weak cross-sectional dependence in the sense that two groups of observations will be approximately independent if their economic distance is large, which is closely related to the settings in jenish2009CentralLimit which assumes the observations form random fields, and kojevnikov2020LimitTheorems that specifies the distance to be the distance on a network. For our results to hold, the observation of the economic distance is not necessary, provided that there exists such a distance for which cross-sectional dependence is weak and the mixing rates apply. In practice, the observation of a network can be useful to check whether the conditions hold and for feasible inference (see methods proposed in kojevnikov2020LimitTheorems and kojevnikov2021BootstrapNetwork).
We handle the cross-sectional dependence with Stein's method. The Stein's method begins with the seminal paper stein1972BoundError, which estimates the error of normal approximation for a sum of random variables with a certain dependence structure. It has since become one of the state-of-the-art ways to estimate the distance between two random variables under various probability metrics. Several instances of the Stein's method have been developed under the general idea of “auxiliary randomization” introduced by stein1986ApproximateComputation, such as the exchangeable pairs approach, zero-bias and size-bias couplings, see ross2011FundamentalsSteins. Using the approach of exchangeable pairs, rinott1997CouplingConstructions and dobler2016QuantitativeJong extend and provide convergence rates for the central limit theorems for U-statistics in hall1984CentralLimit and dejong1987CentralLimit respectively, with an independent sample \(\left\{X_{i}: i \in \mathcal{I}\right\}\). However, their construction of exchangeable pairs depends on the independence condition and is not straightforward to be extended to the case of cross-sectional dependence.
We base our analysis on a variant of Stein's coupling method and the decomposable method barbour1989CentralLimit, chen2010SteinCouplings, as stated in (ref). This approach is tailored to handle U-statistics, particularly degenerate ones, and differs from the other Stein's methods typically used for analyzing sums of dependent random variables. We believe this method is of interest in its own right.
Our results are of interest to other strands of literature related to U-statistics, and we name a few examples here. athey2019GeneralizedRandom shows that a generalised random forest has a U-statistic representation. The estimation and inference for multiway clustered data and exchangeable arrays depend crucially on U-statistics; see, for example, chiang2021InferenceHighdimensional, cattaneo2022UniformInference, chiang2023UsingTwoWay, and chiang2024StandardErrors. In some cases, the objective function in an M-estimation takes the form of a U-statistic; see honore1997PairwiseDifference. In panel data settings, this also opens up a new avenue for showing limiting distributions along the cross-sectional dimension. anatolyev2020LimitTheorems proposes a method that allows for unspecified cross-sectional dependence and utilizes independence or martingale properties along the time series dimension. zaffaroni2019FactorModels considers conditional factor model with short time series. Our results can be used in the settings when the number of time series observations is limited, due to data availability or stationarity considerations.
This paper is organised in the following way. Section (ref) sets up the model and defines the key ingredients required for the theoretical results, with some examples given in Section (ref). Section (ref) presents the main findings where we estimate the error of normal approximation to the non-degenerate and degenerate U-statistics. The convergence rates will depend on the mixing rates, the sparsity of the cross-sectional dependence, and moments of the kernel functions. Section (ref) provides an application of the CLT for degenerate U-statistics on the nonparametric specification tests. Section (ref) concludes. Proofs for the theoretical results and further extensions are collected in the Appendix.
For \(n\in \mathbb{N}\), \([n]\) denotes the set \(\left\{1, 2,\dots, n\right\}\). \(\mathbb{I}(\cdot)\) is the indicator function. For a vector of indices \(\mathbf{i} = (i_{1},\dots, i_{q})\), we will write \(i\in \mathbf{i}\) if \( i\in \left\{i_{1},\dots, i_{q}\right\}\), and we define the subvector \(\mathbf{i}_{-i_{k}} = ({i_{1},\dots, i_{k - 1}, i_{k + 1}, \dots, i_{q}})\). For a collection of random vectors \(\left\{X_{i}:i\in \mathcal{I}\right\}\), \(\sigma\left(X_{i}: i\in \mathcal{I}\right)\) denotes the \(\sigma\)-algebra generated by \(\left\{X_{i}\right\}\), and \(\tilde{X}_{i}\) denotes a random vector that has the same distribution as \(X_{i}\). \(\left\lvert\cdot\right\rvert\) denotes the Euclidean norm and \(\left\lVertX\right\rVert_{p} = \left[\operatorname{\mathrm{E}}\left\lvertX\right\rvert^{p}\right]^{{1}/{p}}\) denotes the \(L^{p}\)-norm of \(X\). We will write \(\operatorname{\mathrm{E}}(\cdot)\) for unconditional expectations and \(\operatorname{\mathrm{E}}^{X} \left(Y\right)\) for conditional expectation \(\operatorname{\mathrm{E}}\left[Y \mid X\right]\). For sequences of positive numbers, \(a_{n} = O(b_{n})\) or \(a_{n}\lesssim b_{n}\) if \(a_{n} \leq C b_{n}\) for some constant \(C > 0\) and sufficiently large \(n\); \(a_{n} =\Theta\left(b_{n}\right)\) if \(a_{n} = O(b_{n})\) and \(b_{n} = O(a_{n})\). A sequence of random variables \(X_{n} \rightsquigarrow X\) if \(X_{n}\) converges in distribution to a random variable \(X\).
Suppose we observe a sample \(\left\{X_{i}: i\in \mathcal{I}_{n}\right\}\), where each \(X_{i} :\left(\Omega, \mathcal{F}, \operatorname{P}\right) \to (\mathbb{R}^{d}, \mathcal{B}_{\mathbb{R}^{d}})\) is a random vector, \(\mathcal{I}_{n}\) is the index set for the entities of interest, and \(\mathcal{B}_{\mathbb{R}^{d}}\) is the Borel \(\sigma\)-algebra on \(\mathbb{R}^{d}\). We will also refer to an index \(i\in \mathcal{I}_{n}\) as a node, to subsets \(\mathcal{I}_{n}'\subset \mathcal{I}_{n}\) as groups of nodes, and we assume the sample size \(\left\lvert\mathcal{I}_{n}\right\rvert = n\). To model cross-sectional dependence, we assume that \(\left(\mathcal{I}_{n}, \mathbf{d}\right)\) is a metric space with distance \(\mathbf{d}: \mathcal{I}_{n}\times \mathcal{I}_{n} \to \left[0, \infty\right]\). For example, \(\mathbf{d}\) can measure physical distance, as in a spatial model, or the strength of economic links between entities. For two groups of nodes \(\mathcal{I}_{1}, \mathcal{I}_{2} \subset \mathcal{I}_{n}\), we naturally define the distance
The strength of cross-sectional dependence among the observations \(\left\{X_{i}:i\in \mathcal{I}_{n}\right\}\) will be measured by the \(\beta\)-mixing coefficients. The \(\beta\)-mixing coefficient between two \(\sigma\)-algebras \(\mathcal{F}_{1}, \mathcal{F}_{2}\) is defined as
where the supremum is taken over the set of all finite partitions \(\mathcal{A} = \left\{A_i: 1 \leq i \leq n_A\right\}\) and \(\mathcal{B}=\left\{B_{j}: 1\leq j \leq n_B\right\}\) of \(\Omega\) for some \(n_A, n_B < \infty\), with \( \mathcal{A}\subset \mathcal{F}_{1}\) and \(\mathcal{B}\subset \mathcal{F}_{2}\). It is well known that \(\beta\left(\mathcal{F}_{1}, \mathcal{F}_{2}\right) = 0\) if and only if \(\mathcal{F}_{1}\) and \(\mathcal{F}_{2}\) are independent doukhan1994Mixing.
We will make use of the following sets, consisting of ordered pairs of groups of nodes that are at least distance \(m\) apart and satisfy restrictions on their sizes. This definition is common in the random-field literature jenish2009CentralLimit.
The restrictions we impose on cross-sectional dependence are stated in the following assumptions.
A few comments are in order. The assumptions on weak dependence in the observations \(\left\{X_{i}\right\}\) are similar to standard assumptions in the random-field literature; see, for example, jenish2009CentralLimit. The first lower-bound condition is imposed to rule out increasing concentration of observations into arbitrarily small neighborhoods as the sample grows. What is essential for our counting arguments is that local multiplicity remains controlled. More generally, the distance here should be interpreted as a device for controlling the strength of dependence. In spatial settings, even observations at the same physical location need not be perfectly dependent, so it is reasonable to allow for a strictly positive “dependence distance” in such cases. The second assumption extends familiar mixing conditions from time series and random fields to the more general settings we consider.
We choose to state our conditions in terms of \(\beta\)-mixing, or absolute regularity, other than some other mixing concepts considered in the literature on limit theorems for sums of random variables (for example, kojevnikov2020LimitTheorems considers \(\psi\)-mixing) because \(\beta\)-mixing random variables have stronger coupling properties, as reflected in Berbee's Lemma ((ref)), which is essential for handling the nonlinear U-statistics kernel. Berbee's Lemma can be extended to the case of \(\alpha\)-mixing under a Lipschitz continuity condition on the kernel, see bradley1983ApproximationTheorems and dehling2010CentralLimit. Therefore, our results can be readily generalised to allow for \(\alpha\)-mixing under Lipschitz continuity. Since the extension to \(\alpha\)-mixing complicates the notation without adding much insight, we maintain the assumption of \(\beta\)-mixing. Our results are first stated for the case where the distance \(\mathbf{d}\) is deterministic, and in Section (ref) we discuss extensions that allow the underlying sample \(\left\{X_{i}:i\in \mathcal{I}\right\}\) to be conditionally \(\beta\)-mixing by providing a conditional version of Berbee's lemma. This extension is relevant in practice.
It is well known from the literature on random fields and network dependence that the topology of the index space \(\mathcal{I}_{n}\) plays an important role in determining the strength of cross-sectional dependence. Indeed, from the definition of \(\mathcal{P}_{n}\left(n_{1},n_{2},m\right)\) in (ref) and (ref) in (ref), we see that cross-sectional dependence is stronger if the mixing coefficients \(\beta(\sigma(\mathcal{I}_{1}), \sigma(\mathcal{I}_{2}))\) are larger for each \((\mathcal{I}_{1}, \mathcal{I}_{2})\), or if more nodes in \(\mathcal{I}_{n}\) lie within distance \(m\) of one another and the indices are more clustered. We therefore need further requirements on the topology of the index space \(\mathcal{I}_{n}\) with respect to the distance \(\mathbf{d}\). In particular, we define the following quantities that measure the “sparsity” of the index space.
The normal approximation bounds below also use counts of index vectors that fall into different dependence configurations. For a vector of \(q\) indices \(\mathbf{i}=(i_1,\ldots,i_q)\) and \(m<\infty\), form a graph on the indices \(\{i_{1},\ldots,i_{q}\}\) by connecting positions \(i_{r}\) and \(i_{s}\) whenever \(d(i_r,i_s)\leq m\). We say that \(\mathbf{i}\) is \(m\)-connected if this graph is connected, and denote by \(\tau_q^m\) the number of vectors of \(q\) indices that are \(m\)-connected. More generally, \(\tau_{q_1,\ldots,q_S}^m\) denotes the number of vectors of \(q\) indices whose \(m\)-connected components have sizes \(q_1\geq \cdots \geq q_S\), where \(\sum_{s=1}^S q_s=q\). For example, \(\tau_{2,2}^m\) counts vectors \((i_1,i_2,i_3,i_4)\) whose four positions can be partitioned into two \(m\)-connected components that are separated from each other by more than \(m\). And \(\tau_{1,1,1,1}^m\) counts vectors \((i_1,i_2,i_3,i_{4})\) where any two of the four indices have distance greater than \(m\). The full definition based on equivalent classes is given in the Appendix.
The distribution of \(\tau^{m}_{q_{1},\dots,q_{S}}/n^{q}\) over different profiles \((q_{1},\dots,q_{S})\) provides important information about the sparsity of the index space \(\mathcal{I}_{n}\). Consider the case \(q = 2\) with fixed \(m>0\). Then \(\tau_{2}^{m}\) is the number of ordered pairs \((i,j)\in \mathcal{I}_{n}^{2}\) that lie in one another's \(m\)-neighbourhoods, \(i\stackrel{m}{\leftrightarrow} j\), while \(\tau_{1,1}^{m}\) is the number of pairs \((i,j)\) such that \(\mathbf{d}(i,j) > m\). When the index space is sparse, the proportion of pairs of nodes that are more than \(m\) units apart, \({\tau_{1,1}^{m}}/{n^{2}}\), is expected to be higher. When the index space is less sparse, the relative proportion \(\tau_{2}^{m} / n^{2}\) increases.
It is easy to compute the profile \(\pi_{m}(\mathbf{i})\) of a given vector of indices \(\mathbf{i}\) and each \(\tau^{m}_{q_{1},\dots,q_{S}}\) when the distance is observed. Furthermore, a simple bound on the sparsity measure \(\tau^{m}_{q_{1},\dots, q_{S}}\) can be obtained if we have a bound on the growth rate of the maximum neighbourhood sizes. Let \(\eta_{m,i} = |\mathcal{N}_{i}^{m}|\) be the size of the \(m\)-neighbourhood of \(i\) for each \(i \in \mathcal{I}_{n}\), and let \(\eta_{m} = \max_{i} \eta_{m,i}\), which measures the maximum size of \(m\)-neighbourhoods. Then we have the following lemma, which is useful in many applications.
Our main results provide bounds on the Wasserstein distance between a standard normal random variable and a suitably normalised U-statistic \(S_{n}\) constructed from a weakly cross-sectionally dependent collection of observations \(\left\{X_{i}:i\in \mathcal{I}_{n}\right\}\). We provide conditions under which the Wasserstein distance converges to \(0\) as the sample size increases. Since convergence in Wasserstein distance implies convergence in distribution, the central limit theorems follow from these general results. The Wasserstein distance between two random variables \(V_{1}\) and \(V_{2}\) is defined by
where \(\mathbb{L}_{1} = \{g: \mathbb{R} \to \mathbb{R} : |g(x) - g(y)| \leq |x - y|\}\) is the class of \(1\)-Lipschitz continuous functions.
We discuss some examples of cross-sectional dependence structures that can be incorporated into our framework and how their sparsity can be measured.
In addition to assumptions on the mixing rate and the sparsity of the index space, we also require moment restrictions on the kernel function $H_{n}$. For simplicity, we omit the subscript $n$ when the dependence of $H_{n}$ on the sample size is clear. Let \(\tilde{X}_{j}\) be an independent copy of \(X_{j}\). We assume the following regularity conditions on the kernel function.
The requirement of a symmetric kernel can be relaxed by considering \(\frac{1}{2}(H(x,y) + H(y,x))\). The square-integrability assumption ensures that the Hoeffding decomposition is well defined in the non-degenerate case and is important for the CLT for degenerate U-statistics, as in hall1984CentralLimit. We define the following quantities, analogous to Hall's condition (ref), with $H_{ij} := H_{n}(X_{i}, X_{j})$ and \(H_{i \tilde{j}} = H(X_{i}, \tilde{X}_{j})\), for each \(i\neq j\),
and let \(\tilde{X}_{k_{2}}\) be independent of \(X_{k_{1}}\) and identically distributed as \(X_{k_{2}}\),
Our bound on the normal approximation error will be a consequence of the interplay between sparsity, mixing rates, and moments.
For non-degenerate U-statistics with cross-sectionally dependent underlying processes, we use the following Hoeffding decomposition. Let \(\hat{H}_{k}(x) = \operatorname{\mathrm{E}} H(X_{k}, x)\) and \(\theta_{ik} = \operatorname{\mathrm{E}} \hat{H}_{k}(X_{i}) = \operatorname{\mathrm{E}} \hat{H}_{i}(X_{k})\), for \(x\) in the support of \(\left\{X_{i} : i\in \mathcal{I}\right\}\). Then
where, for \(h_{i} = (n-1)^{-1} \sum_{k\neq i} (\hat{H}_{k}(X_{i}) - \theta_{ik})\),
It is now straightforward to obtain the following Law of Large Numbers for non-degenerate U-statistics. It generalises existing results for independent observations; see serfling1980ApproximationTheorems and Lemma 3.1 in powell1989SemiparametricEstimation.
Stein's method can then be applied to obtain a central limit theorem by showing that the projection part converges in distribution to a Gaussian random variable and that the remainder term is asymptotically negligible. By constructing a specific Stein's coupling in (ref) and using Berbee's Lemma ((ref)), we obtain the following normal approximation result.
Similar to the classical central limit theorems for time series and random fields, there is a trade-off between the moment restrictions and the mixing rates. When \(\beta(n_{1},n_{2},m)\) tends to \(0\) faster as \(m\to \infty\), we can take a smaller \(\delta\) and thus relax the moment assumptions. There is also a trade-off, through the choice of \(m\), between the sparsity of the cross-sectional dependence and the mixing rates. When we choose a larger \(m\), the mixing rate \(\beta(n_{1},n_{2},m)\) gets smaller at the cost of a larger \(\tau^{m}_{q_{1},\dots,q_{s}}\). In the independent case, one may view distinct nodes as having dependence distance \(\mathrm d(i,j)=\infty\). As a result, for any fixed finite \(m\), the mixing terms vanish, \(\eta_m=1\), and (ref) yields the usual \(n^{-1/2}\) normal-approximation rate.
As a corollary, if we assume bounded third moments and bounds on the maximum neighbourhood sizes, then the following central limit theorem follows by combining (ref) with (ref) and taking \(\delta = 1\).
We make a few remarks about the theorems. Firstly, (ref) provides a Berry--Esseen type result under \(\beta\)-mixing; see Theorem 2.7 of dobler2015NewBerryEsseen for independent observations. Because of (ref) and the decomposable approach we use to construct the coupling variables, only bounds on the third moment are necessary, in contrast to kojevnikov2020LimitTheorems, which requires moments of order higher than \(4\) and \(\psi\)-dependence. Secondly, we have used only \(\eta_{m}\), the maximum \(m\)-neighbourhood size, to describe the sparsity. In the literature on random fields and network dependence, summability conditions involving the sizes of “neighbourhood shells” \(\mathcal{S}_{i}^{m} = \{j: m < \mathbf{d}(i,j) \leq m+1\}\), multiplied by the mixing rates, are often assumed bolthausen1982CentralLimit,kojevnikov2020LimitTheorems. These conditions usually provide a more precise description of the dependence structure, leading to a weaker set of conditions on the mixing rates, since \(|\mathcal{S}_{i}^{m}| \leq |{\mathcal{N}_{i}^{m}}|\). In any case, we have the freedom to choose a suitable sequence \(m_{n}\), so we do not lose much applicability by stating (ref) and (ref) in terms of \(\eta_{m}\) with the benefits of simpler notations.
Next, we consider a HAC estimator of the “long-run” variance \(\nu_{n}^{2}\). For non-degenerate U-statistics, thanks to the Hoeffding decomposition, this problem is similar to estimating the variance of sample means under cross-sectional dependence. For example, kojevnikov2020LimitTheorems and kojevnikov2021BootstrapNetwork propose a network HAC estimator for the variance of the sample mean and a bootstrap procedure under network dependence, provided that knowledge of the underlying network is available.
Consider a tapering function \(\kappa: \mathbb R_{+} \to [0,1]\) that is nonincreasing, with \(\kappa(0)=1\) and \(\kappa(z)=0\) for \(z>1\). For some positive sequence \(b_n\), let \(\kappa_{ij}(b_n)=\kappa\!\left(\frac{\mathbf d(i,j)}{b_n}\right)\), \(Q_i := \frac{1}{n-1}\sum_{k\neq i} H_{ik}\), and \(\bar Q := \frac{1}{n}\sum_{i\in I_n} Q_i\). We consider the following HAC estimator: \[ \hat \nu_n^2 := \sum_{i\in I_n}\sum_{j\in I_n} \kappa_{ij}(b_n) (Q_i-\bar Q)(Q_j-\bar Q). \]
To show consistency of the estimator, we impose the following additional conditions.
These conditions are similar to those in kojevnikov2020LimitTheorems for the network HAC estimator of the variance of the sample mean. The first two conditions are a standard moment restriction on the kernel function \(H\) and a bandwidth condition requiring the underlying dependence structure to be sufficiently sparse and weak. The third condition is the analogue of Assumption 4.1(ii) in kojevnikov2020LimitTheorems, which bounds the contribution to the bias of the HAC estimator from each pair of nodes in a neighbourhood shell \(\{(i,j): m < \mathbf{d}(i,j) \leq m+1\}\). The last condition is a homogeneity condition on the mean, so that we can eliminate the bias arising from variation in the mean across \(i\in I_n\); this condition holds, for example, when the underlying process is stationary.
Hence, under suitable regularity conditions, we can use the HAC estimator when the underlying distance \(\mathbf{d}\) is observed.
For the degenerate case, \(S_{n} = \sum_i \sum_{k\neq i} H_{ik}\) and \(\hat{H}_{k}(x) = \operatorname{\mathrm{E}} H(X_{k},x) = 0\). The following theorem follows from the more general result (ref) in (ref) and (ref), which generalises hall1984CentralLimit's result to the case of cross-sectional dependence for degenerate U-statistics.
There is also a trade-off between the mixing and moment conditions. If the mixing coefficient converges to zero more rapidly, one can choose a smaller value of \(\delta\), thereby weakening the required moment assumptions. The choice of \(m\) reflects a further trade-off: a larger \(m\) improves the mixing coefficient \(\beta(m)\), but at the cost of increasing the neighborhood size \(\eta_m\). In the special case of an independent sample, \(\beta(m)=0\) for all \(m\geq 1\). Taking \(m=1\) gives \(N_i^m=\{i\}\) and hence \(\eta_m=1\). The bound in Theorem 3.5 then recovers the classical independent-sample benchmark conditions for degenerate U-statistics, as in hall1984CentralLimit.
As a corollary, the following CLT holds under a geometric mixing rate and slightly stronger moment conditions than those in hall1984CentralLimit. In practice, once we have a specific choice of the kernel function \(H\), the conditions depend on the mixing rate and moment restrictions, which can be checked empirically, as in classical time-series settings. The sparsity restriction can be checked if the distance function is known; otherwise, it may be justified through economic reasoning.
(ref) and the corollary appear to be suitable for index spaces that are relatively “uniform”, in the sense that the \(\eta_{m,i}\) are of the same order across \(i\in \mathcal{I}\), so that using the maximum cardinalities of neighbourhoods \(\eta_{m}\) in (ref) is a reasonable choice. For more asymmetric networks, it may be more appropriate to consider the general result and use a more precise counting method in the estimation, as in (ref).
In practice, we also need to estimate the long-run variance $s_{n}^{2}$ for degenerate U-statistics. In the time-series setting, fan1999CentralLimit show that a variance estimator constructed as though the sample were independent can nonetheless be consistent, and demonstrate this in the context of nonparametric specification testing. More generally, bootstrap methods for U-statistics have been extensively studied under independence and time-series dependence. For independent samples, bickel1981AsymptoticTheory develop bootstrap theory for U-statistics, while arcones1992BootstrapStatisticsa extend these results to degenerate U-statistics. Under weak dependence, dehling2010CentralLimit establish bootstrap validity for U-statistics of mixing processes, and leucht2013DependentWild propose a dependent wild bootstrap for degenerate U-statistics.
Extending bootstrap methods for U-statistics to settings with cross-sectional dependence is an important direction that we intend to explore in future work, but it seems beyond the scope of this paper. Instead, we provide a complementary result to the central limit theorem for degenerate U-statistics ((ref)) by simplifying the variance expression. This result is used in the proof of (ref) to show that a variance estimator $\hat{s}_n^2$, constructed as if the sample were independent, remains consistent for the variance of the test statistic under suitable moment and sparsity conditions, in the same spirit as fan1999CentralLimit.
In statistics and econometrics, there has always been a trade-off between model complexity and statistical inference. For a statistical model \(\left\{\operatorname{P}_{\gamma}: \gamma\in \Gamma\right\}\), where \(\operatorname{P}_{\gamma}\) denotes the probability law that governs the data generating process indexed by a parameter \(\gamma\) in the parameter space \(\Gamma\), we face the risk of model misspecification if we estimate a restricted model \({\Gamma}_0\), while the true parameter lies outside the models we consider, \(\gamma_{0} \notin {\Gamma}_0\). On the other hand, with a more complex model, we may lose identification power and efficiency. This trade-off necessitates close examination of how we specify our models and calls for the development of specification tests. We will focus on the test for the specification of a regression model, \(Y_{i} = g(Z_{i}) + u_{i}\), where we observe \(\left\{(Z_{i}, Y_{i}): i\in \mathcal{I}\right\}\) allowing for cross-sectional dependence among the observations.
A vast literature exists on such specification tests. Bearing in mind the trade-off mentioned, this paper will consider the specification test where we have the null hypothesis that a regression function lies in a parametric family indexed by some finite-dimensional parameter against a general nonparametric alternative, i.e. \(H_{0}: g(z) = g(z,\gamma)\) for a known parametric function indexed by a parameter \(\gamma \in \Gamma\subset \mathbb{R}^{d}\) against the alternative that \(g(z)\) is a smooth nonparametric function. Several such tests have been proposed based on the kernel smoothing method; for example, hardle1993ComparingNonparametric compares the \(L_{2}\) distance between the parametric and nonparametric fit, while zheng1996ConsistentTest, fan1999CentralLimit, and li2007NonparametricEconometrics propose tests based on the idea that if the parametric model is correctly specified, then a kernel smoothing fit for the residuals should be approximately zero. Such tests avoid the random denominator problem. Additionally, fan1999CentralLimit dealt with absolutely regular time series data, which was later generalised to broader mixing concepts by gao2008CentralLimit and others. Finally, horowitz2001AdaptiveRateOptimal and horowitz2002AdaptiveRateOptimal propose adaptive and rate optimal tests that are uniformly consistent against alternatives by considering multiple bandwidths.
For these test statistics, the dominating parts are second-order U-statistics and the Central Limit Theorems for the test statistics are derived under conditions of hall1984CentralLimit, dejong1987CentralLimit, and their extensions to the time series settings. With the normal approximation theory developed in the previous section, we are now able to generalise the nonparametric specification test to allow for cross-sectional dependence.
Suppose we have a collection of cross-sectional observations \(\left\{(Y_{i}, Z_{i}) : i\in \mathcal{I}\right\}\), where \(\mathcal{I}\) is the index set and \(n =\left\lvert\mathcal{I}\right\rvert\) is the sample size, \(Y_{i}\) is a random scalar, and \(Z_{i}\in \mathbb{R}^{d}\). Suppose for \(i \in \mathcal{I}\), \(\operatorname{\mathrm{E}}(Y_{i}\mid Z_{i}) = g(Z_{i})\). We wish to test the null hypothesis of a parametric model specification for \(g(\cdot, \gamma)\) where \(\gamma\in \Gamma_0\) against a nonparametric alternative,
Let \(u_{i} = Y_{i} - g(Z_i, \gamma_0)\), we will allow for cross-sectional dependence in \(X_{i} = \left(u_{i}, Z_{i}\right)\) under the null hypothesis \(H_{0}\). Let \(\hat{\gamma}\) be the nonlinear regression estimate for \(\gamma\) which is assumed to be \(\sqrt{n}\)-consistent, and \(\hat{u}_{i} = Y_{i} - g\left(Z_{i}, \hat{\gamma}\right)\). Following fan1996ConsistentModel and li2007NonparametricEconometrics, we base our test on the kernel estimator for \(\operatorname{\mathrm{E}}\left[u_{i}\operatorname{\mathrm{E}}\left[u_{i}\mid Z_{i}\right] f_{i}(Z_{i})\right]_{i})}\), which equals zero under the null hypothesis \(H_{0}\), where \(f_{i}(\cdot)\) is the marginal density function of \(Z_{i}\). Consider the following statistic,
where \(b\) is a bandwidth we choose. We will write \(K_{ij} = K\left(\frac{Z_{i} - Z_{j}}{b}\right)\) for a kernel weighting function \(K\). For later use, let $I_{n}^{\circ} = \sum_{i}\sum_{j\neq i} b^{- \frac{d}{2}} u_{i}u_{j} K\left(\frac{Z_{i} - Z_{j}}{b}\right)$. This oracle statistic \(I_{n}^{\circ}\) is a U-statistic of order \(2\) with samples \(\left\{\left(u_i, Z_i\right) : i \in \mathcal{I}\right\}\) and kernel
where \(x = (x_{1}, x_{2}')'\) and \(y = (y_{1}, y_{2}')'\in \mathbb{R}^{1}\times \mathbb{R}^{d}\). The feasible statistic \(I_{n}\) in (ref) is the plug-in version obtained by replacing \(u_i\) with \(\hat u_i\). Let \(\hat{s}_{n}^{2} = 2 b^{- d} \sum_{i}\sum_{j\neq i} \hat{u}^{2}_{i}\hat{u}^{2}_{j} K_{ij}^{2}\), we will consider the following test statistic
We make the following assumptions on the data generating process. Firstly, we allow for cross-sectional dependence in the sample but require the sample to be \(\beta\)-mixing.
And we make the following assumptions on the moments of the error terms, and the smoothness of the density function of \(Z_{i}\)'s and the regression function.
This assumption that the supports of \(Z_{i}\)'s are not disjoint will be satisfied if we assume stationarity. The following assumptions on the kernel smoothing function \(K\) are standard in the literature and satisfied for commonly chosen kernel smoothing functions, such as normal kernel, Epanechnikov kernel, and uniform kernel.
Our main result is the following distribution theory for the test statistics.
In this section, we conduct simulations to show the finite-sample performance of the proposed specification test with cross-sectionally dependent data.\footnote{The codes for implementing the nonparametric specification test and the simulations are available at \href{https://github.com/lwg342/specification-test-csc-code}{https://github.com/lwg342/specification-test-csc-code}.} The simulation specification is close to the one in horowitz2001AdaptiveRateOptimal. For each simulation we generate \(n = 2000\) samples of \((z_{i}, y_{i})\), where \(z_{i} \sim N\left(0,25\right)\) and the null-hypothesis model is
where we take the true \(\gamma_{0} = \gamma_{1} = 1\), and the alternative model is
where \(\phi\) is the standard normal density function and \(\tau\) is the bandwidth parameter where \(\tau = 1\) or \(0.25\), it controls the shape of deviation from the null model with a smaller \(\tau\) implying a more “spiked” deviation. The \(\psi\) parameter controls the overall magnitude of deviation of the true model from the null models. The error terms \(u_{i}\) are generated according to three types of models. The first is the benchmark i.i.d. model, where \(u_{i} \sim N\left(0,1\right)\). The second is an AR(1) model with \(u_{i} = \rho u_{i - 1} + \upsilon_{i}\) where \(\upsilon_{i} \sim N\left(0,1\right)\). The third is a two-way clustering structure. Assume there is a bijective mapping \(\pi: [n] \to [n_{1}]\times [n_{2}]\), such that \(i,j\) are independent if \(\pi(i) = \left(\pi_{1}(i), \pi_{2}(i)\right)\) and \(\pi(j)\) share no common coordinates. In particular the model we consider is \(u_{i} = \frac{1}{\sqrt{2}} \left[\lambda_{\pi_{1}(i)} + F_{\pi_{2}(i)}\right]\) where each \(\lambda, F\) are independent standard normal random variables. For the two-way clustering design, we use \((n_{1}, n_{2}) = (40, 50)\) and \((100, 20)\), so that the total sample size remains \(n = 2000\) in each case.
For each simulated dataset \(\left(z_{i}, y_{i}\right)\), we compute the residual \(\hat{u}_{i}\) from the linear regression of \(y_{i}\) on \(z_{i}\). We then compute the test statistic \(\mathbf{T}_{n}\) using the residuals \(\hat{u}_{i}\), with fixed bandwidth \(b_{n} = n^{-1/4}\), and reject the null when \(\mathbf{T}_{n} > 1.645\). We repeat the simulation 2000 times and compute the empirical size and power of the test. The results are shown in Table (ref).
It can be seen from the table that when the null hypothesis \(\mathbb{H}_{0}\) is true, the empirical probability of the test statistic rejecting the null hypothesis is close to the nominal level of 0.05, even when we allow for cross-sectional dependence. When the alternative hypothesis \(\mathbb{H}_{1}\) is true, the power of the test statistic depends on the deviation from the null as well as the cross-sectional dependence structure. In general when the deviation from the null is large, then the power of the test is reasonably good. When the deviation gets closer to being local, the effect of cross-sectional dependence is more noticeable. With stronger cross-sectional dependence, the power generally decreases, as we use the consistent variance estimator in \autoref*{thm:degenerate u variance} which requires weak dependence and sparsity. In this case, perhaps we could consider more sophisticated variance estimators that take into account the cross-sectional dependence structure.
We apply the nonparametric specification test to a dataset on Airbnb listings in Dublin to test whether the log of listing prices is linear in the log-distance to Dublin city center, conditional on a set of listing and host characteristics. Airbnb listings are likely to exhibit local spatial dependence because nearby listings may share unobserved neighbourhood quality, common demand shocks, and local amenities. The data are obtained from Inside Airbnb\footnote{The data can be downloaded from \url{https://insideairbnb.com}} and are based on a Dublin extract that contains only the quarterly listings snapshot scraped on September 16, 2025. We restrict attention to listings classified as “Entire home or apartment.” After retaining observations with positive prices, the sample contains 2,958 listings. Excluding observations with missing regression controls yields a final testing sample of 2,494 listings.
The outcome variable is \(Y_i=\log(\text{price}_i)\), and the main regressor of interest is the log-distance to the city center, as in the literature combes2019CostsAgglomeration,gupta2022FlatteningCurve. Specifically, we define \(Z_i=\log(d_i+0.25)\), where \(d_i\) denotes the straight-line distance, measured in kilometers, from listing \(i\) to a central Dublin reference point at \((53.35^\circ\text{N},\; 6.26^\circ\text{W})\), which is the location of the Spire of Dublin. We add the small constant $0.25$, following gupta2022FlatteningCurve, to regularize the transformation for listings located very close to the city center. Robustness checks with alternative distance transformations are reported in the Appendix, and the main conclusions are not sensitive to this choice.
The null hypothesis is \(H_0:\ \mathbb{E}[Y_i\mid Z_i,W_i]=\alpha+\beta Z_i+W_i'\delta\), where \(W_i\) is a vector of observable listing and host characteristics. We estimate the null model by OLS and then compute the proposed test statistic with Epanechnikov kernel over \(Z_i\). As a benchmark choice, we use the rule-of-thumb bandwidth \(b_0=1.06\,\hat{\sigma}(Z)\,n^{-1/5}\), which equals \(b_0=0.212\) in our sample.
Table (ref) reports the OLS estimation results (Panel A) and the specification test statistics (Panel B). In Panel B, we report the benchmark test statistic computed as in (ref) together with its asymptotic $p$-value under the simple variance estimator. Because Airbnb listings are likely to be spatially dependent, we also report results from dependence-robust bootstrap procedures. We implemented a neighbourhood block wild bootstrap, grid block wild bootstraps, and Leucht-style spatial dependent multiplier bootstraps evaluated at several dependence scales. The neighbourhood and grid procedures are in the spirit of wild and block bootstrap methods for heteroskedastic and dependent data mammen1993BootstrapWild,zhu2007BootstrappingEmpirical,cameron2011RobustInference. The multiplier procedure is instead motivated by leucht2013DependentWild, although our spatial implementation is an adaptation to cross-sectional proxy distance rather than a literal application of their theorem. The residual dependence diagnostics reported in the Appendix indicate that dependence is local. In the Dublin dataset, the mean residual cross-product within 1 km is about 0.008, while beyond 2 km it is essentially zero. This pattern informs our choice of the hyperparameters for the bootstrap procedures.
First, we compute a neighbourhood block wild bootstrap. Let \(g(i)\) denote the neighbourhood block of listing \(i\), based on the variable neighbourhood. In each bootstrap replication, we draw a Rademacher multiplier \(\xi_{g(i)}\in\{-1,+1\}\), and the bootstrap outcome is \(Y_i^*=\hat m_0(X_i)+\xi_{g(i)}\hat u_i\), where \(\hat m_0(X_i)\) is the fitted value under the null. We then refit the same null model using \(Y_i^*\), obtain bootstrap residuals \(\hat u_i^*\), and recompute the studentized bootstrap statistic. The bootstrap \(p\)-value is the empirical upper-tail probability of \(T_n^*(b)\) relative to the observed \(T_n(b)\). In the Dublin file, however, the variable neighbourhood yields only four non-missing blocks, so this procedure should be interpreted only as a coarse diagnostic.
Second, we compute finer grid block wild bootstraps. The bootstrap construction is identical, namely \(Y_i^*=\hat m_0(X_i)+\xi_{g(i)}\hat u_i\), but the block indicator \(g(i)\) is now defined by projected square grid cells. This gives a much finer partition of space: the 0.5 km grid yields 713 blocks and the 1.0 km grid yields 359 blocks.
Third, we report a Leucht-style spatial dependent multiplier bootstrap. Unlike the two block bootstraps, this procedure does not resample the outcomes. Instead, it applies listing-level multipliers \(\xi_i\) whose dependence decays with projected geographic distance, with the decay controlled by the scale parameter \(\ell\). The bootstrap statistic is computed analogously to the original test statistic, but with \(I_n^{*}(b) = \sum_i \sum_{k\neq i}b^{-1/2}\hat u_i\hat u_k K\!\left(\frac{Z_i-Z_k}{b}\right)\xi_i\xi_k\), where the multiplier field is generated as follows. Let \(d(i,j)\) denote projected geographic distance, let \(\eta_j \overset{iid}{\sim} N(0,1)\), and define \[ \xi_i=\sum_{j\in \mathcal N_i(\ell)} w_{ij}(\ell)\eta_j, \qquad w_{ij}(\ell)= \frac{\rho\!\left(d(i,j)/\ell\right)} {\left(\sum_{m\in \mathcal N_i(\ell)}\rho\!\left(d(i,m)/\ell\right)^2\right)^{1/2}}, \] where \(\mathcal N_i(\ell)=\{j:\rho(d(i,j)/\ell)>0\}\). In the reported calculations, we use the Bartlett decay function \(\rho(u)=\max\{1-u,0\}\), so only listings within distance \(\ell\) receive positive weight. This creates a spatially structured perturbation field. We report the results for \(\ell\in\{0.1,0.25,0.5,1.0\}\) km in (ref).
Table (ref) reports the OLS estimates under the null hypothesis together with the corresponding specification test results. In the Dublin dataset, both the asymptotic test and all dependence-robust bootstrap procedures yield \(p\)-values well above conventional significance levels. For Dublin listings, these findings suggest that a linear model provides an adequate approximation to the relationship between price and location in logarithms. The non-rejection is consistent across all dependence-robust procedures. Additional robustness checks, reported in the Appendix, corroborate the same qualitative conclusion.
This paper develops normal approximation results for second-order U-statistics when the underlying observations exhibit cross-sectional dependence. Using Stein's coupling method together with a decomposable argument, we obtain Wasserstein bounds for both non-degenerate and degenerate U-statistics. The results extend the classical theory of U-statistics beyond independent and time-series settings, and also complement recent limit theorems for sums of cross-sectionally or network-dependent random variables.
The main difficulty is that U-statistics combine two sources of dependence: the dependence generated by the nonlinear kernel and the dependence already present in the underlying sample. Our results show that asymptotic normality can still be obtained under conditions on the mixing rate, the sparsity of the dependence structure, and the moments of the kernel. For non-degenerate U-statistics, the projection component is dominant and feasible inference can be conducted using a HAC-type variance estimator when the underlying distance is observed. For degenerate U-statistics, we provide sufficient conditions under which normal approximation remains valid and discuss variance estimation in the spirit of existing specification-testing procedures.
We also illustrate the usefulness of the theory through a nonparametric specification test that allows for cross-sectional dependence. The empirical application to Airbnb pricing demonstrates how the proposed framework can be used in settings where spatial or network dependence is a natural concern.
Several directions remain for future research. First, it would be useful to develop uniform laws of large numbers for U-statistics under cross-sectional dependence. Such results would be relevant for estimators based on pairwise comparisons, such as the pairwise-difference estimators studied by honore1997PairwiseDifference. Second, bootstrap inference for U-statistics under cross-sectional dependence is an important open direction. While bootstrap methods are well understood in classical settings for non-degenerate U-statistics and sample averages, and network block bootstrap methods have been proposed for sample means under network dependence by kojevnikov2021BootstrapNetwork, the bootstrap theory for degenerate U-statistics under general cross-sectional dependence appears substantially more challenging. We leave this problem for future work.