Empirical Process Results for Exchangeable Arrays

Empirical Process Results for Exchangeable Arrays

Introduction

Taking into account dependence between observations is crucial for making correct inference. For instance, different observations may face common shocks, tending to correlate them positively and thus leading to overly optimistic inference when ignored bertrand2004. Such common shocks may arise if the data are polyadic (e.g., dyadic), namely they involve interactions between several units of a given population. An example is international trade, where each observation corresponds to a pair of countries, one exporting and the other importing. We can then expect that two such pairs may be dependent whenever they share at least one country, because of that country's specificities in terms of international trade. Common shocks may also correspond to aggregate fluctuations that affect all units sharing some characteristics. For instance, wages of two individuals may be correlated either because they live in the same geographical area, or because they work in the same sector. We refer to multiway clustering when there are several dimensions along which units may be correlated.

holland1976local, fafchamps2007formation derived variance formulas for linear regressions with dyadic data, while cameron2011 propose similar formulas for multiway clustering. The Stata command ivreg2 and the R package multiwaycov are now used routinely to report standard errors accounting for multiway clustering. However, theory has lagged behind this practice. Tabor2019 shows the asymptotic validity of inference based on holland1976local's suggestion for dyadic data, but for OLS estimators only. graham2018 and graham2019kernel study respectively parametric regressions and density estimation with dyadic data. Regarding multiway clustering, the only papers we are aware of are the recent works of menzel2017 and mackinnon2017. Again, they focus on linear parameters.\footnote{On the other hand and interestingly, menzel2017 studies inference both with and without asymptotically normality. He also shows that refinements in asymptotic approximations are possible using the wild bootstrap.}

In this paper, we establish uniform laws of large numbers (LLN) and central limit theorems (CLT) for such type of data. Uniform LLNs and CLTs are key in showing consistency and asymptotic normality of nonlinear estimators under weak regularity conditions. As such, they have been studied extensively with i.i.d. but also dependent data. We refer to, e.g., vanderVaartWellner1996 and GineNickl2015 for overviews with i.i.d. data, and dehling2002empirical for the case of time series bertail2017, han2019complex. Noteworthy, we obtain these uniform LLNs and CLTs under the same moment restrictions and conditions on the class of functions as those usually considered with i.i.d. data. Thus, statistical results deducted from the uniform LLNs and CLTs with i.i.d. data directly extend to the exchangeable arrays we consider. As a proof of concept, we consider Z-estimators and smooth functionals of the empirical cumulative distribution function (cdf).

We also study consistency of a direct generalization of the standard bootstrap for i.i.d. data to polyadic data. A related bootstrap scheme for multiway clustering is the so-called pigeonhole bootstrap, suggested by mccullagh2000resampling and studied by owen2007pigeonhole, but for which no uniform result has been established so far. For both, we establish weak convergence of the corresponding process. These results imply the validity of the corresponding bootstrap schemes in a wide range of setting, including Z-estimators and smooth functionals of the empirical cdf.

To prove these results, we first argue that polyadic data correspond to dissociated, jointly exchangeable arrays. Similarly, multiway clustering corresponds to dissociated separately exchangeable arrays. We then rely extensively on the so-called Aldous-Hoover-Kallenberg representation Hoover1979, Aldous1981,kallenberg1989 for such arrays. This representation allows us in particular to prove a symmetrization lemma, which is very useful to derive the uniform LLNs and CLTs. This lemma generalizes a similar result for i.i.d. data, but also for U-processes delapena1999. Note that simple LLNs and CLTs have been already proved, or are direct consequences of known results on dissociated, jointly exchangeable arrays. For LLNs, we refer to Eaglesonweber78 and Lemma 7.35 in kallenberg05. For CLTs, see silverman1976limit. But to our knowledge, no abstract uniform LLNs and CLTs have been proved so far for such arrays.

Finally, we illustrate our results with two applications to international trade. In the first, we test whether international trade remains stable from one year to another, using a Kolmogorov-Smirnov test. Given the dependence structure over pairs of countries and through time, the asymptotic distribution of the test under the null is complicated, making the bootstrap attractive. We show that neglecting the dependence between dyads leads to important overrejection of the null hypothesis. Next, we estimate the so-called gravity equation, a very popular model for explaining trade between countries. Since silva2006log, this equation has often been estimated with Poisson pseudo maximum likelihood, an estimator for which our results apply. Again, much fewer explanatory variables are significant at usual levels when accounting for dependence between pairs of countries than when considering such pairs to be i.i.d. observations silva2006log.

The paper is organized as follows. Section (ref) describes the set-up and gives our main results for jointly exchangeable arrays. In addition to uniform LLNs and CLTs, we prove weak convergence of our bootstrap scheme. We also show results for Z-estimators and smooth functionals of the empirical cdf. Section (ref) considers a few extensions. In particular, we study separately exchangeable arrays. An important difference for such arrays is that the multiple dimensions, corresponding to different sources of clustering, may not grow at the same rate. We show that our results still hold in this case. We also study “degenerate” cases (in the same sense as with U-processes) and consider another bootstrap scheme. The two applications to international trade are developed in Section (ref). The appendix presents three key lemmas. In the supplementary material, we present additional extensions. In particular, we generalize our main results to cases where the number of observations for each $k$-tuple (e.g., the number of matches between two sport players) varies. We also display Monte Carlo simulations and all the proofs of our results.

The set up and main results

Set up

Before formally defining our data generating process, we introduce some notation. For any $A\subset \mathbb R$ and $B\subset \mathbb R^k$ for some $k\geq 2$, we let $A^+=A \cap (0,\infty)$ and $$\overline{B}=\left\{b=(b_1,...b_k)\in B: \; \forall (i,j)\in \{1,...,k\}^2, i\neq j, b_i \neq b_j\right\}.$$ We then let $\mathbb{I}_k=\overline{\mathbb N^{+k}}$ denote the set of $k$-tuples of $\mathbb{N}^+$ without repetition. Similarly, for any $n\in \mathbb N^+$, we let $\mathbb I_{n,k}=\overline{\{1,...,n\}^k}$. For any $\bm{i}=(i_1,...,i_k)$ and $\bm{j}=(j_1,...,j_k)$ in $\mathbb N^k$, we let $\bm{i} \odot \bm{j} = (i_1 j_1,...,i_k j_k)$. With a slight abuse of notation, we also let, for any $\bm{i}=(i_1,...,i_k)\in \mathbb N^k$, $\{\bm{i}\}$ denote the set of distinct elements of $(i_1,...i_k)$. For any $r\in\{1,...,k\}$, we let $$\mathcal{E}_r=\left\{(e_1,...,e_k) \in \{0,1\}^k: \sum_{j=1}^k e_j = r\right\}.$$ Finally, for any $A\subset \mathbb N^+$, we let $\mathfrak{S}(A)$ denote the set of permutations on $A$. For any $\bm{i}=(i_1,...,i_k)\in \mathbb N^{+k}$ and $\pi\in \mathfrak{S}(\mathbb N^+)$, we let $\pi(\bm{i})=(\pi(i_1),...,\pi(i_k))$.

We are interested in polyadic data, that is to say random variables $Y_{\bm{i}}$ (whose support is denoted by $\mathcal{Y}$) indexed by $\bm{i} \in \mathbb I_k$. Dyadic data, which are the most common case, correspond to $k=2$. For instance, when considering trade data, $Y_{i_1, i_2}$ corresponds to export flows from country $i_1$ to country $i_2$. In network data, $Y_{i_1, i_2}$ could be a dummy for whether there is a link from $i_1$ to $i_2$. In directed networks, $Y_{i_1,i_2}\neq Y_{i_2,i_1}$, while $Y_{i_1,i_2}=Y_{i_2,i_1}$ in undirected networks. Similarly, $Y_{i_1, i_2, i_3}$ could capture whether $(i_1,i_2,i_3)$ forms a triad or not wasserman1994social. $Y_{\bm{i}}$ could also correspond to data subject to multiway clustering. Then $i_1$,..., $i_k$ are the indexes corresponding to the different dimensions of clustering, for instance geographical areas and sectors of activity. In such cases, however, adaptations of our set-up are needed, and we postpone this discussion to Section (ref) below.

We assume that the random variables are generated according to a jointly exchangeable and dissociated array, defined formally as follows:

The first part imposes that the labelling conveys no information: the joint distribution of the data remains identical under any possible permutation of the labels. The second part states that the array is dissociated: the variables are independent if they share no unit in common. For instance, $Y_{(i_1,i_2)}$ must be independent of $Y_{(j_1,j_2)}$ if $\{i_1,i_2\}\cap \{j_1,j_2\}=\emptyset$. On the other hand, Assumption (ref) does not impose independence otherwise. This is important in many applications. In the international trade example, $Y_{i_1, i_2}$ and $Y_{i_1,i_3}$ are likely to be dependent because if $i_1$ is open to international trade, it tends to export more than the average to any other country. It may also import more from other countries, meaning that $Y_{i_1, i_2}$ and $Y_{i_3,i_1}$ could also be dependent.

Lemma (ref) below is very helpful to better understand the dependence structure imposed by joint exchangeability and dissociation. It may be seen as an extension of de Finetti's theorem to arrays satisfying such restrictions. It is also key in establishing our asymptotic results below.

This result is due to kallenberg1989 but a weaker version, where the equality only holds in distribution, is known as Aldous-Hoover representation Aldous1981,Hoover1979. Accordingly, we refer to (ref) as the AHK representation hereafter. To illustrate it, let us consider dyadic data ($k=2$). Then, according to Lemma (ref), we have, for every $i_1<i_2$,

Thus, in the example of trade flows, the volume of exports from $i_1$ to $i_2$ depends on factors specific to $i_1$ and $i_2$, such as their own GDP, but also on factors relating both, such as the distance between the two countries. (ref) has been also used by bickel2009nonparametric and bickel2011method to model network formation (in which case $Y_{i_1,i_2}=1$ if there is a link between $i_1$ and $i_2$, 0 otherwise). Note also the link between (ref) and U-statistics: $Y_{i_1,i_2}$ would correspond to such a statistic if $\tau$ did not depend on its third argument.

Under Assumption (ref), the $(Y_{\bm{i}})_{\bm{i}\in \mathbb I_k}$ have a common marginal probability distribution, which we denote by $P$. We are interested in estimating and making inference on features of this distribution, such as its expectation or a quantile, based on observing the first $n$ units only, namely the sample $(Y_{\bm{i}})_{\bm{i} \in \mathbb I_{n,k}}$, with $n\geq k$.

Uniform laws of large numbers and central limit theorems

Let $\mathcal{F}$ denote a class of real-valued functions admitting a first moment with respect to the distribution $P$ and let $Pf$ denote the corresponding moment $\mathbb E\left[f(Y_{\boldsymbol{1}})\right]$ (with $\boldsymbol{1}$ the $k-$tuple $(1,...,k)$). To avoid measurability issues and the use of outer expectations subsequently, we maintain the following assumption:

Assumption (ref) is not necessary but often imposed cherno2014,Kato2019. We refer to Kosorok2006 (Kosorok2006, pp.137-140) for further discussion.

In this section, we study the empirical measure $\mathbb{P}_{n}$ and the empirical process $\mathbb{G}_{n}$ defined on $\mathcal{F}$ by $$\mathbb{P}_{n}f=\frac{(n-k)!}{n!}\sum_{\bm{i} \in \mathbb I_{n,k}}f(Y_{\bm{i}}),$$ $$\mathbb{G}_{n}f=\sqrt{n}\left(\mathbb{P}_{n}f - Pf\right).$$ Let $\ell^\infty(\mathcal{F})$ denote the set of bounded functions on $\mathcal{F}$. We prove below that under restrictions on $\mathcal{F}$, $\mathbb{P}_{n}f$ converges almost surely to $Pf$ uniformly over $f\in \mathcal{F}$, while $\mathbb{G}_{n}$ converges weakly in $\ell^\infty(\mathcal{F})$ to a Gaussian process. We refer to, e.g., vanderVaartWellner1996 for a formal definition of weak convergence of empirical processes. These results, stronger than pointwise convergence of $\mathbb{P}_{n}f$ and $\mathbb{G}_{n}f$, are key in establishing the consistency and asymptotic normality of, e.g., smooth functionals of the empirical cdf or Z- and M-estimators. We consider briefly applications in Section (ref) below, and refer to Part 3 of vanderVaartWellner1996 for a more comprehensive review of statistical applications of empirical process results.

We use the rate $\sqrt{n}$ to normalize $\mathbb P_nf -Pf$, though we have $n!/(n-k)!$ different random variables. In general, we cannot expect a better rate of convergence. To see this, let $(X_i)_{i\in\mathbb N^+}$ be i.i.d. random variables and let $Y_{\bm{i}}=\sum_{j\in \{\bm{i}\}} X_j$. Then $(Y_{\bm{i}})_{\bm{i}\in \mathbb I_k}$ satisfies Assumption (ref), and $\mathbb P_nf$ boils down to an average over $n$ i.i.d. terms only. In some cases, however, for instance if the $(Y_{\bm{i}})_{\bm{i}\in \mathbb I_k}$ are i.i.d., the convergence rate is faster than $\sqrt{n}$.\footnote{ As with U-statistics, we expect different rates depending on the degree of “degeneracy”.} Theorem (ref) below remains valid in such cases, but the limit Gaussian process is then degenerate. We come back in more details to such cases in Section (ref) below.

Let us now introduce the restrictions on $\mathcal{F}$ that we use to obtain uniform laws. We require additional notation for that purpose. For any $\eta>0$ and any seminorm $||\cdot||$ on a space containing $\mathcal{F}$, $N(\eta,\mathcal{F},||\cdot||)$ denotes the minimal number of $||\cdot||$-closed balls of radius $\eta$ with centers in $\mathcal{F}$ needed to cover $\mathcal{F}$. $N_{[\;]}(\eta,\mathcal{F},||\cdot||)$ denotes the minimal number of $\eta$-brackets needed to cover $\mathcal{F}$, where an $\eta$-bracket for $f\in\mathcal{F}$ is a pair of functions $(\ell, u)$ such that $\ell \leq f \leq u$ and $||u-\ell||<\eta$. The seminorms we consider hereafter are $\|f\|_{\mu,r}=(\int |f|^rd\mu)^{1/r}$ for any $r\geq 1$ and probability measure or cdf $\mu$. Hereafter, an envelope of $\mathcal{F}$ is a measurable function $F$ satisfying $F(u)\geq\sup_{f\in \mathcal{F}}|f(u)|$. Finally, we let $\mathcal{Q}$ denote the set of probability measures with finite support on $\mathcal{Y}$.

Assumptions (ref) and (ref) are exactly the same as the conditions often imposed with i.i.d. data to show uniform LLNs and CLTs vanderVaart2000.\footnote{In vanderVaart2000, the supremum in Assumptions (ref) and (ref) is taken over the set of probability measures $Q$ with finite support on $\mathcal{Y}$ and such that $||F||_{Q,2}>0$. This additional restriction is simply due to a different convention in constructing covering numbers, as vanderVaart2000 considers open balls while we use closed balls, following, e.g., Kato2019.} In particular, Assumption (ref)-(i) (resp. (ii)) imposes a condition on what is usually referred to as the uniform (resp. bracketing) entropy integral, see, e.g., vanderVaartWellner1996. Finiteness of the uniform entropy integral is satisfied by any VC-type class of functions cherno2014, or by the convex hull of such classes under some restrictions. The bracketing entropy integral is finite for instance for classes of monotone or H\"older continuous functions vanderVaartWellner1996.

The following theorem establishes uniform LLNs and CLTs under these two conditions. We denote by $\bm{1}'$ the $k-$tuple $(1,k+1,...,2k-1)$.

The proof is in Section (ref) of the supplement. When Assumption (ref)-(ii) holds, Part 1 can be proved by essentially combining Theorem 3 in Eaglesonweber78 and Lemma 7.35 in kallenberg05. Part 2 was also proved for a finite $\mathcal{F}$ by silverman1976limit. But the weak convergence result under the bracketing entropy condition, and the uniform laws under the uniform entropy conditions, do not follow from such results. To prove the former, we adapt a maximal inequality in GineNickl2015 (2015, see their Lemma 3.5.12) to our context. To this end, we show that Hoeffding's bound on U-statistic Hoeffding1963 still applies to our context.

To prove the results under the uniform entropy conditions, the key ingredient, as with i.i.d. data, is a symmetrization lemma stated in Appendix (ref) below and proved in the supplement. Its proof relies extensively on Lemma (ref) and a decoupling inequality that may be of independent interest (see Lemma (ref)). The latter result generalizes a similar inequality for U-processes de1992decoupling. In the proofs of both lemmas, we follow similar strategies as with U-processes, with two complications. First, even with $k=2$, $Y_{\bm{i}}$ does not only depend on $U_{i_1}$ and $U_{i_2}$, but also on $U_{\{i_1,i_2\}}$. Second, when $k\geq 3$, dependence between observations arises not only because of single-unit terms such as $U_{i_1}$ or $U_{i_2}$, but also because of multiple-unit terms such as $U_{\{i_1,i_2\}}$.

As in the i.i.d. case, Assumption (ref) is actually stronger than necessary to obtain the uniform law of large numbers. The following proposition gives an exact characterization, where, for simplicity, we restrict to $k=2$. It is similar to the characterization for i.i.d. data GineNickl2015 or for U-processes delapena1999. Let us introduce the following norms:

Proposition (ref) emphasizes the two aspects of dissociated, exchangeable arrays. The first is i.i.d. variations, through the random entropy term related to $||\cdot||_{1,2}$, which only involves $(U_{\{i_1,i_2\}})_{\bm{i}\in\mathbb I_{n,2}}$. The second is U-statistic like variations, through the random entropy term related to $||\cdot||_{1,1}$: up to negligible terms, $||f||_{1,1}$ only depends on $(U_{i_1})_{1\leq i_1\leq n}$. Key in establishing the necessity of these two conditions is a weak converse of the symmetrization lemma for $k=2$, see Equation (ref) in the supplement.

Convergence of the bootstrap process

We now study the properties of the following bootstrap sampling scheme, which extends the pigeonhole bootstrap mccullagh2000resampling,owen2007pigeonhole to jointly separable arrays:

Then we consider $\mathbb{P}^{\ast}_{n}$ and $\mathbb{G}_{n}^{\ast}$, defined on $\mathcal{F}$ by $$\mathbb{P}_{n}^{\ast}f=\frac{(n-k)!}{n!}\sum_{\bm{i} \in \mathbb I_{n,k}} W_{\bm{i}} f(Y_{\bm{i}}),$$ $$\mathbb{G}^{\ast}_{n}f=\sqrt{n} \left(\mathbb{P}_{n}^{\ast}f - \mathbb{P}_{n}f \right).$$ Asymptotic validity of the bootstrap amounts to showing that conditional on the data $(Y_{\bm{i}})_{\bm{i} \in \mathbb I_k}$, $\mathbb{G}_n^{\ast}$ converges weakly to the process $\mathbb{G}$ defined in Theorem (ref).\footnote{For the sake of brevity, we focus afterwards on convergence results under the sole uniform entropy condition (Assumption (ref)-(i)).} As discussed in, e.g., vanderVaartWellner1996 (1996, Chapter 3.6), the outer almost-sure conditional weak convergence boils down to proving

where $\text{BL}_1$ is the set of bounded and Lipschitz functions from $\ell^{\infty}(\mathcal{F})$ to $[0,1]$ and “$\stackrel{\text{as}*}{\longrightarrow}$” denotes outer almost-sure convergence.

This theorem ensures the asymptotic validity of the bootstrap above not only for sample means, but also for smooth functionals of the empirical cdf and nonlinear estimators, as we shall see below. The proof of Theorem (ref), in Section (ref) of the supplement, follows the same lines as that of Theorem (ref), though some of the corresponding steps are more involved, as often with the bootstrap. In particular, to prove pointwise convergence, we use arguments in Lindeberg's proof of the CLT for triangular arrays, Theorem (ref).1 and Urysohn's subsequence principle, combined with Prohorov's theorem.

Note that in contrast with the standard bootstrap for i.i.d. data, $$\mathbb E\left(\mathbb{P}^{\ast}_{n}(f)\big| (Y_{\bm{i}})_{\bm{i} \in \mathbb I_k}\right)= \frac{1}{n^k}\sum_{\bm{i} \in \mathbb I_{n,k}}f(Y_{\bm{i}})\neq \mathbb P_nf.$$ However, the difference between $\mathbb P_n$ and $\mathbb P'_n$, the empirical measure with weights $1/n^k$, becomes negligible as $n\rightarrow\infty$. Accordingly, we also show in the proof of Theorem (ref) the almost-sure conditional convergence of $\sqrt{n} \left(\mathbb{P}_{n}^{\ast}f - \mathbb{P}'_{n}f \right)$, in addition to that of $\mathbb G^*_n$.

Application to nonlinear estimators

Theorem (ref) ensures consistency and asymptotic normality of a large class of estimators. In turn, Theorem (ref) shows that using the bootstrap for such estimators is asymptotically valid. To illustrate these points, we consider here two popular classes of estimators, namely Z-estimators and smooth functionals of the empirical cdf. Similar results could be obtained for, e.g., M-estimators cheng2010 or generalized method of moments estimators hansen1982large.

Let us first consider Z-estimators. Let $\Theta$ denote a normed space, endowed with the norm $\|\cdot\|_{\Theta}$ and let $(\psi_{\theta,h})_{(\theta,h)\in \Theta\times \mathcal{H}}$ denote a class of real, measurable functions. Let $\Psi(\theta)(h)=P\psi_{\theta,h}$, $\Psi_n(\theta)(h)=\mathbb P_n\psi_{\theta,h}$ and $\Psi^*_n(\theta)(h)=\mathbb P^*_n\psi_{\theta,h}$. We let, for any real function $g$ on $\mathcal{H}$, $\|g\|_{\mathcal{H}}=\sup_{h\in \mathcal{H}} |g(h)|$. The parameter of interest $\theta_0$, which satisfies $\Psi(\theta_0)=0$, is estimated by $\widehat{\theta}=\arg\min_{\theta\in\Theta} \|\Psi_n(\theta)\|_{\mathcal{H}}$. We also define $\widehat{\theta}^*=\arg\min_{\theta\in\Theta} \|\Psi^*_n(\theta)\|_{\mathcal{H}}$ as the bootstrap counterpart of $\widehat{\theta}$. The following theorem extends Theorem 13.4 in Kosorok2006 to jointly exchangeable and dissociated arrays. For related results on Z-estimators in the i.i.d. case, see Section 3.2 in vanderVaartWellner1996 and wellner1996.

Next, we consider smooth functionals of $F_Y$, the cdf of $Y_{\bm{i}}$. Suppose that $\mathcal{Y}\subset \mathbb R^p$ for some $p\in \mathbb N^+$ and $\theta_0=g(F_Y)$, where $g$ is Hadamard differentiable vanderVaartWellner1996. We estimate $\theta_0$ with $\widehat{\theta}=g(\widehat{F_Y})$, where $\widehat{F_Y}$ denotes the empirical cdf of $(Y_{\bm{i}})_{\bm{i} \in\mathbb I_{n,k}}$. Finally, we let $\widehat{\theta}^*$ denote the bootstrap counterpart of $\widehat{\theta}$.

In practice, $\mathbb{D}_0$ often corresponds to the set of functions that are continuous everywhere or at a certain point $y_0$. This is the case for instance with $g:F_Y\mapsto F_Y^{-1}(\tau)$ for $\tau \in (0,1)$. In such cases, one can show that $\mathbb{G} \in \mathbb{D}_0$ under the same condition as for i.i.d. data, namely that $F_Y$ is continuous everywhere or at the point $F_Y^{-1}(\tau)$.

Extensions

We now consider several extensions to our main results. First, we study the asymptotic behavior of the properly normalized empirical process in degenerate cases where $K(f,f)=0$. Second, we establish additional results on the bootstrap. Third, we study separately, rather than jointly, separable arrays. Other extensions to arrays with multiple observations per $k$-tuple and arrays where $Y_{\bm{i}}$ is defined even if there are identical indices in $\bm{i}$ are considered in the supplement. We also develop therein a test that the data are in fact i.i.d.

Degenerate cases

We consider here situations where $K(f,f)=0$ for all $f\in\mathcal{F}$, focusing for simplicity on $k=2$.\footnote{If $K(f,f)=0$ for only some $f\in\mathcal{F}$, we focus on $\mathcal{F}'=\{f\in \mathcal{F}: K(f,f)=0\}$.} Such a degeneracy appears for instance if the variables in the array are actually i.i.d., in which case $\sqrt{n}\mathbb G_n$ converges to a Gaussian process with covariance kernel $K(f_1,f_2)=\mathbb{C}ov(f_1(Y_{1,2}),f_2(Y_{1,2}))$. As another example menzel2017,bretagnolle1983, suppose that $Y_{i_1,i_2}=X_{i_1}X_{i_2}$, with $(X_i)_{i\in\mathbb N^+}$ i.i.d. variables with $\mathbb E(X_1)=0$, $\mathbb V(X_1)=1$. Let also $\mathcal{F}=\{f_\lambda(x)=\lambda x, \lambda\in I\}$ for a compact $I\subset\mathbb R$. Then one can easily see that $\sqrt{n}\mathbb G_n$ converges weakly in $\ell^\infty(\mathcal{F})$ to $\mathbb G(f_\lambda)=\lambda (Z^2 - 1) $, with $Z$ a standard normal variable.

More generally and as with U-processes ArconesGine1993, when $K(f,f)=0$, the rate of convergence of $\mathbb P_nf-Pf$ is $n^{-1}$ rather than $n^{-1/2}$ and the asymptotic distribution may not be normal. For any $(i_1,i_2)\in\mathbb I_2$, let $Y_{i_1,i_2}=\tau(U_{i_1},U_{i_2}, U_{\{i_1,i_2\}})$ be the Aldous-Hoover-Kallenberg representation where, without loss of generality, the variables in $\tau(\cdot,\cdot,\cdot)$ are assumed to be uniform on $[0,1]$. Let $\psi_m(u)=\left(1+\mathds{1}_{\{m\geq 2\}}\right)^{1/2}\cos\left(m\pi u\right)$ for $m$ even and $\psi_m(u)=\sqrt{2}\sin((m+1)\pi u)$ for $m$ odd. Then $(\psi_m)_{m\in\mathbb N}$ forms an orthonormal basis of $L^2[0,1]$. For all $\bm{m}\in\mathbb N^3$ and any $f\in\mathcal{F}$, we define $\mu_{\bm{m}}(f)$ by $$\mu_{\bm{m}}(f)=\mathbb E\left[\left[f(Y_{1,2})-\mathbb E\left(f(Y_{1,2})\right)\right]\psi_{m_1}(U_{1})\psi_{m_2}(U_{2}) \psi_{m_3}(U_{\{1,2\}})\right].$$ Let $(Z_m)_{m \in \mathbb N^+}$, $(Z_{m_1,m_2})_{(m_1,m_2)\in\mathbb N\times\mathbb N^+}$ and $(Z_{\{m_1,m_2\},m_3})_{(m_1,m_2,m_3)\in\mathbb N^2\times\mathbb N^+:m_1<m_2}$ denote independent standard normal variables. We then define the process $\mathbb G^d$ on $\mathcal{F}$ by

To prove the convergence of $\sqrt{n}\mathbb G_n$, we consider a condition on $\mathcal{F}$ that slightly differs from Assumption (ref)-(i).

Assumption (ref) is more stringent than Assumption (ref)-(i). A similar condition was also imposed by ArconesGine1993 for degenerate U-processes of order 1, see their condition (5.1).

As with degenerate U-processes ArconesGine1993, the limit process is a Gaussian chaos process. The result is based in particular on a symmetrization lemma and a maximal inequality taylored to these degenerate cases. Specifically, the symmetrized process only includes Rademacher variables at the pair $\{i_1,i_2\}$ level, or products $\varepsilon_{i_1}^{(1)}\varepsilon_{i_1}^{(2)}$ of Rademacher variables. We refer to Lemmas S(ref) and S(ref) in the supplement for more details.

Finally, we note that the bootstrap process considered above does not generally converge to $\mathbb G^d$.\footnote{The same holds true for the multiplier bootstrap process considered below.} With i.i.d. data, for instance, one can show that the variance of the bootstrapped mean converges to $3\mathbb V(Y_{i_1,i_2})$. We expect similar phenomena as with U statistics, where the bootstrap is known to fail in degenerate cases Arcones1992, Arcones1994. In the close case of separately exchangeable arrays (see Section (ref) below), menzel2017 shows that a suitable wild bootstrap is consistent for the sample average, whether or not we have degeneracy. Whether such a result generalizes to the empirical process is left for future research.

Further results on the bootstrap

Theorem (ref) shows convergence of the bootstrap process under conditions on $\mathcal{F}$ that ensure the convergence of the initial process $\mathbb G_n$. The following result shows that under moment conditions, convergence of $\mathbb G_n$ is actually necessary for the convergence of $\mathbb G_n^*$ to a Gaussian process.

Theorem (ref) may be seen as a partial extension to jointly exchangeable arrays of Theorem 2.4 in gine1990, which, with i.i.d. data, establishes the equivalence between the convergence of the bootstrap process and $PF^2<\infty$ together with convergence of the initial process.

With i.i.d. data, several other bootstrap schemes than the multinomial bootstrap are possible: see, e.g., barbe1995 for an extensive review. The situation is probably no different with jointly exchangeable arrays. To illustrate this, we consider a version of the multiplier bootstrap adapted to such data Kosorok2003. Specifically, let $(\xi_i)_{i=1}^n$ be a sequence of i.i.d. random variables that are centered, have unit variance and are independent from the original data $(Y_{\bm{i}})_{\bm{i}\in\mathbb I_{n,2}}.$ We then consider the following process: $$\mathbb{G}_n^{m*} : f \mapsto \frac{1}{\sqrt{n}}\sum_{i_1=1}^n\xi_{i_1}\left(\frac{1}{n-1}\sum_{1\leq i_2\neq i_1 \leq n} \left[f(Y_{i_1,i_2})+ f(Y_{i_2,i_1})\right]- 2 \mathbb{P}_nf\right).$$

The next theorem shows the conditional weak convergence of $\mathbb{G}_n^{m*} $ under the same conditions on $\mathcal{F}$ as previously.

Separately exchangeable arrays

Up to now, we have considered cases where the $n$ units that interact stem from the same population. In some cases, however, they do not, because the $k$ populations differ. For instance, we may be interested only in relationships between men and women. In that case, the symmetry condition in Assumption (ref) has to be strengthened: both the labelling of men and the labelling of women should be irrelevant. This corresponds to so-called separately exchangeable arrays, defined formally in Assumption (ref) below. Another important motivation for considering separately exchangeable arrays is multiway clustering, namely dependence arising through different dimensions of clustering. For instance, wages of workers may be affected by local shocks or sector-of-activity shocks. In such cases, we observe $Y_{i_1,i_2}$, the wage of a worker in geographical area $i_1$ and sector of activity $i_2$.\footnote{Oftentimes, we actually have several observations per cell, and the number varies from one cell to another. This extension is discussed in Section (ref) of the supplement.}

More generally, we consider in this section random variables $Y_{\bm{i}}$ where $\bm{i}=(i_1,...,i_k)\in \mathbb N^{+k}$, implying that repetitions (e.g. $\bm{i}=(1,...,1)$) are allowed. We impose the following condition on these random variables.

This condition is stronger than Assumption (ref) since it implies in particular equality in distribution for $\pi_1=...=\pi_k$.

Let us redefine $\boldsymbol{1}$ here as $(1,...,1)$ and let $\bm{n}=(n_1,...,n_k)$, where $n_j\geq 1$ denotes the number of units observed in population $j$ (or cluster $j$ with multiway clustering). Note that in general, $n_j \neq n_{j'}$ for $j\neq j'$. The sample at hand is then $(Y_{\bm{i}})_{\boldsymbol{1}\leq \bm{i}\leq \bm{n}}$, where $\bm{i}\geq \bm{i}'$ means that $i_j\geq i'_j$ for all $j=1,...,k$. Let $\underline{n}=\min(n_1,...,n_k)$. The empirical measure and empirical process that we consider for separately exchangeable arrays are:

We also consider the “pigeonhole bootstrap”, suggested by mccullagh2000resampling and studied, in the case of the sample mean and for particular models, by owen2007pigeonhole. This bootstrap scheme is very close to the one we considered in Section (ref) for jointly exchangeable arrays, except that the weights are now independent from one coordinate to another:

The bootstrap process $\mathbb G_{\bm{n}}^{\ast}$ is thus defined on $\mathcal{F}$ by $$\mathbb G^{\ast}_{\bm{n}}f=\sqrt{\underline{n}} \left(\frac{1}{\prod_{j=1}^k n_j}\sum_{\boldsymbol{1} \leq \bm{i} \leq \bm{n}}\left(W_{\bm{i}} - 1\right)\sum_{\ell=1}^{N_{\bm{i}}}f(Y_{\bm{i},\ell})\right).$$

Henceforth, we consider the convergence of $\mathbb P_{\bm{n}}$, $\mathbb G_{\bm{n}}$ and $\mathbb G^*_{\bm{n}}$ as $\underline{n}$ tends to infinity. More precisely, as with multisample U-statistics vanderVaart2000, we assume that there is an index $m\in \mathbb N^+$, left implicit hereafter, and increasing functions $g_1,...,g_k$ such that for all $j$, $n_j=g_j(m)\rightarrow \infty$ as $m\rightarrow \infty$ (we also assume without loss of generality that for all $m\in\mathbb N^+$, $g_j(m+1)>g_j(m)$ for some $j$). The following theorem extends Theorems (ref) and (ref) to this set-up.

Theorem (ref) includes the case where $\lambda_j=0$ for some $j$, corresponding to “strongly unbalanced” designs with different rates of convergence to $\infty$ along the different dimensions of the array. In that case, only the dimensions with the slowest rate of convergence contribute to the asymptotic distribution, as can be seen in (ref).

Because the $(n_j)_{j=1...k}$ are not all equal in general, Theorem (ref) does not follow directly from Theorem (ref), even if Assumption (ref) is stronger than Assumption (ref). We prove the result by showing a simpler and convenient version of the symmetrization lemma in this setting. We refer to Lemma S(ref) in the supplement for more details.

Applications to international trade

Finally, we illustrate the importance of accounting for dependence in real dyadic data, through two applications to international trade data.

Evolution of international trade

There is a large interest in economics on the evolution of international trade. But before analyzing the causes and consequences of such an evolution, one must check that there is indeed some significant changes. In this first application, we test whether the distribution of exports remains the same between two consecutive years, using Comtrade data on all countries from 2012 to 2018. We use for that purpose the Kolmogorov-Smirnov (KS) test statistic $$KS_t=\sup_{u\in\mathbb R}\left|\frac{1}{n(n-1)}\sum_{(i_1,i_2)\in\mathbb I_{n,2}} \mathds{1}_{\{T_{i_1,i_2,t}\leq u\}} - \mathds{1}_{\{T_{i_1,i_2,t+1}\leq u\}}\right|.$$ where $T_{i_1,i_2,t}$ denotes the trade volume from country $i_1$ to country $i_2$ in year $t$. Let us assume that Assumption (ref) holds, with $Y_{\bm{i}}=(T_{\bm{i},t}, T_{\bm{i},t+1})$. Then, under the null hypothesis that the distributions of $T_{\bm{i},t}$ and $T_{\bm{i},t+1}$ are equal, we have, by Theorem (ref), $\sqrt{n} KS_t \stackrel{d}{\longrightarrow} \|\mathbb{G}\|_{\mathcal{F}}$, with $\mathcal{F}=\{f_u(x,y)=\mathds{1}_{\{x\leq u\}}-\mathds{1}_{\{y\leq u\}}\}$. Given the dependence structure both between pairs of countries and across time, the distribution of $\|\mathbb{G}\|_{\mathcal{F}}$ depends on the true data generating process. To estimate it, we rely on the recentered bootstraped test statistic: $$KS^*_t=\sup_{u\in\mathbb R}\left|\frac{1}{n(n-1)}\sum_{(i_1,i_2)\in\mathbb I_{n,2}} (W_{\bm{i}} -1) \left(\mathds{1}_{\{T_{i_1,i_2,t}\leq u\}}- \mathds{1}_{\{T_{i_1,i_2,t+1}\leq u\}}\right)\right|.$$ We compute the p-value of the test by $\mathbb P\left(KS^*_t>KS_t \big| (Y_{\bm{i}})_{\bm{i} \in\mathbb I_{n,k}}\right)$. For the sake of comparison, we also compute p-values based on alternative forms of dependence that have been considered in applied work on similar data. Specifically, we also assume that the variables $(Y_{\bm{i}})_{\bm{i}}$ are i.i.d. We then assume pairwise clustering, where $Y_{i_1,i_2}$ and $Y_{i_2,i_1}$ may be dependent, but $Y_{\bm{i}}$ and $Y_{\bm{j}}$ are independent if $\bm{j}$ is not a permutation of $\bm{i}$. We also consider one-way clustering according to $i_1$ (and, similarly, according to $i_2$). In this case, $Y_{i_1,i_2}$ and $Y_{i_1,i_3}$ may be dependent, but $Y_{i_1,i_2}$ and $Y_{i'_1,i_3}$ are independent as soon as $i_1\neq i'_1$, whether or not $i_2=i_3$. For each of these cases, we use the bootstrap, but with different bootstrap schemes accounting for these different dependence structures.

The results are displayed in Table (ref). They suggest significant changes in export volumes in some years but not all. In particular, international trade seems very stable between 2015 and 2017. There is some evidence of changes between 2012 and 2015 but we still do not reject the null hypothesis at the 1% level for the years 2013-2014. The other columns of the table shows the importance of accounting for dependence along both dimensions. In particular, assuming i.i.d. data or pairwise dependence always leads to a strong rejection of the null, except for 2015-2016.\footnote{ A concern is that if the data are actually i.i.d. (or, more generally, pairwise dependent), our bootstrap is conservative, which would explain the discrepancy between the p-values under pariwise dependence and non-degenerate joint exchangeability. Using the methodology in Section (ref) of the supplement, we test for pairwise dependence. For the eight years we consider, the null hypothesis is rejected at all standard levels, with p-values always smaller than $10^{-4}$.} Clustering along exporters also leads to artificially small p-values, in particular for the pairs 2013-2014, 2014-2015 and 2016-2017. In this context, clustering along importers leads to results that are closer to those based on dyadic data.

Estimation of a gravity equation

Second, we revisit silva2006log, who estimate the so-called gravity equation for international trade. Omitting the year index, this gravity equation states that $T_{i_1,i_2}$ satisfies

where $G_{i}$ denotes country $i$'s GDP, which would correspond to the mass of $i$ in a traditional gravity equation, $D_{i_1,i_2}$ denotes the distance between $i_1$ and $i_2$, $A_{i_1,i_2}$ are additional control variables and $\eta_{i_1,i_2}$ is an unobserved term.

To estimate $\theta_0 = (\alpha_0,...,\alpha_3, \beta')'$, silva2006log suggest to use the Poisson pseudo maximum likelihood (PPML for short) estimator $\widehat{\theta}$. The idea, formalized in gourieroux1984pseudo, is that with i.i.d data, the PPML estimator is consistent and asymptotically normal for $\theta_0$ even if $T_{\bm{i}}$ does not follow a Poisson model, provided that $\mathbb E\left[\eta_{\bm{i}}|X_{\bm{i}}\right]=1$, with $X_{\bm{i}}=(1,\ln(G_{i_1}),\ln(G_{i_2}),\ln(D_{\bm{i}}),A_{\bm{i}})$. This is because the PPML estimator is based on the empirical counterpart of

and this equality holds true if $\mathbb E\left[\eta_{\bm{i}}|X_{\bm{i}}\right]=1$.

Now, assuming as in silva2006log that the variables $(Y_{\bm{i}})_{\bm{i} \in \mathbb I_2}$ (with $Y_{\bm{i}}=(T_{\bm{i}},X_{\bm{i}})$) are i.i.d. is restrictive. We suppose instead that Assumption (ref) holds. Then Theorem (ref) applies to this setting, implying that $\widehat{\theta}$ is still consistent and asymptotically normal in this case.\footnote{In this case, $\mathcal{H}=\{1,...,\text{dim}(X_{\bm{i}})\}$ and $\psi_{\theta,h}(Y_{\bm{i}})=X_{h,\bm{i}}(T_{\bm{i}} - \exp(X_{\bm{i}}\theta_0))$. Then the key conditions 2 and 3 in Theorem (ref) are satisfied as soon as $\Theta$ is bounded, see e.g. Example 19.7 in vanderVaart2000.} Nonetheless, the rates of convergence and asymptotic variance are different in the two cases, resulting in different inference on $\theta_0$.\footnote{The same application has been considered by graham2018, who shows, assuming convergence of a certain sample average, the asymptotic normality of the PPML estimator under the same dependence structure as ours. On the other hand, he neither considers bootstrap-based inference nor proves the consistency of his (asymptotic) variance estimator.}

We use the same dataset as silva2006log, which covers 136 countries for year 1990, and consider the exact same specification as the one they use in their Table 3. In this specification, the additional control variables $A_{\bm{i}}$ include exporter- and importer-level variables, namely their GDP per capita, a dummy variable equal to one if countries are landlocked and a remoteness index, which is the log of GDP-weighted average distance to all other countries. It also includes variables at the pair level, namely dummy variables for contiguity, common language, colonial tie, free-trade agreement and openness. This openness dummy is equal to one if at least one country is part of a preferential trade agreement. We refer to silva2006log for additional details.

Table (ref) below presents the results. The first column displays the point estimates, which, as expected, are identical to those in silva2006log. The other columns display the p-values for the null hypothesis that $\theta_{0j}$, the $j$-th component of $\theta_0$, is equal to 0. We consider the same forms of dependence as with the KS test above. Under joint exchangeability, we compute the p-value $p_j$ for $\theta_{0j}=0$ using $p_j=\mathbb{P}\left(|\widehat{\theta}_j^* - \widehat{\theta}_j|>|\widehat{\theta}_j| \big| (Y_{\bm{i}})_{\bm{i} \in\mathbb I_{n,k}}\right)$. For other forms of dependence, we follow the usual practice of computing the p-values using the asymptotic normality of $\widehat{\theta}_j$ and estimators of the asymptotic variance under these various dependence structures.

Using our bootstrap leads to much larger p-values than under the i.i.d. assumption. Only the log of distance and the log of GDP of the exporter and the importer appear to be significant at the $10^{-3}$ levels, whereas five additional control variables are significant at that level under the i.i.d. assumption. In particular, common language and importer's remoteness are not even significant at the usual 5% level.\footnote{ As in Footnote (ref) above, we test for pairwise dependence, to see whether our results could be driven by the fact that our bootstrap is conservative in such cases. We obtain a p-value smaller than $10^{-4}$ and thus reject this hypothesis at all usual levels.} Interestingly, there is also a gap between assuming one-way clustering, either at the exporter or at the importer level, and assuming to have a jointly exchangeable and dissociated array. In the former case, we still have seven variables that are significant at the $10^{-3}$ levels. Confidence intervals, not displayed here, lead to similar conclusions. In particular, compared to the average length of i.i.d.-based 95% confidence intervals, those based on pairwise clustering are only 8% wider. Those based on one-way clustering on exporters (resp. importers) are 20% (resp. 17%) larger. On the other hand, those based on Assumption (ref) are 136% wider.

Conclusion

While polyadic data are increasingly used in applied work, and empirical researchers routinely account for multiway clustering when computing standard errors, the statistical theory behind these forms of dependence has lagged behind. Following bickel2009nonparametric and menzel2017, we link these dependence structures to jointly and separately exchangeable arrays. Using representation results for such arrays, we then prove uniform laws of large numbers and central limit theorems. These results imply consistency and asymptotic normality of various nonlinear estimators under such dependence. We also establish the general validity of natural extensions of the standard nonparametric bootstrap to such arrays. Our application shows that using those bootstrap schemes may make a large difference compared to assuming i.i.d. data or clustering along a single dimension, as has often been done.

One caveat is that for the bootstrap confidence intervals to be valid, the asymptotic variance of the estimator should be positive. This may not be the case, for instance if the data $(Y_{\bm{i}})_{\bm{i}\in \mathbb I_k}$ are actually i.i.d. Inference based on the wild bootstrap without this positivity condition has been studied for sample averages under multiway clustering by menzel2017. How to conduct inference on nonlinear estimators under joint exchangeability or multiway clustering without this positivity condition remains an avenue for future research.