# Inference in Linear Dyadic Data Models
with Network Spillovers††thanks: We thank Aimee Chin, Hugo Jales, Taisuke
Otsu, Bent E. Sørensen, Vitor Possebom and seminar participants at the
University of Houston and at the 2022 Texas Econometrics Camp for valuable
comments and suggestions. The authors report there are no competing interests
to declare.

Nathan Canen Department of Economics, University of Houston, Science Building
3581 Cullen Boulevard Suite 230, Houston, 77204-5019, Texas, United States.
Email<EMAIL_ADDRESS>Ko Sugiura Department of Economics, University of
Houston, Science Building 3581 Cullen Boulevard Suite 230, Houston,
77204-5019, Texas, United States. Email<EMAIL_ADDRESS>
When using dyadic data (i.e., data indexed by pairs of units), researchers
typically assume a linear model, estimate it using Ordinary Least Squares and
conduct inference using “dyadic-robust” variance estimators. The latter
assumes that dyads are uncorrelated if they do not share a common unit (e.g.,
if the same individual is not present in both pairs of data). We show that
this assumption does not hold in many empirical applications because indirect
links may exist due to network connections, generating correlated outcomes.
Hence, as we formally prove and illustrate in Monte Carlo simulations,
“dyadic-robust” estimators can be severely biased in such situations. We
develop a consistent variance estimator appropriate for such contexts by
leveraging results in network econometrics. Our estimator has good finite
sample properties in numerical simulations. We then illustrate our message
with an application to politicians’ voting behavior when they are seating
neighbors in the European Parliament.

Keywords: Dyadic data, Networks, Inference, Cross-sectional dependence,
Congressional Voting.

JEL Codes: C12, C21, C40, D72

##  1 Introduction

Dyadic data is categorized by the dependence between two sets of sampled units
(dyads). For example, exports between the U.S. and Canada depend on both
countries (and, plausibly, their characteristics), while arrangements between
a seller and a buyer show similar dependence. This contrasts to classical data
in the social sciences that only depends on a single unit of observation
(e.g., the GDP of the U.S.).

The empirical relevance of dyadic data is showcased by its widespread use,
which has increased over the past two decades (Graham (2020a) provides an
extensive review). For example, applications are found in international trade
(export-import outcomes across countries, Anderson and van Wincoop (2003)),
macroeconomics (exchange rate determination, Lustig and Richmond (2020)),
development (gift giving and insurance between two members in a village,
Fafchamps and Gubert (2007)), political economy (correlation in voting
behavior in Parliament across seating neighbors, Harmon et al. (2019)), among
many others. In all of these examples, applied researchers model the
dependence between dyadic outcomes and observable characteristics using a
linear model, which they then estimate using Ordinary Least Squares (OLS).
However, inference on such estimators for the linear parameters is more
complex.

The main approach in applied work has been the use of the so-called “dyadic-
robust” estimators (e.g., Cameron et al. (2011), Aronow et al. (2015), and
Tabord-Meehan (2019), among others). Such estimators build on the widely-used
assumption in dyadic data that the error terms for dyad $(i,j)$ and for dyad
$(k,l)$ can only be correlated if they share a unit (see Aronow et al. (2015)
and Tabord-Meehan (2019) for a discussion; Cameron and Miller (2014) for a
review). In the context of international trade, this is equivalent to
assuming, for example, that Canada-U.S. and France-U.S. trade may be
correlated because the U.S. is present in both, but there is no correlation
between Canada-U.S. exports and France-Germany ones (conditional on observable
characteristics).

In this paper, we first argue that such an assumption does not hold in many
applications using dyadic data where dyads may be indirectly connected along a
network.111This is a concrete class of applied examples where the assumption
fails. The possibility that cross-sectional dependence in dyadic data might be
more extensive than assumed has been pointed out by Cameron and Miller (2014).
Figure 1 (A) presents a simple example in the context of international trade.
That is, it is completely possible that Canada-U.S. and France-Germany trades
might be correlated along unobservables because they have many indirect
connections (in the figure, through Canada trading with China which then
trades with France).222Such spillovers could be rationalized as individual-
level unobserved heterogeneity: e.g., an unmeasured export orientation of an
economy and unmodeled tastes for consumption in the literature of
international trade (see Graham (2020b) and references therein for details).
This implies that the use of dyadic-robust variance estimators may be
inappropriate for inference in linear models with dyadic data when spillovers
exist along a network.

Figure 1: Hypothetical Example of Trade Network

(a) Trade network across countries

(b) Trade network across pairs of countries

  * Note: The left figure shows a hypothetical example of trade networks across countries. Each double-headed arrow signifies the existence of a bilateral trade flow between countries. (The absence of a double-headed arrow is interpreted as the nonexistence of any type of trade flow between the countries.) In this example, though the Canada-U.S. and France-Germany trades (denoted by black arrows) do not have a country in common, they might still be correlated through two indirect links (denoted by gray arrows): namely, the Canada-China-Germany trade, and the Canada-China-France trade. The right figure illustrates the network across pairs of countries (i.e., dyads). Round squares stand for bilateral trade flows between countries and the adjacency of these trade flows are represented by turned square arrows.

We make this point in two ways. First, we prove that existing “dyadic-robust”
variance estimators are biased for the true asymptotic variance when there are
network spillovers (i.e., when dyads may be correlated along a network even
when they do not share a common unit). This result is then verified in Monte
Carlo simulations with standard set-ups, where we show that such variance
estimators lead to a bias of approximately $14\%$ and its associated
confidence interval undercovers the true parameter by over 9 percentage
points.333These numbers are obtained from the simulation results for the case
where the underlying network is modestly complex (the Barabási and Albert’s
(1999) random network with parameter $\nu=3$), the sample size is large
($N=5000$) and the spillovers are significant ($\gamma=0.8$). See Section 5
for details. This issue may be prevalent in many applied examples. In the case
of international trade (e.g., Cameron et al. (2011) and Cameron and Miller
(2014)), a pair of trade flows may be correlated even if they do not share a
member country because there are many indirect connections among them. In
political economy, those seating far apart from one another in Parliament may
still have correlated behavior through the indirect links of sharing seating
neighbors. In Macroeconomics, currency exchange rates between two countries
may depend on other countries beyond those involved. Furthermore, we show that
the performance worsens with the density of the underlying network (e.g., in
political economy and international trade, networks are often dense, see
Cameron and Miller (2014)) and with sample size, wherein the role of indirect
spillovers appears to be prominent.

To deal with these issues, we develop a consistent variance estimator that can
account for such network spillovers even with dyadic data. This complements
existing dyadic-robust variance estimators (e.g., Aronow et al. (2015)) by
explicitly accounting for the network structure of spillovers.444We provide an
extensive comparison of the relative benefits of each approach in the next
section. We note here, though, that neither approach subsumes the other, as
they depend on different assumptions, and may be more appropriate for
different applications. We prove that our proposed variance estimator is
consistent for the true variance of the OLS estimator in linear models with
dyadic data when the cross-sectional dependence follows an observed
(exogenous) network. Our main insight is that the dependence across all dyads,
including indirect spillovers, can be rewritten as correlations across a
specific network over dyads. This allows us to apply the framework of
Kojevnikov et al. (2021) to such network random variables, although here it is
a network over dyads, rather than individuals.

Finally, we illustrate the extent to which neglecting network spillovers with
dyadic data may bias inference results. Beyond the Monte Carlo simulations
described above, we revisit the application in Harmon et al. (2019) of voting
in the European Parliament. The authors study whether random seating
arrangements (based on naming conventions) induce neighboring politicians to
agree with one another in policy votes. The outcome, whether politicians $i$
and $j$ vote the same way on a policy, is dyadic in nature. However, $i$ and
$j$’s votes may be positively correlated even if they are not neighbors: for
instance, $i$ and $j$ may sit on either side of common neighbors $k$ and $l$,
who influence them both, and this seating arrangement is observed. This chain
of influences is sufficient to induce strong positive correlation across non-
dyads. We show that neglecting such higher order spillovers has significant
empirical consequences: their estimated variance using Aronow et al. (2015) is
roughly $22\%$ smaller than using our consistent estimator accounting for such
spillovers; while the estimate based on the Eicker-Huber-White estimator
ignoring spillovers is approximately 73% smaller than our proposal.

###  1.1 Related Literature

Our approach directly relates to the literature on the regression analysis
based on dyadic random variables and data. Aronow et al. (2015) and Tabord-
Meehan (2019) consider OLS estimation and inference in a linear dyadic
regression model. Meanwhile, Graham (2020a) and Graham (2020b) explore a
likelihood-based approach to dyadic regression models, while Graham et al.
(2022) and Chiang and Tan (2022) provide results for kernel density estimation
in dyadic regression models. As alluded to previously, all such papers assume
that dyads are uncorrelated unless they share a common unit, thereby ignoring
dependence across indirectly linked dyads. On the other hand, our network-
robust variance estimator explicitly captures the higher order correlation
across dyads along the network over the dyads.

However, we emphasize that neither approach subsumes the other. The papers
cited above leave the dependence within “clusters” (groups of dyads that share
units) unrestricted, but assume independence across such clusters. This is
akin to the literature with one-way clustering (e.g., Hansen and Lee (2019)).
By comparison, our approach restricts such dependence among groups of dyads
that share units (i.e., dependence still follows the observed network), but
allows for dependence across such “clusters” of dyads along the dyadic
network. As a result, our approach complements the existing methods based on
dyadic clustering, similar to how inference with spatial data (e.g., Conley
(1999) and Jenish and Prucha (2009))555Our approach still differs from such
inference with spatial data, since the latter routinely assumes the index set
to be a Euclidean metric space, whose metric relies solely on the nature of
the space, and uses it to define the dependence between variables. See also
Ibragimov and Müller (2010) and references therein. In our context, however,
the index set of dyads alone does not suffice to dictate the dependence
structure because indices themselves do not inform us of the network topology;
instead, we first introduce a metric on a network among dyads and our mixing
condition is based on dependence as dyads grow further apart along the
network. complements one-way clustering inference. Which approach to pursue
using dyadic data depends on the researchers’ applications and how they fit
such assumptions.

Our main insight in accommodating such indirect spillovers is that we can
rewrite the correlation structure among dyads as a dyadic network, where links
denote whether they share a common member. As a result, this dyadic network
describes how close/far certain dyads are from sharing members with other
dyads. In doing so, the transformed problem is amenable to appropriate
applications of recent developments in the econometrics of random variables
which are correlated along an (observable, exogenous) network. In particular,
we apply asymptotic results for network-dependent random variables developed
by Kojevnikov et al. (2021)666Vainora (2020) provides another such theoretical
contribution. to an appropriately defined dyadic network, with assumptions
imposed on the latter. Leung (2021) and Leung (2022) also apply the framework
of Kojevnikov et al. (2021) to study, respectively, cluster-robust inference
and causal inference for the case of individual-specific random variables.
These papers focus on the correlation along a network over individuals, rather
than over dyads. Meanwhile, Leung and Moon (2021) derive an asymptotic theory
for dyadic variables in the context of networks, primarily for endogenous
network formation models.

Finally, our paper is also naturally related to the recent literature on
inference for multiway clustering, whether OLS estimation with multiway
clustering (Cameron et al. (2011)); a clustering method in high-dimensional
set-ups (Chiang et al. (2021)); clustering within the time dimension (Chiang
et al. (2022)); clustered inference with empirical likelihood (Chiang et al.
(2022)); bootstrap methods in multiway clustering (Davezies et al., 2021;
Menzel, 2021; MacKinnon et al., 2022b); clustering in the context of average
treatment effects (Abadie et al. (2022)), to name but a few (see MacKinnon et
al. (2022a) for review). Again, one of the common key assumptions that
underlies these studies is that individual observations are divided into
disjoint groups – clusters – and observations in different clusters are not
correlated.777MacKinnon et al. (2022c) propose empirically computable measures
of the cluster-level influence and leverage that can be used to judge to what
extent the underlying assumption of cluster-robust variance estimation is met.
While useful to allow for correlations along time and within such groups
(e.g., countries, or in microlevel data observing multiple villages over time,
as in Banerjee et al. (2013)), this separable structure may be inappropriate
for other environments where spillovers follow a complex form of dependence
along a network (e.g., as those in Figure 1 and our empirical application
below).

The rest of the paper is organized as follows. In the next section, we define
the linear model and additional notation necessary for our results. Section 3
presents results on the consistency and asymptotic normality of the OLS
estimator in this environment, necessary to prove the consistency of our
variance estimator. The latter is our main result, presented in Section 4 and
proved in Appendix. Section 5 presents Monte Carlo simulation results
comparing the finite sample performance of our estimator to the Eicker-Huber-
White estimator and the dyadic-robust one when there are network spillovers,
while Section 6 presents the empirical application. Section 7 concludes. In
Appendix, we provide further details on our set-up, proofs and additional
Monte Carlo exercises, and details on the empirical application.

##  2 Set-up

Assume that we observe a cross section of $N\in\mathbbm{N}$ individuals
located along a network – the latter interpretable as consumers, firms, and
other observation units depending on the context. The dyads present in the
$N$-individual network (i.e., among the $\tbinom{N}{N-1}$ possible dyads), are
called active dyads (e.g., importer-exporter, employer-employee
relationships). Let each active dyad between the units be indexed by
$m\in\mathcal{M}_{N}\subseteq\mathbbm{N}$ where $\mathcal{M}_{N}$ denotes an
index set of the active dyads among $N$ individuals with
$M\coloneqq|\mathcal{M}_{N}|\in\mathbbm{N}$. That is, the dyad for two units
$i$ and $j$ (e.g., firms, individuals) is denoted as some $m$.

We assume that each dyad $m$ is endowed with a triplet of the dyad-specific
variables, forming a triangular array
$\\{(y_{M,m},x_{M,m},\varepsilon_{M,m})\\}_{m\in\mathcal{M}_{N}}$ with respect
to $M$, where $y_{M,m}\in\mathbbm{R}$ is a one-dimensional observable outcome,
$x_{M,m}\in\mathbbm{R}^{K}$ is a $K$-dimensional vector of observable
characteristics with $K\in\mathbbm{N}$, and $\varepsilon_{M,m}\in\mathbbm{R}$
is a one-dimensional random error term that is not observable to the
econometrician. Here, it is assumed that the network structure among
individuals is exogenously determined and so is the network among dyads, which
is furthermore assumed to be conditionally independent of
$\\{\varepsilon_{M,m}\\}_{m\in\mathcal{M}_{N}}$. This means that we only
consider exogenous network formation. Note that the network is assumed to be
observed.

The subsequent arguments will require us to distinguish between a pair of
dyads who share a member (i.e., who are directly linked – which we call,
adjacent) and a pair of dyads who are directly or indirectly linked (which we
call, simply, connected).

######  Definition 2.1 (Adjacent & Connected Dyads).

Two active dyads $m$ and $m^{\prime}$ are said to be adjacent if they have an
individual in common; and they are called connected if there are other pairs
of adjacent dyads between them.

In Figure 1 (A), U.S.-Canada trade is adjacent to Canada-China, and connected
with, though not adjacent to, France-Germany. We emphasize that the adjacency
relationship constitutes a network structure among active dyads, and thus a
network over individuals can be transformed to the one over active dyads: for
example, Figure 1 (B) provides a network over pairs of trading countries
(i.e., active dyads).888This corresponds to thinking about the line graph of
the original graph over individuals. We define the geodesic distance between
two connected dyads $m$ and $m^{\prime}$ to be the smallest number of adjacent
dyads between them. Note that adjacent dyads are a special case of connected
dyads with the geodesic distance equal one.

It is worth noting that there are two different units at play: the number of
sampling units (i.e., individuals), $N$, and the number of dyads, $M$. The
latter constitutes the basis of statistical analysis. To bridge these two
units, we assume that the number of active dyads also goes to infinity as the
number of sampling units is taken to infinity. This is formalized by the
following assumption, which allows us to link large $N$ asymptotics to large
$M$ ones.

######  Assumption 2.1.

There exists a positive constant $c$ such that $M\geq cN$ whenever $N$ is
sufficiently large.

Assumption 2.1 requires the number of active dyads to grow at least at some
constant rate relative to the number of individuals, eliminating the
possibility of extremely sparse networks among the individuals. This is
similar in spirit to Assumption 2.3 of Tabord-Meehan (2019) in which the
minimum degree is assumed to grow at some constant rate relative to the number
of individuals.999Our assumption is milder than Assumption 2.3 of Tabord-
Meehan (2019) since the latter does not allow any individual to be isolated,
while Assumption 2.1 merely constrains the average degree. Similarly, this
assumption is weaker than the assumption that the maximum degree in a network
is bounded even when $N\rightarrow\infty$ (e.g., Penrose and Yukich (2003) and
de Paula et al. (2018)). We choose this viewpoint – large $N,M$ asymptotics –
in light of the empirical relevance of such situations. For example, in
international trade, entry of a new country/firm to a market will mostly
likely increase the number of trade flows in the economy; in political
economy, the more members of parliaments (MEPs) there are, the more pairs of
the MEPs sitting next to each other there will be (see Section 6). Finally,
this assumption by itself does not rule out the possibility that $M$ grows
faster than $N$.

###  2.1 The Linear Model

####  2.1.1 Set-up & Identification

The cross-sectional model of interest takes the form of the linear network-
regression model: for any $N\in\mathbbm{N}$,

$\displaystyle y_{M,m}=x_{M,m}^{\prime}\beta+\varepsilon_{M,m}\ \ \ \ \
\forall m\in\mathcal{M}_{N},$ (1)

where

$\displaystyle Cov(\varepsilon_{M,m},\varepsilon_{M,m^{\prime}}\mid
X_{M})=0\quad\text{unless $m$ and $m^{\prime}$ are connected},$ (2)

and $\beta$ is a $K\times 1$ vector of the regression coefficients and $X_{M}$
denotes the $M\times K$ matrix that records the observed dyad-specific
characteristics, i.e., $X_{M}\coloneqq[x_{M,1},\ldots,x_{M,M}]^{\prime}$.

We note that equation (2) allows for there to be spillovers across the error
terms even when dyads $m$ and $m^{\prime}$ are not adjacent, as long as they
are connected through indirect links. This specification has been of
substantial empirical interest (see, e.g., Harmon et al. (2019) and Lustig and
Richmond (2020)).

We impose the following standard assumptions for identification of $\beta$.

######  Assumption 2.2.

For each $N\in\mathbbm{N}$:
(a) $\sup_{m\in\mathcal{M}_{N}}E\big{[}|\varepsilon_{M,m}|^{2}\big{]}$ exists
and is finite;
(b) $\sup_{m\in\mathcal{M}_{N}}E\big{[}\|x_{M,m}\|\big{]}$ exists and is
finite;
(c) $E\big{[}x_{M,m}x_{M,m}^{\prime}\big{]}$ exists with finite elements and
positive definite for all $m\in\mathcal{M}_{N}$;
(d) $E\big{[}\varepsilon_{M,m}\mid X_{M}\big{]}=0$ for all
$m\in\mathcal{M}_{N}$.

Assumption 2.2 (a) and (b) are standard and jointly imply the finite existence
of the second moment of $y_{M,m}$ for all $m\in\mathcal{M}_{N}$, which in turn
implies the finite existence of the cross moment of $y_{M,m}$ and $x_{M,m}$
for all $m\in\mathcal{M}_{N}$. The third and fourth assumptions are also
standard in the context of the linear regression models and require no
multicolinearity and strict exogeneity, respectively.

Identification of the linear parameter in equation (1) follows from Assumption
2.2 (see Proposition B.1 in Appendix B.1).

Our setting encompasses dyadic models commonly studied in the econometrics
literature. To see this, it is convenient to write $m=d(i,j)$ to denote that
dyad $m$ is the dyad related to units $i$ and $j$.

######  Example 2.1 (Dyadic-Regression Models).

Aronow et al. (2015) and Tabord-Meehan (2019) consider a variant of the linear
regression (1) under the assumption

$\displaystyle Cov(\varepsilon_{M,d(i,j)},\varepsilon_{M,d(k,l)}\mid
x_{M,d(i,j)},x_{M,d(k,l)})=0\hskip 14.22636pt\text{unless
\quad$\\{i,j\\}\cap\\{k,l\\}\neq\emptyset$},$ (3)

with $m=d(i,j)$ representing the dyad between $i$ and $j$. Their model can be
embedded in our formulation by restricting (2) to be

$\displaystyle Cov(\varepsilon_{M,m},\varepsilon_{M,m^{\prime}}\mid
X_{M})=0\hskip 14.22636pt\text{unless $m$ and $m^{\prime}$ are adjacent}.$ (4)

####  2.1.2 Estimation

Throughout this paper, we focus on the Ordinary Least Squares (OLS) estimator
of $\beta$, denoted by $\hat{\beta}$. Under the assumption above, we can write

$\displaystyle\hat{\beta}-\beta$
$\displaystyle=\Bigg{(}\sum_{j\in\mathcal{M}_{N}}x_{M,j}x_{M,j}^{\prime}\Bigg{)}^{-1}\sum_{m\in\mathcal{M}_{N}}x_{M,m}\varepsilon_{M,m}.$
(5)

It is straightforward to verify that $\hat{\beta}$ is unbiased for $\beta$ due
to Assumption 2.2. However, a consistency result is by no means trivial due to
the dependence along the network which induces a complex form of cross-
sectional dependence, hindering a naïve application of the standard theory for
independently and identically distributed ($i.i.d.$) random vectors. This is
clarified in Section 2.2.

We now provide a heuristic overview of the main point of our paper in this
setting, before proving it formally in Section 4. In doing so, we omit the
subscript $M$ (as we did in $\hat{\beta}$) for notational convenience.

####  2.1.3 Inference

Inference about $\beta$ is based on a normal approximation of the distribution
of $\hat{\beta}$ around $\beta$. We focus on hypothesis testing conducted
using the expression:

$\displaystyle\big{(}\widehat{Var}(\hat{\beta})\big{)}^{-\frac{1}{2}}(\hat{\beta}-\beta),$
(6)

where $\widehat{Var}(\hat{\beta})$ is a consistent estimator of the asymptotic
variance of $\hat{\beta}$. Our main result in Section 4 is providing such an
appropriate estimator, which takes the form:

$\displaystyle\widehat{Var}(\hat{\beta})\coloneqq\Bigg{(}\sum_{k\in\mathcal{M}_{N}}x_{k}x_{k}^{\prime}\Bigg{)}^{-1}\Bigg{(}\sum_{m\in\mathcal{M}_{N}}\sum_{m^{\prime}\in\mathcal{M}_{N}}\kappa_{m,m^{\prime}}h_{m,m^{\prime}}\hat{\varepsilon}_{m}\hat{\varepsilon}_{m^{\prime}}x_{m}x_{m^{\prime}}^{\prime}\Bigg{)}\Bigg{(}\sum_{k\in\mathcal{M}_{N}}x_{k}x_{k}^{\prime}\Bigg{)}^{-1},$
(7)

where $\kappa_{m,m^{\prime}}$ is a kernel function that will formally be
defined in Section 4; $h_{m,m^{\prime}}$ represents an indicator function that
takes one if dyads $m$ and $m^{\prime}$ are connected and zero otherwise; and
$\hat{\varepsilon}_{m}\coloneqq y_{m}-x_{m}^{\prime}\hat{\beta}$. This paper
derives conditions under which $\widehat{Var}(\hat{\beta})$ is consistent for
the asymptotic variance of $\hat{\beta}$. Before doing so, let us compare the
variance estimator (7) with a dyadic-robust variance estimator often used in
the empirical literature, which uses one-way clustering based on dyad
groupings.

######  Example 2.2 (Dyadic-Robust Variance Estimator).

An increasing number of applied researchers, such as Harmon et al. (2019) and
Lustig and Richmond (2020), estimate model (1) under specification (2). To
conduct inference, these works employ dyadic-robust variance estimators
proposed by Aronow et al. (2015) and Tabord-Meehan (2019), which are given by:

$\displaystyle\widehat{Var}(\hat{\beta})\coloneqq\Bigg{(}\sum_{k\in\mathcal{M}_{N}}x_{k}x_{k}^{\prime}\Bigg{)}^{-1}\Bigg{(}\sum_{m\in\mathcal{M}_{N}}\sum_{m^{\prime}\in\mathcal{M}_{N}}\mathbbm{1}_{m,m^{\prime}}\hat{\varepsilon}_{m}\hat{\varepsilon}_{m^{\prime}}x_{m}x_{m^{\prime}}^{\prime}\Bigg{)}\Bigg{(}\sum_{k\in\mathcal{M}_{N}}x_{k}x_{k}^{\prime}\Bigg{)}^{-1},$
(8)

where $\mathbbm{1}_{m,m^{\prime}}$ equals one if dyads $m$ and $m^{\prime}$
are adjacent and zero otherwise. Note that the use of the dyadic-robust
variance estimator sets cases in which two dyads are not adjacent but
connected to zero. Note that our estimator (7) accounting for network
spillovers accommodates the correlation across both adjacent and connected
dyads,101010See Definition 2.1 and the subsequent discussion. The choice of
kernel and lag-truncation is discussed in Section 4. while the structure of
the variance estimator (8) may not be coherent with the model specification
(2). This suggests that the dyadic-robust variance estimator may be
inconsistent when there is cross-sectional dependence along a network over
dyads: i.e., when non-adjacent dyads can still affect the correlation
structure and outcomes of dyad $m$.111111Clustering estimators may be
inappropriate when the correlation structure has network spillovers as in (2),
because each agent has a complex (i.e., non-separable) structure of
connections, reflected in a non-separable network across dyads. If the network
model features positive spillovers, then the dyadic-robust variance estimator
will likely underestimate the true variance, leading to conservative
hypothesis testing. Meanwhile, it is likely to overstate the true variance
when there are negative spillovers. We expand on this point in our numerical
simulations. This conjecture is formally proven in Corollary 4.1 and
illustrated in Monte Carlo simulations in Section 5. We note that this is a
feature of applying such dyadic-robust variance estimators to network
spillovers, and not a feature of those estimators per se.

In order to validate the test statistic (6) for inference, this paper
establishes the consistency and asymptotic normality of $\hat{\beta}$ and
develops a consistent network-robust estimator for the asymptotic variance
when the correlation structure is given by (2).

###  2.2 Network Dependent Processes

Let $Y_{M,m}$ be a random vector defined as

$\displaystyle Y_{M,m}\coloneqq x_{M,m}\varepsilon_{M,m},$ (9)

and denote
$\mathcal{C}_{M}\coloneqq\\{x_{M,m}\\}_{m\in\mathcal{M}_{N}}$.131313For the
case of stochastic networks, it is defined to include information about the
network topology as well as the collection of the dyad-specific attributes
$\\{x_{M,m}\\}_{m\in\mathcal{M}_{N}}$. See Kojevnikov et al. (2021). From
equations (5) and (9) we can write

$\displaystyle\hat{\beta}-\beta=\Bigg{(}\frac{1}{M}\sum_{j\in\mathcal{M}_{N}}x_{M,j}x_{M,j}^{\prime}\Bigg{)}^{-1}\frac{1}{M}\sum_{m\in\mathcal{M}_{N}}Y_{M,m}.$
(10)

Our interest lies in proving the asymptotic properties of $\hat{\beta}$ taking
advantage of the expression (10). However, the presence of $\varepsilon_{M,m}$
in $Y_{M,m}$, which is allowed to be correlated along the network over active
dyads, renders our approach nonstandard. The canonical results about
independently and identically distributed ($i.i.d.$) random variables are by
no means applicable due to the complex form of dependence inherent to
$\varepsilon_{M,m}$. Dyadic-robust estimators (e.g., Aronow et al. (2015) and
Tabord-Meehan (2019)) only account for correlations between adjacent dyads.
This may be appropriate in some settings, but not when there may be indirect
spillovers, as in specification (1) – (2). Moreover, the asymptotic theory of
spatially-correlated and time-series data builds on the assumption that the
index set is a metrizable space, an assumption that may not necessarily be
true in the network context.141414Typically the literature assumes a Euclidean
space. See Section 1.1. In this light, the right hand side of (10) cannot
simply be embedded in other widely-studied variants of dependent random
vectors.

The main insight of this paper is that the spillovers across connected – even
if not adjacent – dyads, can be rewritten as the dependence of $Y_{M,m}$’s
along the network of active dyads (hereby, referred to as the “network”). This
transformation allows us to embrace such complex cross-sectional dependence
and appropriately rewrite the problem so that recent results on network
dependent random variables (Kojevnikov et al. (2021)) can be applied. Our
results then allow for correlations of random variables given a common shock
as explicitly dictated by distances between them on the dyadic network. This
dependence is then gauged by an appropriate notion of “network dependence” and
of distance across network dependent random variables, which is introduced
below. Compared to the other existing frameworks, the network-dependent random
variables circumvent the aforementioned limitations, so that the network-
regression model (1) and (2) become amenable to formal statistical inference.

To formally define our notion of network dependence, we first introduce
notation intentionally set close to that of Kojevnikov et al. (2021). By
construction, $Y_{M,m}$ resides in the $K$-dimensional Euclidean space
$\mathbbm{R}^{K}$. First, let $A$ and $B$ be sets of dyads whose cardinalities
are $a$ and $b$, respectively. We denote by $\mathcal{P}_{M}(a,b;s)$ the
collection of pairs $A$ and $B$ with distance $s$.151515The distance here is
meant to be the minimum geodesic distance between these two sets, which is
formally defined in Appendix A. Next, we consider a collection of bounded
Lipschitz functions $f:\mathbbm{R}^{K\times c}\rightarrow\mathbbm{R}$ and
denote it by $\mathcal{L}_{K,c}$. Let $\mathcal{L}_{K}$ be a collection that
gathers $\mathcal{L}_{K,c}$ for all $c\in\mathbbm{N}$.

Now we are finally in a position to formally define the notion of conditional
$\psi$-dependence, first introduced in Kojevnikov et al. (2021), as Definition
2.2.

######  Definition 2.2 (Conditional $\psi$-Dependence given
$\\{\mathcal{C}_{M}\\}$).

A triangular array $\\{Y_{M,m}\in\mathbbm{R}^{K}:M\geq
1,m\in\\{1,\ldots,M\\}\\}$ is called conditionally $\psi$-dependent given
$\\{\mathcal{C}_{M}\\}$, if for each $M\in\mathbbm{N}$, there exist a
$\mathcal{C}_{M}$-measurable sequence
$\theta_{M}\coloneqq\\{\theta_{M,s}\\}_{s\geq 0}$ with $\theta_{M,0}=1$, and a
collection of nonrandom function $(\psi_{a,b})_{a,b\in\mathbbm{N}}$ where
$\psi_{a,b}:\mathcal{L}_{K,a}\times\mathcal{L}_{K,b}\rightarrow[0,\infty)$,
such that for all $(A,B)\in\mathcal{P}_{M}(a,b;s)$ with $s>0$ and all
$f\in\mathcal{L}_{K,a}$ and $g\in\mathcal{L}_{K,b}$,

$\displaystyle\big{|}Cov(f(Y_{M,A}),g(Y_{M,B})\mid\mathcal{C}_{M})\big{|}\leq\psi_{a,b}(f,g)\theta_{M,s}\
\ \ \ \ a.s.$

We call $\\{\theta_{M,s}\\}$ the dependence coefficients of $\\{Y_{M,m}\\}$.
Intuitively, $\psi_{a,b}(f,g)$ dictates the magnitude of covariation of
$f(Y_{M,A})$ and $g(Y_{M,B})$. Meanwhile, $\theta_{M,s}$ captures the strength
of covariation net of the effect of scaling, analogous to correlation
coefficients. As the minimum distance $s$ between $A$ and $B$ grows, the
dependence $\\{\theta_{M,s}\\}$ between $Y_{M,A}$ and $Y_{M,B}$, tends to
zero.161616Note that $\\{\theta_{M,s}\\}$ are random because Definition 2.2 is
a conditional statement. We maintain Assumption 2.1 of Kojevnikov et al.
(2021) on the upper bounds of both $\psi_{a,b}(f,g)$ and $\theta_{M,s}$
throughout the paper.171717Details of the mathematical environment are
relegated to Appendix A.

We base our theory on Definition 2.2, which conditions on $\mathcal{C}_{M}$.
This allows us to focus only on the dependence of the random vectors $Y_{M,m}$
intrinsic to the error terms $\varepsilon_{M,m}$ (i.e., after conditioning out
those from covariates). This is of practical importance when $Y_{M,m}$ depends
on other covariates with the error terms being the sole source of the
systematic dependence structured over the network.

##  3 Preliminary Asymptotic Results

This section establishes the consistency and asymptotic normality of
$\hat{\beta}$ tailored for our network dependent random variables. This is
necessary before providing a consistent estimator for the variance of
$\hat{\beta}$. Let us first introduce some relevant notation.

As will become transparent shortly, asymptotic theories for $\hat{\beta}$ rest
on tradeoffs between the correlation of the network-dependent random vectors
(i.e., the dependence coefficients) and the denseness of the underlying
network. To measure the denseness, we first define two concepts of
neighborhoods: for each $m\in\mathcal{M}_{N}$ and
$s\in\mathbbm{N}\cup\\{0\\}$,

$\displaystyle\mathcal{M}_{N}(m;s)\coloneqq\\{m^{\prime}\in\mathcal{M}_{N}:\rho_{M}(m,m^{\prime})\leq
s\\},$
$\displaystyle\mathcal{M}_{N}^{\partial}(m;s)\coloneqq\\{m^{\prime}\in\mathcal{M}_{N}:\rho_{M}(m,m^{\prime})=s\\},$

where $\rho_{M}(m,m^{\prime})$ denotes the geodesic distance between dyads $m$
and $m^{\prime}$.181818Recall that we define the geodesic distance between two
connected dyads $m$ and $m^{\prime}$ to be the smallest number of adjacent
dyads between them. The former set collects all the $m$’s neighbors whose
distance from $m$ is no more than $s$ (which we call a neighborhood), whilst
the latter registers all the $m$’s neighbors whose distance from $m$ is
exactly $s$ (which we call a neighborhood shell).

Next, we define two types of density measures of a network: for $k,r>0$,

$\displaystyle\Delta_{M}(s,r;k)\coloneqq\frac{1}{M}\sum_{m\in\mathcal{M}_{N}}\max_{m^{\prime}\in\mathcal{M}_{N}^{\partial}(m;s)}|\mathcal{M}_{N}(m;r)\backslash\mathcal{M}_{N}(m^{\prime};s-1)|^{k},$
$\displaystyle\delta_{M}^{\partial}(s;k)\coloneqq\frac{1}{M}\sum_{m\in\mathcal{M}_{N}}|\mathcal{M}_{N}^{\partial}(m;s)|^{k},$

where it is assumed that $\mathcal{M}_{N}(m^{\prime};-1)=\emptyset$. The
former measure gauges the denseness of a network in terms of the average size
of a version of the neighborhood. The latter expresses the denseness of a
network as the average size of the neighborhood shell raised by $k$.
Kojevnikov et al. (2021) show that controlling the asymptotic behavior of an
appropriate composite of these two measures (denoted by $c_{M}$ and defined in
Assumption 4.1) is sufficient for the Law of Large Numbers (LLN) and Central
Limit Theorem (CLT) of the network dependent random variables (Condition ND).

###  3.1 Asymptotic Properties of $\hat{\beta}$

####  3.1.1 Consistency

To prove that $\hat{\beta}$ is consistent for $\beta$ requires two additional
conditions.

######  Assumption 3.1.

There exists $\eta>0$ such that $\sup_{N\geq
1}\max_{m\in\mathcal{M}_{N}}E\big{[}|\varepsilon_{M,m}|^{1+\eta}\mid\mathcal{C}_{M}\big{]}<\infty$.

######  Assumption 3.2.

$\frac{1}{M}\sum_{s\geq
1}\delta_{M}^{\partial}(s;1)\theta_{M,s}\stackrel{{\scriptstyle
a.s.}}{{\rightarrow}}0$ as $M\rightarrow\infty$.

Assumption 3.1 requires that the errors are not too large once conditioned on
common shocks. Together with the full rank assumption above, this implies
Assumption 3.1 of Kojevnikov et al. (2021) for each $u$-th element of
$Y_{M,m}$, denoted by $Y_{M,m}^{u}$ with $u\in\\{1,\ldots,K\\}$.

Assumption 3.2 controls the tradeoff between the denseness of the underlying
network and the covariability of the random vectors. If the network becomes
dense, then the dependence of the associated random variables has to decay
much faster. This requirement is consistent with the typical assumption in the
spatial econometrics literature upon which the empirical analyses of the
aforementioned examples (i.e., international trade, legislative voting, and
exchange rates) are based, because it embodies the idea that spillovers decay
as they propagate farther (see, e.g., Kelejian and Prucha (2010)). For
instance, Acemoglu et al. (2015) assume that network spillovers are zero if
agents are sufficiently distantly connected on a geographical network.

######  Theorem 3.1 (Consistency of $\hat{\beta}$).

Under Assumptions 2.1, 2.2, 3.1 and 3.2,
$\|\hat{\beta}-\beta\|_{2}\stackrel{{\scriptstyle p}}{{\rightarrow}}0$ as
$N\rightarrow\infty$.

###### Proof.

See Appendix B.2. ∎

Theorem 3.1 establishes the consistency of $\hat{\beta}$ under the scenario
where the number of sampling units $N$ goes to infinity. When Assumption 2.1
is dropped, the result continues to hold in terms of the number of active
dyads $M$.

####  3.1.2 Asymptotic Normality

Before we can prove the consistency of our variance estimator, we must derive
the asymptotic distribution of $\hat{\beta}$. To do so, we study the
asymptotic distribution of the sum of network dependent random vectors,
$S_{M}\coloneqq\sum_{m\in\mathcal{M}_{N}}Y_{M,m}$, which is present in
$\hat{\beta}$ in equation (10).

Let $S_{M}^{u}$ be the $u$-th entry of $S_{M}$ for $u\in\\{1,\ldots,K\\}$ and
denote the unconditional variance of $S_{M}^{u}$ by $\tau_{M}^{2}\coloneqq
Var(S_{M}^{u})$. Since $S_{M}^{u}$ is not a sum of independent variables, its
variance cannot be simply expressed as a sum of variances. We thus need to
explicitly take into account covariance between the random variables
$\\{Y_{M,m}^{u}\\}_{m\in\mathcal{M}_{N}}$. We study the CLT for the normalized
sum of $Y_{M,m}^{u}$, which is given by $\frac{S_{M}^{u}}{\tau_{M}}$.

As usual, the CLT for this normalized sum requires us to strengthen conditions
(3.1) and (3.2) to control the variability of higher moments of our random
variables. This is done in Assumptions A.2 and A.4 in Appendix A. These
assumptions are written in terms of conditional expectations, whereas we are
interested in the unconditional distribution of $\frac{S_{M}}{\tau_{M}}$.
Assumption A.3 in Appendix A then bridges the conditional and unconditional
variances by requiring that the conditional variance grows at the same rate as
the unconditional one, yielding a rate comparable to Lemma 1 in Aronow et al.
(2015).191919In Appendix B.4, we prove asymptotic normality under a milder
condition.

Finally, we require two regularity conditions to guarantee that the asymptotic
variance is well-defined. The first assumption, Assumption 3.3, is necessary
for one of the matrices in the expression to be well-defined. The uniform
boundedness of arbitrary moments of $x_{M,m}$ of Assumption 3.3 strengthens
Assumption 2.2 (b). The second assumption, Assumption 3.4, assures that the
middle part of the asymptotic variance is non-trivial.

######  Assumption 3.3.

(a) For all $N\geq 1$, $\\{x_{M,m}\\}_{m\in\mathcal{M}_{N}}$ have uniformly
bounded support.
(b)
$\lim_{N\rightarrow\infty}\Big{(}\frac{1}{M}\sum_{k\in\mathcal{M}}E\big{[}x_{M,k}x_{M,k}^{\prime}\big{]}\Big{)}$
is positive definite.

As pointed out in Tabord-Meehan (2019), the bounded support assumption can be
relaxed at the cost of imposing an alternative condition pertain to higher-
order moments. For instance, the boundedness of the 16th order moment of
$x_{M,m}$ is sufficient for the subsequent results. The second condition
guarantees the invertibility of an appropriate matrix for non-$i.i.d.$ random
variables in the asymptotic variance expression.

######  Assumption 3.4.

$\lim_{N\rightarrow\infty}\frac{N}{M^{2}}\sum_{m\in\mathcal{M}_{N}}\sum_{m^{\prime}\in\mathcal{M}_{N}}E\big{[}\varepsilon_{M,m}\varepsilon_{M,m^{\prime}}x_{M,m}x_{M,m^{\prime}}^{\prime}\big{]}$
exists with finite elements.

Under these assumptions, we obtain the asymptotic distribution of
$\hat{\beta}$:

######  Theorem 3.2 (Asymptotic Normality of $\hat{\beta}$).

Under Assumptions 2.1, 3.3, 3.4, A.2, A.3 and A.4,
$\sqrt{N}\big{(}\hat{\beta}-\beta\big{)}\stackrel{{\scriptstyle
d}}{{\rightarrow}}\mathcal{N}(0,AVar(\hat{\beta}))$ as $N\rightarrow\infty$,
where $AVar(\hat{\beta})$ is given by

$\displaystyle
AVar(\hat{\beta})\coloneqq\lim_{N\rightarrow\infty}N\Big{(}\sum_{k\in\mathcal{M}_{N}}E\big{[}x_{M,k}x_{M,k}^{\prime}\big{]}\Big{)}^{-1}\Big{(}\sum_{m\in\mathcal{M}_{N}}\sum_{m^{\prime}\in\mathcal{M}_{N}}E\big{[}\varepsilon_{M,m}\varepsilon_{M,m^{\prime}}x_{M,m}x_{M,m^{\prime}}^{\prime}\big{]}\Big{)}\Big{(}\sum_{k\in\mathcal{M}_{N}}E\big{[}x_{M,k}x_{M,k}^{\prime}\big{]}\Big{)}^{-1},$
(11)

which is positive semidefinite with finite elements.

###### Proof.

See Appendix B.4. ∎

Theorem 3.2 demonstrates the asymptotic normality of $\hat{\beta}$ in terms of
$N$, the number of sampling units.

##  4 Consistent Estimation of the Asymptotic Variance of $\hat{\beta}$ under
Network Spillovers

This section studies consistent estimation of $AVar(\hat{\beta})$ given that
the consistency and asymptotic normality of $\hat{\beta}$ are established
(Theorems 3.1 and 3.2). In light of Assumption 2.2 (iv), $Y_{M,m}$ is
centered, i.e., $E\big{[}Y_{M,m}\big{]}=0$ for each $m\in\mathcal{M}_{N}$.
Define $Var(\hat{\beta})$ as

$\displaystyle Var(\hat{\beta})$ $\displaystyle\coloneqq
N\Big{(}\frac{1}{M}\sum_{k\in\mathcal{M}_{N}}E\big{[}x_{M,k}x_{M,k}^{\prime}\big{]}\Big{)}^{-1}\Big{(}\frac{1}{M^{2}}\sum_{s\geq
0}\sum_{m\in\mathcal{M}_{N}}\sum_{m^{\prime}\in\mathcal{M}_{N}^{\partial}(m;s)}E\big{[}Y_{M,m}Y_{M,m^{\prime}}^{\prime}\big{]}\Big{)}\Big{(}\frac{1}{M}\sum_{k\in\mathcal{M}_{N}}E\big{[}x_{M,k}x_{M,k}^{\prime}\big{]}\Big{)}^{-1}.$
(12)

Clearly, $AVar(\hat{\beta})=\lim_{N\rightarrow\infty}Var(\hat{\beta})$. The
estimator developed in this section is a type of a kernel estimator. Define a
kernel function $\omega:\overline{\mathbbm{R}}\rightarrow[-1,1]$ such that
$\omega(0)=1$, $\omega(z)=0$ whenever $|z|>1$, and $\omega(z)=\omega(-z)$ for
all $z\in\overline{\mathbbm{R}}$. We denote by $b_{M}$ the bandwidth, or the
lag truncation. The feasible variance estimator of our interest (7) can be
written as

$\displaystyle\widehat{Var}(\hat{\beta})=\Big{(}\frac{1}{M}\sum_{k\in\mathcal{M}_{N}}x_{M,k}x_{M,k}^{\prime}\Big{)}^{-1}\Big{(}\frac{1}{M^{2}}\sum_{s\geq
0}\sum_{m\in\mathcal{M}_{N}}\sum_{m^{\prime}\in\mathcal{M}_{N}^{\partial}(m;s)}\omega_{M}(s)\hat{Y}_{M,m}\hat{Y}_{M,m^{\prime}}^{\prime}\Big{)}\Big{(}\frac{1}{M}\sum_{k\in\mathcal{M}_{N}}x_{M,k}x_{M,k}^{\prime}\Big{)}^{-1},$
(13)

with $\omega_{M}(s)\coloneqq\omega\big{(}\frac{s}{b_{M}}\big{)}$ for all
$s\geq 0$, and $\hat{Y}_{M,m}\coloneqq x_{M,m}\hat{\varepsilon}_{M,m}$, where
$\hat{\varepsilon}_{M,m}\coloneqq y_{M,m}-x_{M,m}^{\prime}\hat{\beta}$.

There are two approaches for the choice of the associated lag truncation
parameter. First, the researcher may already know (or is willing to impose)
the truncation, perhaps due to a theoretical/institutional motivation. For
instance, Acemoglu et al. (2015) set the lag to two in a related problem.
Alternatively, the researcher could use a data-driven choice. Assumption 4.1
(c) below suggests it should depend on both the sample size and the network
topology, including the average degree of the dyadic network. One such
selection rule is suggested in Kojevnikov et al. (2021) based on their proofs:
$b_{M}=2\log(M)/\log(\max(average~{}degree,1.05))$.

The consistency of the variance estimator requires two sets of additional
assumptions. The first set is Assumption 4.1 of Kojevnikov et al. (2021), but
stated here in terms of the network over dyads.

######  Assumption 4.1 (Kojevnikov et al. (2021), Assumption 4.1).

There exists $p>4$ such that (a) $\sup_{N\geq
1}\max_{m\in\mathcal{M}_{N}}E\big{[}|\varepsilon_{m}|^{p}\mid\mathcal{C}_{M}\big{]}<\infty$
$a.s.$;
(b) $\lim_{M\rightarrow\infty}\sum_{s\geq
1}|\omega_{M}(s)-1|\delta_{M}^{\partial}(s)\theta_{M,s}^{1-\frac{2}{p}}=0$
$a.s.$;
(c) $\lim_{M\rightarrow\infty}\frac{1}{M}\sum_{s\geq
0}c_{M}(s,b_{M};2)\theta_{M,s}^{1-\frac{4}{p}}=0$ $a.s.$, where

$\displaystyle
c_{M}(s,m;k)\coloneqq\inf_{\alpha>1}\big{(}\Delta_{M}(s,r;k\alpha)\big{)}^{\frac{1}{\alpha}}\Big{(}\delta_{M}^{\partial}\big{(}s;\frac{\alpha}{\alpha-1}\big{)}\Big{)}^{\frac{\alpha-1}{\alpha}}.$

Assumption 4.1 (a) is a stronger counterpart to Assumption A.2, as it requires
that a higher-order (i.e., higher than fourth order) conditional moment be
well-defined. Assumption (b) posits a tradeoff between the kernel function,
the denseness of a network and the dependence coefficients. Specifically, the
kernel function $\omega_{M}$ is required to converge to one sufficiently fast.
Kojevnikov et al. (2021) demonstrate primitive conditions under which this
requirement is fulfilled (Proposition 4.2). Assumption (c) indicates a
condition between the bandwidth, the denseness and the dependence
coefficients, which is preserved in the suggested choice for $b_{M}$ above.
For a fixed $p>4$, (b) and (c) preserve the same interpretation as Assumption
3.2: namely, the correlation must decay much faster relative to the growth of
the density of the network, or vice versa.

Another set of conditions refers to the denseness of the network.

######  Assumption 4.2.

(a) $\sup_{N\geq 1}\sum_{s\geq 0}\delta_{M}^{\partial}(s;1)<\infty$; (b)
$\lim_{M\rightarrow\infty}\frac{1}{M}\sum_{s\geq 0}c_{M}(s,b_{M};2)=0$.

Both requirements of Assumption 4.2 restrict the denseness of the network
among dyads: the requirement (a) means that the total density measured by the
neighborhood shell is bounded; and the requirement (b) embodies the idea that
the average density, evaluated in terms of the composite measure of the
denseness, vanishes as $M$ goes to infinity (see Assumption 4.1 (c) for the
definition of $c_{M}$). These conditions rule out the situation where the
network becomes progressively dense: most notably, the case where every single
individual unit is directly linked to every other individual. Note that, while
Assumption 4.1 constrains the denseness of the network relative to the
correlations of the associated random variables, Assumption 4.2 only refers to
the network configuration.

The following theorem is the main theoretical contribution of this paper.

######  Theorem 4.1 (Consistency of the Network-Robust Variance Estimator).

Under the conditions for Theorem 3.2, and Assumptions 4.1 and 4.2,
$\big{\|}N\widehat{Var}(\hat{\beta})-Var(\hat{\beta})\big{\|}_{F}\stackrel{{\scriptstyle
p}}{{\rightarrow}}0$ as $N\rightarrow\infty$, where $\|\cdot\|_{F}$ indicates
the Frobenius norm.

###### Proof.

See Appendix C.2. ∎

Theorem 4.1 establishes the consistency of our proposed variance estimator
accounting for network spillovers across dyads in the sense of the Frobenius
norm.

It follows from Theorem 4.1 that the dyadic-robust variance estimator (8) is
inconsistent for the true variance when the underlying network involves a non-
negligible degree of far-away correlations.202020If the network is such that
there are only adjacent dyads (i.e., when equation (4) is true), then the
result above implies consistency of this estimator for dyadic dependence. By
comparison, Lemma 1 of Aronow et al. (2015) and Proposition 3.1 and 3.2 of
Tabord-Meehan (2019) also provide consistent variance estimators for the
dyadic dependence case without higher-order network spillovers. However, these
results and ours do not subsume one another. Indeed, their estimators can
accommodate flexible dependence within clusters of dyads that share common
units, while we assume that the network of spillovers is observed even if
there are only adjacent connections. Specifically, the following corollary
states that the dyadic-robust variance estimator of Aronow et al. (2015) may
not necessarily be consistent when it is naïvely applied to the network-
regression model with non-zero correlations beyond direct neighbors.

######  Corollary 4.1 (Inconsistency of Dyadic-Robust Estimators with Network
Spillovers).

Suppose that the assumptions required in Theorem 4.1 hold. Assume, in
addition, that

$\displaystyle\inf_{N\geq 1}\frac{1}{M}\sum_{s\geq
2}\sum_{m\in\mathcal{M}_{N}}\sum_{m^{\prime}\in\mathcal{M}_{N}^{\partial}(m;s)}\big{\|}E\big{[}\varepsilon_{M,m}\varepsilon_{M,m^{\prime}}x_{M,m}x_{M,m^{\prime}}^{\prime}\big{]}\big{\|}_{F}>0.$
(14)

Then, the dyadic-robust estimator (8) applied to the network-regression model
(1) and (2) is inconsistent.

###### Proof.

See Appendix C.3. ∎

The added assumption (14) in Corollary 4.1 pertains to both the network
configuration of active dyads and the regression variables. It represents a
setting where the spillovers from far-away neighbors are not negligible even
when $N$ is large. For instance, the condition (14) can still hold even if
there are not so many neighbors, as long as the covariances between the error
terms are sufficiently large. On the other hand, if far-away neighbors in the
network have only a negligible effect on the cross-sectional dependence, the
dyadic-robust variance estimator remains a good approximation for the
asymptotic variance of linear dyadic data models with network spillovers
across dyads. These insights are investigated further in Section 5 using
numerical simulations.

##  5 Monte Carlo Simulation

###  5.1 Simulation Design

We compare three types of variance estimators across different specifications
and network configurations. We use the Eicker-Huber-White estimator as a
benchmark, the dyadic-robust estimator of Tabord-Meehan (2019) as a comparison
accounting for the dyadic nature of the data (when inappropriately used in the
presence of network spillovers), and our proposed estimator which is robust to
network spillovers across dyads.

We first generate networks on which random variables are assigned. We follow
Canen et al. (2020) among others by employing two models of random graph
formations. They are referred to as Specification 1 and 2. Specification 1
uses the Barabási and Albert’s (1999) model of preferential attachment, with
the fixed number of edges $\nu\in\\{1,2,3\\}$ being established by each new
node.212121In generating the Barabási-Albert random graphs, we follow Canen et
al. (2020) by choosing the seed to be the Erdös-Renyi random graph with the
number of nodes equal the smallest integer above $5\sqrt{N}$, where $N$
denotes the number of nodes. Specification 2 is based on the Erdös-Renyi
random graph (Erdös and Rényi, 1959, 1960) with probability
$p=\frac{\lambda}{N}$ for $N$ denoting the number of nodes and
$\lambda\in\\{1,2,3\\}$ being a parameter that governs the probability
relative to the node size. The summary statistics for the networks generated
by Specification 1 and 2 are given in Appendix D.1. By construction, the
maximum degree and the average degree increase monotonically as we increase
the parameters in both specifications. The number of active edges (i.e.,
dyads) also increases with the sample size regardless of the specification.
This reflects the fact that each node tends to have more far-away neighbors as
the network becomes denser.

For each of the randomly generated networks, the simulation data is generated
from the following simple network-linear regression:

$\displaystyle y_{m}=x_{m}\beta+\varepsilon_{m},$

with $m\coloneqq d(i,j)$ representing the dyad between agent $i$ and $j$. The
dyad-specific regressor $x_{m}$ is defined as $x_{m}\coloneqq|z_{i}-z_{j}|$,
where both $z_{i}$ and $z_{j}$ are drawn independently from
$\mathcal{N}(0,1)$. The regression coefficient is fixed to $\beta=1$ across
specifications.

The dyad-specific error term $\varepsilon_{m}$ is constructed to exhibit non-
zero correlation with $\varepsilon_{m^{\prime}}$ as long as dyads $m$ and
$m^{\prime}$ are connected (i.e., in the network terminology, there exists a
path in the simulated network), while the strength of the correlation is
assumed to decay as they are more distant. To that end, we draw
$\varepsilon_{m}\coloneqq\sum_{m^{\prime}}\gamma_{m,m^{\prime}}\eta_{m,m^{\prime}}$,
where $\gamma_{m,m^{\prime}}$ equals $\gamma^{s}$ if the distance between $m$
and $m^{\prime}$ is $s$, and $0$ otherwise, for $\gamma\in[0,1]$232323In this
simulation, we focus on cases of positive spillovers. Cases of negative
spillovers can be analysed in an analogous manner. and $s\in\\{1,\ldots,S\\}$
with $S$ being the maximum geodesic distance that the spillover propagates to.
Each $\eta_{m,m^{\prime}}$ is drawn $i.i.d.$ from $\mathcal{N}(0,1)$. If
$\gamma=1$, then spillover effects are the same no matter how far the agents
are apart, i.e., the spillover effects do not decay. If $\gamma=0$, there are
no spillover effects, so the dyadic-robust variance estimator should be
consistent. The case of $S=2$ corresponds to a situation where up to friends
of friends may matter for spillovers.

We consider three scenarios for each type of network. In the main text, we set
$S=2$ and $\gamma=0.8$. The results for $S=2$ with $\gamma=0.2$ are given in
Appendix D.3, and the ones for $S=1$ with $\gamma=0.8$ are in Appendix D.4.
Throughout the experiments, we employ the rectangular kernel with the lag
truncation equal two for the sake of comparison.242424This choice of the
kernel function and the lag truncation corresponds to placing unit weights
both on pairs of adjacent dyads and pairs of connected dyads of distance equal
two, and zero weights otherwise. As far as expositional purpose is concerned,
this is sufficient for us to accommodate all the correlations across the
network since it is constructed to be $S=2$. This set-up is equivalent to
assuming that the researcher has correct information regarding how far the
network spillovers propagate to. In practice, however, it is rare to be given
such knowledge prior to analysis, and one can use the choice of lag truncation
suggested in the previous section. The nominal size is set at $\alpha=0.05$
throughout the $5000$ iterations of the Monte Carlo simulation.

###  5.2 Results

In Table 1 we present the coverage probability for $\beta$ and the average
length of the confidence interval across simulations. To do so, we compute the
$t$-statistic using the OLS estimator for $\beta$ and different variance
estimators under a Normal distribution approximation.252525It is well known
that estimates of a variance-covariance matrix may be negative semidefinite
when the sample size is very small. This occurs in four out of 5000
simulations when $N=500$. Rather than dropping such observations, we follow
Cameron et al. (2011) and augment the eigenvalues of the matrix by adding a
small constant, say $0.005$, thereby obtaining a new variance estimate that is
more conservative. See the note below Table D.4 for detail of the procedure.
The finite-sample properties of the three variance estimators are further
illustrated in Figure D.2 in Appendix D.2.

The results for the empirical coverage probabilities depend on two dimensions:
the sample size ($N$) and the denseness of the underlying network
(parametrized by $\nu$ and $\lambda$). The coverage probability for each
estimator improves with the sample size. However, when spillovers are high
($\gamma=0.8$), only our proposed network-robust variance estimator has
coverage close to 95%, as predicted in Theorem 4.1. Meanwhile, both the
Eicker-Huber-White and the dyadic-robust variance estimators perform poorly as
the underlying network becomes denser, no matter which specification of the
network is involved. For example, in Specification 1 with $\nu=3$ and the
largest sample size ($N=5000$), the confidence intervals based on the Eicker-
Huber-White and the dyadic-robust variance estimators do not cover the true
parameter 615 and 455 times out of 5000 simulations (12.3% and 9.1%),
respectively. This is because as the network gets denser, each agent tends to
have more direct or indirect links to other agents, leading to larger higher-
order correlations, which is neither accounted by the Eicker-Huber-White nor
by the dyadic-robust variance estimator. On the other hand, the network-robust
variance estimator is designed to capture such higher-order correlations and,
thus, its coverage remains stable across network configurations.

A similar conclusion is drawn from the average length of the confidence
intervals. The confidence intervals for the Eicker-Huber-White and dyadic-
robust variance estimators are typically shorter than those for our proposed
estimator when $\gamma$ is large and $S=2$. This is because the former
undercovers the true parameter (in the presence of positive spillovers),
confirming our theoretical results (Corollary 4.1). The difference is often
around 10-20% of the length of the network-robust estimator.

However, as the magnitude of spillovers decreases (i.e. $\gamma$ tends to 0),
higher-order spillovers are less pronounced, so that the biases from using the
Eicker-Huber-White and dyadic-robust variance estimators disappear. This is
shown in Table D.4 of Appendix D.3 for the case of $S=2$ and $\gamma=0.2$.
When $S=1$, the dyadic-robust variance estimator coincides with our proposed
estimator (i.e., there are no spillovers from non-adjacent links). This is
shown in Table D.5 of Appendix D.4.

Table 1: The empirical coverage probability and average length of confidence
intervals for $\beta$ at $95\%$ nominal level: $S=2$, $\gamma=0.8$.

|  | Specification 1 |  | Specification 2
---|---|---|---|---
| $N$ | $\nu=1$ | $\nu=2$ | $\nu=3$ |  | $\lambda=1$ | $\lambda=2$ | $\lambda=3$
|  | Coverage Probability |
Eicker-Huber-White | 500 | 0.8774 | 0.8680 | 0.8712 |  | 0.8910 | 0.8700 | 0.8750
| 1000 | 0.8798 | 0.8732 | 0.8728 |  | 0.8924 | 0.8810 | 0.8876
| 5000 | 0.8792 | 0.8708 | 0.8770 |  | 0.8928 | 0.8890 | 0.8790
Dyadic-robust | 500 | 0.9216 | 0.8976 | 0.8938 |  | 0.9320 | 0.9212 | 0.9172
| 1000 | 0.9286 | 0.9134 | 0.9006 |  | 0.9366 | 0.9274 | 0.9242
| 5000 | 0.9342 | 0.9118 | 0.9090 |  | 0.9388 | 0.9332 | 0.9224
Network-robust | 500 | 0.9300 | 0.9166 | 0.9146 |  | 0.9370 | 0.9370 | 0.9410
| 1000 | 0.9392 | 0.9338 | 0.9328 |  | 0.9456 | 0.9446 | 0.9476
| 5000 | 0.9486 | 0.9440 | 0.9426 |  | 0.9468 | 0.9472 | 0.9480
|  | Average Length of the Confidence Intervals |
Eicker-Huber-White | 500 | 0.3678 | 0.4093 | 0.4819 |  | 0.2874 | 0.2853 | 0.2962
| 1000 | 0.2663 | 0.3024 | 0.3309 |  | 0.2048 | 0.2013 | 0.2068
| 5000 | 0.1321 | 0.1585 | 0.1756 |  | 0.0916 | 0.0898 | 0.0938
Dyadic-robust | 500 | 0.4260 | 0.4543 | 0.5204 |  | 0.3280 | 0.3285 | 0.3365
| 1000 | 0.3124 | 0.3393 | 0.3611 |  | 0.2359 | 0.2323 | 0.2366
| 5000 | 0.1580 | 0.1775 | 0.1920 |  | 0.1060 | 0.1043 | 0.1077
Network-robust | 500 | 0.4413 | 0.4932 | 0.5681 |  | 0.3370 | 0.3487 | 0.3655
| 1000 | 0.3263 | 0.3733 | 0.4079 |  | 0.2435 | 0.2482 | 0.2585
| 5000 | 0.1666 | 0.1994 | 0.2215 |  | 0.1097 | 0.1117 | 0.1183

  * Note: The upper-half of the table displays the empirical coverage probability of the asymptotic confidence interval for $\beta$, and the lower-half showcases the average length of the estimated confidence intervals. As the sample size ($N$) increases, the empirical coverage probability for our estimator accounting for network spillovers approaches $0.95$, the correct nominal level. However, that is not the case for alternative estimators.

##  6 Empirical Illustration: Legislative Voting in the European Parliament

We now turn to an empirical application to demonstrate the performance of our
variance estimator with real-world data. In doing so, we revisit the work of
Harmon et al. (2019) on whether legislators who sit next to each other in
Parliament tend to vote more alike on policy proposals.

In their original article, the authors study the existence of peer effects in
legislative voting. They focus on the European Parliament, whose Members
(MEPs) are voted in through elections in each European Union (EU) member
country every five years. The Parliament convenes once or twice a month, in
either Brussels or Strasbourg, to debate and vote on a series of proposals.

Once elected to the European Parliament (EP), these MEPs are organized into
European Political Groups (EPGs), which aggregate similar ideological
members/parties across countries. As Harmon et al. (2019) describe, these EPGs
function as parties for many of the traditional party-functions in other
legislatures, including coordination on policy and policy votes. Most
importantly, MEPs sit within their EPG groups in the chamber. However, within
each EPG group, non-party leaders traditionally sit in alphabetical order by
last name. See Figure E.4 in Appendix E.1 for an example. Hence, to the extent
that last names are exogenous (i.e., politicians are not selected or choose
last names anticipating seating arrangements in the European Parliament),
seating arrangements are used as a quasi-natural experiment: it is sometimes
by chance that one legislator sits next to another. And across different
periods, as politicians get replaced, they find themselves sitting with new
neighbors.

###  6.1 Data

We adopt the same dataset used in Harmon et al. (2019), which collects the
MEP-level data on votes cast in the EP. The dataset records what each MEP
voted for (Yes or No), where she was seated, and a number of individual
characteristics (e.g. country, age, education, gender, tenure, etc). Its
sample period is the plenary sessions held in both Brussels and Strasbourg
between October 2006 and November 2010. This corresponds to the sixth and
seventh terms of the European Parliament.

For our empirical illustration, we restrict the sample to the policies voted
in Strasbourg during the seventh term and we focus on the seating pattern
between July 14th - July 16th, 2009 (which involved 116 different proposals
being voted on). The resulting sample has 2,431,261 observations, which are
split into 422 politicians forming 26,099 pairs (i.e., dyads) of MEPs over 116
proposals.262626There are 334 pairs of adjacent dyads and 591 pairs of
connected dyads. Further information on the construction of our sample is
detailed in Appendix E.2. These restrictions keep the main set-up in Harmon et
al. (2019), while allowing us to evaluate variance estimators with smaller
sample sizes.

###  6.2 Empirical Setup

The empirical analysis in Harmon et al. (2019) studies whether MEPs who are
seated together are more likely to vote similarly. To make it consistent with
our notation above, let $m=d(i,j)$ denote the dyad between legislator $i$ and
$j$, with $t$ being used to index the proposal. Let $Agree_{d(i,j),t}$ be an
indicator that takes one if MEP $i$ and $j$ cast the same vote on proposal
$t$, and zero otherwise. Likewise, $SeatNeighbors_{d(i,j),t}$ is a binary
variable that equals one if MEP $i$ and $j$ are seated next to each other when
the vote for proposal $t$ is taken place, and zero otherwise. We view the
seating arrangement as a network among the MEPs: i.e., two MEPs who are seated
next to each other within the same political group are treated as an active
dyad, which in turn constitutes a network over pairs of MEPs within each
political group (see the note below Figure E.4 for details).272727This
accommodates the row-by-EP-by-EPG clustering implemented by the authors. This
is already an extension to Harmon et al. (2019) which assume that only seating
neighbors would have correlated error termsß. One could extend this further by
allowing connections across parties.

We follow Harmon et al. (2019) in assuming that such seating arrangements are
exogenously determined, i.e., the adjacency relation is exogenously formed.
Their main specification is a linear model:

$\displaystyle
Agree_{d(i,j),t}=\beta_{0}+\beta_{1}SeatNeighbors_{d(i,j),t}+\varepsilon_{d(i,j),t}.$
(15)

The authors originally conducted inference using the estimator in Aronow et
al. (2015), assuming that dyads $m=d(i,j)$ cannot be correlated with
$m^{\prime}=d(k,l)$ unless they share a common member (i.e., there is no
correlation across errors if $i,j$ and $k,l$ do not include a common unit). We
will now compare this approach to using the variance estimator introduced in
Section 4, which allows the error terms to exhibit non-zero correlations as
long as they are connected on the network over dyads represented by the
adjacency relation of seating arrangements in Parliament. In this analysis, we
use the rectangular kernel with the lag truncation equal the longest path in
the constructed network.282828This combination of the kernel density and the
lag truncation enables us to accommodate all the possible correlations across
connected dyads (i.e., pairs of MEPs).

Inspired by Harmon et al. (2019), we consider three specifications: (I) a
simple linear regression model as given in (15); (II) the model (15) augmented
with a flexible set of other demographic variables;292929Following Harmon et
al. (2019), we include indicators whether country of origins, quality of
education, freshman status and gender, respectively, are the same, as well as
differences in ages and tenures. See the note below Table 2 for details. and
(III) the model (15) with both a flexible set of other demographic variables
and day-specific fixed effects. When fixed effects are present in their
original estimation, we estimate a within-difference model via OLS.

###  6.3 Results

The main results of our empirical analysis are summarized in Table 2. Panel A
displays the parameter estimates for the three different specifications. It
can be seen in this panel that our point-estimates are consistent in magnitude
with the original estimates of Harmon et al. (2019) (columns 6 and 7 of Table
4), as they are close to 0.006 (their original results) and stable across
specifications.303030Note that our dependent variable is equal to one if two
MEPs vote the same and zero otherwise, while Harmon et al. (2019) code it as
one if MEPs vote differently. Hence, to compare our estimates with theirs, the
signs on the estimates of $SeatNeighbors$ must be flipped. Hence, changes to
point-estimates are not due to sample selection. The positive coefficient for
SeatNeighbors indicates that the MEPs sitting together tend to vote more
similarly than those sitting apart, providing evidence in favor of their
original hypothesis. The coefficients on the covariates (displayed in Panel C
of Appendix Table E.7) are also quantitatively and qualitatively similar to
those in their original paper. For instance, our estimates for SameCountry are
0.056, while their estimates are around 0.051, suggesting politicians from the
same country are more likely to vote similarly on policies.

Panel B shows the standard errors for the regression coefficient of
SeatNeighbors using different variance estimators. Building on Section 5, the
panel compares three different types of heteroskedasticity robust estimators:
namely, the Eicker-Huber-White, the dyadic-robust estimator (used in their
original work), and our proposed network-robust estimator. We note that, due
to the smaller sample size, the standard errors in our exercise are larger
than the original authors’. Hence, the results in this section are not
directly comparable to those in Harmon et al. (2019). Rather, we compare the
different estimators within our own sample.

As foreshadowed in the Monte Carlo simulations, the Eicker-Huber-White
estimates are the smallest, followed by the dyadic-robust estimates, which, in
turn, are smaller than the network-robust estimates. In fact, for
Specification (III), the Eicker-Huber-White estimate is roughly 73% smaller
than using the estimator accounting for network spillovers across dyads, while
the dyadic-robust one is 22% smaller. This fact entails two implications.
First, our finding provides empirical evidence in support of the existence of
indirect positive spillovers among the MEPs: even distant connections may
indirectly generate correlated behavior among politicians $i$ and $j$. Second,
the use of alternative estimators not accounting for such spillovers
undercovers the true parameter and may generate biased hypothesis testing
about the regression coefficient of SeatNeighbors. The difference in estimates
appears quantitatively meaningful in this empirical example.

Table 2: Spillovers in Legislative Voting – Main Analysis

| Specification (I) | Specification (II) | Specification (III)
---|---|---|---
Panel A: Parameter estimates for Seat neighbors |  |  |
Seat neighbors | 0.0069 | 0.0060 | 0.0060
Panel B: Standard errors for Seat neighbors |  |  |
Eicker-Huber-White | 0.0031 | 0.0030 | 0.0030
Dyadic-robust | 0.0075 | 0.0082 | 0.0087
Network-robust | 0.0095 | 0.0104 | 0.0112

  * Note: Panel A displays the parameter estimates for the three different specifications, and Panel B shows the standard errors for the regression coefficient of SeatNeighbors using different variance estimators. Adjacency of MEPs is defined at the level of a row-by-EP-by-EPG. (See the note below Figure E.4.) Independent variables are as follows: Seat neighbors is an indicator variable denoting whether both MEPs sit together; Same country represents an indicator for whether both MEPs are from the same country; Same quality education is an indicator showing whether both MEPs have the same quality of education background, measured by if both have the degree from top 500 universities; Same freshman status encodes whether both MEPs are freshman or not; Age difference is the difference in the MEPs’ ages; and Tenure difference measures the difference in the MEPs’ tenures. A full description of the result is provided in Table E.7.

##  7 Conclusion

Dyadic data have widespread applications in social sciences. Researchers
typically assume that dyads are uncorrelated if they do not share a common
unit, an assumption that is leveraged in inference. We showed that this is not
the case in many models whereby strategic interactions occur on a network:
while data may be dyadic, the cross-sectional dependence may be much more
complex and spillover beyond pairwise interactions. For instance, trade
between countries may depend on trade between auxiliary partners, even those
beyond the dyads being considered. In political economy, whether a politician
votes with a colleague may depend on all intermediate seating neighbors beyond
one’s own immediate ones. We verified this using both theoretical results and
Monte Carlo simulations. A new consistent variance estimator for parameters in
a linear model with dyadic data was derived, which has correct asymptotic
coverage and good finite sample properties. This estimator can account for
network spillovers while existing ones (e.g., Aronow et al. (2015) and Tabord-
Meehan (2019)) cannot, although it requires different assumptions. It is also
simple to implement and provided stable inference in an application to
concordance in legislative voting across pairs of members using real-world
data. We found that, when accounting for network spillovers across dyads,
standard errors grew by at least 20% compared to the original specification.

To conclude, we would like to clarify that our goal in this exercise is
neither to criticize dyadic-robust variance estimators, which are a
fundamental part of the empiricist’s toolkit, nor to suggest our approach
should always be used. Rather, we wish to draw attention that researchers
should fully specify the cross-sectional dependence in their model. If the
conventional assumption of dyadic dependence correctly specifies the
environment in question, or when spillovers beyond immediate neighbors might
be negligible, then previous approaches suffice. However, as we have discussed
above, many existing applications seemingly apply the latter method even if it
is seemingly inappropriate to their setting. Hence, we recommend researchers
to continue to fully specify their model, including full specification of
their covariance structure, thereby clarifying what type of inference
procedure is most appropriate for their environment.

## References

  * Abadie et al. (2022) Abadie, A., S. Athey, G. W. Imbens, and J. M. Wooldridge (2022).  When should you adjust standard errors for clustering?  The Quarterly Journal of Economics Forthcoming.
  * Acemoglu et al. (2015) Acemoglu, D., C. García-Jimeno, and J. A. Robinson (2015).  State capacity and economic development: A network approach.  American Economic Review 105(8), 2364–2409.
  * Anderson and van Wincoop (2003) Anderson, J. E. and E. van Wincoop (2003).  Gravity with gravitas: A solution to the border puzzle.  American Economic Review 93(1), 170–192.
  * Aronow et al. (2015) Aronow, P. M., C. Samii, and V. A. Assenova (2015).  Cluster–robust variance estimation for dyadic data.  Political Analysis 23(4), 564–577.
  * Banerjee et al. (2013) Banerjee, A., A. G. Chandrasekhar, E. Duflo, and M. O. Jackson (2013).  The diffusion of microfinance.  Science 341(6144).
  * Barabási and Albert (1999) Barabási, A. L. and R. Albert (1999).  Emergence of scaling in random networks.  Science 286(5439), 509–512.
  * Cameron et al. (2011) Cameron, A. C., J. B. Gelbach, and D. L. Miller (2011).  Robust inference with multiway clustering.  Journal of Business & Economic Statistics 29(2), 238–249.
  * Cameron and Miller (2014) Cameron, A. C. and D. L. Miller (2014).  Robust inference for dyadic data.  Working Paper.
  * Canen et al. (2020) Canen, N., J. Schwartz, and K. Song (2020).  Estimating local interactions among many agents who observe their neighbors.  Quantitative Economics 11(1), 917–956.
  * Chiang et al. (2022) Chiang, H. D., B. E. Hansen, and Y. Sasaki (2022).  Standard errors for two-way clustering with serially correlated time effects.  Working Paper.
  * Chiang et al. (2021) Chiang, H. D., K. Kato, Y. Ma, and Y. Sasaki (2021).  Multiway cluster robust double/debiased machine learning.  Journal of Business & Economic Statistics, 1–11.
  * Chiang et al. (2022) Chiang, H. D., Y. Matsushita, and T. Otsu (2022).  Multiway empirical likelihood.  Working Paper.
  * Chiang and Tan (2022) Chiang, H. D. and B. Y. Tan (2022).  Empirical likelihood and uniform convergence rates for dyadic kernel density estimation.  Journal of Business & Economic Statistics, 1–9.
  * Conley (1999) Conley, T. G. (1999).  GMM estimation with cross sectional dependence.  Journal of Econometrics 92(1), 1–45.
  * Davezies et al. (2021) Davezies, L., X. D’Haultfœuille, and Y. Guyonvarch (2021).  Empirical process results for exchangeable arrays.  The Annals of Statistics 49(2), 845–862.
  * de Paula et al. (2018) de Paula, A., S. Richards-Shubik, and E. Tamer (2018).  Identifying preferences in networks with bounded degree.  Econometrica 86(1), 263–288.
  * Erdös and Rényi (1959) Erdös, P. and A. Rényi (1959).  On random graphs i.  Publicationes Mathematicae 6, 290–297.
  * Erdös and Rényi (1960) Erdös, P. and A. Rényi (1960).  On the evolution of random graphs.  Publications of the Mathematical Institute of the Hungarian Academy of Sciences 5, 17–61.
  * Fafchamps and Gubert (2007) Fafchamps, M. and F. Gubert (2007).  The formation of risk sharing networks.  Journal of Development Economics 83(2), 326 – 350.
  * Graham (2020a) Graham, B. S. (2020a).  Dyadic regression.  In B. Graham and Á. d. Paula (Eds.), The Econometric Analysis of Network Data (1st ed.)., Chapter 2, pp. 23–40. Academic Press.
  * Graham (2020b) Graham, B. S. (2020b).  Network data.  In S. N. Durlauf, L. P. Hansen, J. J. Heckman, and R. L. Matzkin (Eds.), Handbook of Econometrics, Volume 7A, Volume 7 of Handbook of Econometrics, Chapter 2, pp. 111–218. Elsevier.
  * Graham et al. (2022) Graham, B. S., F. Niu, and J. L. Powell (2022).  Kernel density estimation for undirected dyadic data.  Working Paper.
  * Hansen and Lee (2019) Hansen, B. E. and S. Lee (2019).  Asymptotic theory for clustered samples.  Journal of Econometrics 210(2), 268–290.
  * Harmon et al. (2019) Harmon, N., R. Fisman, and E. Kamenica (2019).  Peer effects in legislative voting.  American Economic Journal: Applied Economics 11(4), 156–80.
  * Ibragimov and Müller (2010) Ibragimov, R. and U. K. Müller (2010).  t-statistic based correlation and heterogeneity robust inference.  Journal of Business & Economic Statistics 28(4), 453–468.
  * Jenish and Prucha (2009) Jenish, N. and I. R. Prucha (2009).  Central limit theorems and uniform laws of large numbers for arrays of random fields.  Journal of Econometrics 150(1), 86–98.
  * Kelejian and Prucha (2010) Kelejian, H. H. and I. R. Prucha (2010).  Specification and estimation of spatial autoregressive models with autoregressive and heteroskedastic disturbances.  Journal of Econometrics 157(1), 53–67.
  * Kojevnikov et al. (2021) Kojevnikov, D., V. Marmer, and K. Song (2021).  Limit theorems for network dependent random variables.  Journal of Econometrics 222(2), 882–908.
  * Leung (2021) Leung, M. P. (2021).  Network cluster-robust inference.  Working Paper.
  * Leung (2022) Leung, M. P. (2022).  Causal inference under approximate neighborhood interference.  Econometrica 90(1), 267–293.
  * Leung and Moon (2021) Leung, M. P. and H. R. Moon (2021).  Normal approximation in large network models.  Working Paper.
  * Lustig and Richmond (2020) Lustig, H. and R. J. Richmond (2020).  Gravity in the exchange rate factor structure.  The Review of Financial Studies 33(8), 3492–3540.
  * MacKinnon et al. (2022a) MacKinnon, J. G., M. Ø. Nielsen, and M. D. Webb (2022a).  Cluster-robust inference: A guide to empirical practice.  Journal of Econometrics Forthcoming.
  * MacKinnon et al. (2022b) MacKinnon, J. G., M. Ø. Nielsen, and M. D. Webb (2022b).  Fast and reliable jackknife and bootstrap methods for cluster-robust inference.  Working Paper.
  * MacKinnon et al. (2022c) MacKinnon, J. G., M. Ø. Nielsen, and M. D. Webb (2022c).  Leverage, influence, and the jackknife in clustered regression models: Reliable inference using summclust.  Working Paper.
  * Menzel (2021) Menzel, K. (2021).  Bootstrap with cluster-dependence in two or more dimensions.  Econometrica 89(5), 2143–2188.
  * Penrose and Yukich (2003) Penrose, M. D. and J. E. Yukich (2003).  Weak laws of large numbers in geometric probability.  Annals of Applied Probability 13(1), 277–303.
  * Tabord-Meehan (2019) Tabord-Meehan, M. (2019).  Inference with dyadic data: Asymptotic behavior of the dyadic-robust t-statistic.  Journal of Business & Economic Statistics 37(4), 671–680.
  * Vainora (2020) Vainora, J. (2020).  Network dependence.  Working Paper.

##  Appendix A Mathematical Set-up

This section lays out the mathematical set-up of our model in more detail,
heavily drawing from Kojevnikov et al. (2021).

We first define a collection of pairs of sets of dyads. For any positive
integers $a$, $b$ and $s$, define

$\displaystyle\mathcal{P}_{M}(a,b;s)\coloneqq\\{(A,B):A,B\subset\mathcal{M}_{N},|A|=a,|B|=b,\rho_{M}(A,B)\geq
s\\},$

where

$\displaystyle\rho_{M}(A,B)\coloneqq\min_{m\in A}\min_{m^{\prime}\in
B}\rho_{M}(m,m^{\prime}),$ (16)

with $\rho_{M}(m,m^{\prime})$ denoting the geodesic distance between dyads $m$
and $m^{\prime}$, i.e., the smallest number of adjacent dyads between dyads
$m$ and $m^{\prime}$. In words, the set $\mathcal{P}_{M}(a,b;s)$ collects all
two distinct sets of active dyads whose sizes are $a$ and $b$ and that have no
dyads in common.

Next we consider a collection of bounded Lipschitz functions. Define

$\displaystyle\mathcal{L}_{K}\coloneqq\\{\mathcal{L}_{K,c}:c\in\mathbbm{N}\\},$

where

$\displaystyle\mathcal{L}_{K,c}\coloneqq\\{f:\mathbbm{R}^{K\times
c}\rightarrow\mathbbm{R}:\|f\|_{\infty}<\infty,\text{Lip}(f)<\infty\\},$

with $\|\cdot\|_{\infty}$ representing the supremum norm and $\text{Lip}(f)$
being the Lipschitz constant.313131It is immediate to see that $\mathbbm{R}$
is a normed space with respect to the Euclidean norm, while the
$\mathbbm{R}^{K\times c}$ can be equipped with the norm
$\rho_{c}(x,y)\coloneqq\sum_{\ell=1}^{c}\|x_{\ell}-y_{\ell}\|$ where
$x,y\in\mathbbm{R}^{K\times c}$ and
$\|z\|\coloneqq(z^{\prime}z)^{\frac{1}{2}}$, thereby the Lipschitz constant is
defined as $\text{Lip}(f)\coloneqq\min\\{w\in\mathbbm{R}:|f(x)-f(y)|\leq
w\rho_{c}(x,y)\ \forall x,y\in\mathbbm{R}^{K\times c}\\}$. In words, the set
$\mathcal{L}_{K,c}$ collects all the bounded Lipschitz functions on
$\mathbbm{R}^{K\times c}$ and the set $\mathcal{L}_{K}$ moreover gathers such
sets with respect to $c\in\mathbbm{N}$.

Lastly, we write

$\displaystyle Y_{M,A}\coloneqq(Y_{M,m})_{m\in A},$

and $Y_{M,B}$ is analogously defined. Let $\\{\mathcal{C}_{M}\\}_{M\geq 1}$
denote a sequence of $\sigma$-algebras and be suppressed as
$\\{\mathcal{C}_{M}\\}$.

The network dependent random variables are characterized by the upper bound of
their covariances, first defined in Definition 2.2 of Kojevnikov et al.
(2021).

######  Definition A.1 (Conditional $\psi$-Dependence given
$\\{\mathcal{C}_{M}\\}$).

A triangular array $\\{Y_{M,m}\in\mathbbm{R}^{K}:M\geq
1,m\in\\{1,\ldots,M\\}\\}$ is called conditionally $\psi$-dependent given
$\\{\mathcal{C}_{M}\\}$, if for each $M\in\mathbbm{N}$, there exist a
$\mathcal{C}_{M}$-measurable sequence
$\theta_{M}\coloneqq\\{\theta_{M,s}\\}_{s\geq 0}$ with $\theta_{M,0}=1$, and a
collection of nonrandom function $(\psi_{a,b})_{a,b\in\mathbbm{N}}$ where
$\psi_{a,b}:\mathcal{L}_{K,a}\times\mathcal{L}_{K,b}\rightarrow[0.\infty)$,
such that for all $(A,B)\in\mathcal{P}_{M}(a,b;s)$ with $s>0$ and all
$f\in\mathcal{L}_{K,a}$ and $g\in\mathcal{L}_{K,b}$,

$\displaystyle\big{|}Cov\big{(}f(Y_{M,A}),g(Y_{M,B})\mid\mathcal{C}_{M}\big{)}\big{|}\leq\psi_{a,b}(f,g)\theta_{M,s}\
\ \ \ \ a.s.$

Intuitively, this definition states that the upper bound must be decomposed
into two components. The first part $\psi_{a,b}(f,g)$ is deterministic and
depends on nonlinear Lipschitz functions $f$ and $g$. The other component
$\theta_{M,s}$ is stochastic and depends only on the distance of the random
variables on the underlying network. The former, nonrandom component reflects
the scaling of the random variables as well as that of the Lipschitz
transformations, while the latter random part stands for the covariability of
the two random variables. We call $\theta_{M,s}$ the dependence coefficient.
We follow Kojevnikov et al. (2021) in assuming boundedness for these two
components.

######  Assumption A.1 (Kojevnikov et al. (2021), Assumption 2.1).

The triangular array $\\{Y_{M,m}\in\mathbbm{R}^{K}:M\geq
1,m\in\\{1,\ldots,M\\}\\}$ is conditionally $\psi$-dependent given
$\\{\mathcal{C}\\}$ with the dependence coefficients $\\{\theta_{M,s}\\}$
satisfying the following conditions: (a) there exists a constant $C>0$ such
that $\psi_{a,b}(f,g)\leq C\times
ab\big{(}\|f\|_{\infty}+\text{Lip}(f)\big{)}\big{(}\|g\|_{\infty}+\text{Lip}(g)\big{)}$;
(b) $\sup_{M\geq 1}\max_{s\geq 1}\theta_{M,s}<\infty$ a.s.

Assumption A.1 is maintained throughout the paper and employed to show
asymptotic properties of our estimators such as the consistency and asymptotic
normality, and the consistency of the network-robust variance estimator for
dyadic data.

To prove Theorem 3.2, we make use of the following assumptions which
strengthen their counterparts in the proof of consistency.

######  Assumption A.2 (Conditional Finite Moment of $\varepsilon_{m}$).

There exists $p>4$ such that
$\sup_{N\geq
1}\max_{m\in\mathcal{M}_{N}}E\big{[}|\varepsilon_{m}|^{p}\mid\mathcal{C}_{M}\big{]}<\infty$
a.s.

Assumption A.2 allows for the same interpretation as Assumption 3.1, i.e., the
random error term $\varepsilon_{m}$ cannot be too large, conditional on a
common component. This assumption, however, is more stringent than the
previous one because it now requires the finiteness of a higher moment of
$\varepsilon_{m}$. When coupled with the full rank condition (Assumption 2.2),
Assumption A.2 implies Assumption 3.3 of Kojevnikov et al. (2021) for each
$Y_{M,m}^{u}$ with $u\in\\{1,\ldots,K\\}$.

Our paper considers two types of variances: namely, unconditional
($\tau^{2}_{M}$, introduced in the main text) and the conditional variance of
$S_{M}^{u}$, $u\in\\{1,\ldots,K\\}$, given $\mathcal{C}_{M}$, denoted by
$\sigma_{M}^{2}\coloneqq Var(S_{M}^{u}\mid\mathcal{C}_{M})$, which is based on
the conditional assumptions described in the text. The following assumption
links the growth rates of the two of them.

######  Assumption A.3 (Growth Rates of Variances).

$\frac{\sigma_{M}^{2}}{\tau_{M}^{2}}\stackrel{{\scriptstyle
a.s.}}{{\rightarrow}}1$ as $N\rightarrow\infty$.

Moreover, we maintain Assumption 3.4 of Kojevnikov et al. (2021), which
regulates the conditional variance and the dependence coefficients, as
appropriate.

######  Assumption A.4 (Kojevnikov et al. (2021), Assumption 3.4).

There exists a positive sequence $r_{M}\rightarrow\infty$ such that for
$k=1,2$,

$\displaystyle\frac{M^{2}\theta_{M,r_{M}}^{1-1/p}}{\sigma_{M}}\stackrel{{\scriptstyle
a.s.}}{{\rightarrow}}0,$ $\displaystyle\frac{M}{\sigma_{M}^{2+k}}\sum_{s\geq
0}c_{n}(s,r_{M};k)\theta_{M,s}^{1-\frac{2+k}{p}}\stackrel{{\scriptstyle
a.s.}}{{\rightarrow}}0,$

as $M\rightarrow\infty$, where $p>4$ is the same as the one that appears in
Assumption A.2.

Similar to Assumption 3.2, this assumption binds the covariance of the random
variables, the dependence reflected in the dependence coefficients, and the
underlying network. The first part of this assumption requires that
$\sigma_{M}$ grows at least as fast as $M$ does. Kojevnikov et al. (2021)
claim that it is satisfied if the long-run variance $\frac{\sigma_{M}^{2}}{M}$
is bounded away from $c^{2}$ for all $M\geq 1$ for some constant $c$
independent of $M$, i.e., $\sigma_{M}\geq c\sqrt{M}$. The second part can be
interpreted in an analogous manner to Assumption 3.2. That is, given
$\sigma_{M}$ growing at least at the same rate of $M$, the composite of the
density of the network and the magnitude of the correlations of the random
variables must decay fast enough.

##  Appendix B Proofs of Auxiliary Results

###  B.1 Identification of $\beta$

######  Proposition B.1 (Identification).

Under Assumption 2.2, the regression parameter $\beta$ in (1) is identified.

###### Proof.

For each $m\in\mathcal{M}_{N}$, premultiply the model (1) by $x_{M,m}$ to
obtain

$\displaystyle
x_{M,m}y_{M,m}=x_{M,m}x_{M,m}^{\prime}\beta+x_{M,m}\varepsilon_{M,m}\hskip
28.45274pt\forall m\in\mathcal{M}_{N}.$

Taking the expectation with respect to
$\\{(x_{M,m},y_{M,m},\varepsilon_{M,m})\\}_{m\in\mathcal{M}_{N}}$ implies:

$\displaystyle E\Big{[}x_{M,m}y_{M,m}\Big{]}$
$\displaystyle=E\Big{[}x_{M,m}x_{M,m}^{\prime}\Big{]}\beta+E\Big{[}x_{M,m}\varepsilon_{M,m}\Big{]}.$

The second term on the right hand side is equal to 0, due to Assumption 2.2
(d). Next, Assumption 2.2 (c) ensures existence of the inverse of the
expectation term in the first term of the right hand side, ensuring
identification. ∎

###  B.2 Consistency of $\hat{\beta}$

Proof of Theorem 3.1: From (5), (9) and (10), we can write

$\displaystyle\hat{\beta}-\beta$
$\displaystyle=\Big{(}\sum_{j\in\mathcal{M}_{N}}x_{M,j}x_{M,j}^{\prime}\Big{)}^{-1}\sum_{m\in\mathcal{M}_{N}}x_{M,m}\varepsilon_{M,m}$
$\displaystyle=\Big{(}\frac{1}{M}\sum_{j\in\mathcal{M}_{N}}x_{M,j}x_{M,j}^{\prime}\Big{)}^{-1}\frac{1}{M}\sum_{m\in\mathcal{M}_{N}}Y_{M,m}$
$\displaystyle=\frac{1}{M}\sum_{m\in\mathcal{M}_{N}}\Big{(}\frac{1}{M}\sum_{j\in\mathcal{M}_{N}}x_{M,j}x_{M,j}^{\prime}\Big{)}^{-1}Y_{M,m}.$

Define
$\tilde{Y}_{M,m}\coloneqq\Big{(}\frac{1}{M}\sum_{j\in\mathcal{M}_{N}}x_{M,j}x_{M,j}^{\prime}\Big{)}^{-1}Y_{M,m}$
and let $\tilde{Y}_{M,m}^{u}$ be the $u$-th entry of $\tilde{Y}_{M,m}$. That
is,

$\displaystyle\tilde{Y}_{M,m}^{u}$ $\displaystyle=D^{u}Y_{M,m}$
$\displaystyle=D^{u}x_{M,m}\varepsilon_{M,m},$

where $D^{u}$ stands for the $u$-th row of the matrix
$\big{(}\frac{1}{M}\sum_{j\in\mathcal{M}_{N}}x_{M,j}x_{M,j}^{\prime}\big{)}^{-1}$.
Moreover, let $\hat{\beta}^{u}$ and $\beta^{u}$, respectively, denote the
$u$-th entry of $\hat{\beta}$ and $\beta$, so that we can write

$\displaystyle\hat{\beta}^{u}-\beta^{u}$
$\displaystyle=\frac{1}{M}\sum_{m\in\mathcal{M}_{N}}\tilde{Y}_{M,m}^{u},$

In light of Assumption 2.2 (d), it holds that for any $N>0$ and for each
$m\in\mathcal{M}_{N}$

$\displaystyle E\big{[}\tilde{Y}_{M,m}^{u}\mid\mathcal{C}_{M}\big{]}$
$\displaystyle=D^{u}x_{M,m}\underbrace{E\big{[}\varepsilon_{M,m}\mid\mathcal{C}_{M}\big{]}}_{0}$
$\displaystyle=0.$

By Theorem 3.1 of Kojevnikov et al. (2021),
$\Big{\|}\frac{1}{M}\sum_{m\in\mathcal{M}_{N}}\Big{(}\tilde{Y}_{M,m}^{u}-\underbrace{E\big{[}\tilde{Y}_{M,m}^{u}\mid\mathcal{C}_{M}\big{]}}_{0}\Big{)}\Big{\|}_{\mathcal{C}_{M},1}\stackrel{{\scriptstyle
a.s.}}{{\rightarrow}}0$. Hence,

$\displaystyle\Big{\|}\frac{1}{M}\sum_{m\in\mathcal{M}_{N}}\tilde{Y}_{M,m}^{u}\Big{\|}_{\mathcal{C}_{M},1}\stackrel{{\scriptstyle
a.s.}}{{\rightarrow}}0\qquad M\rightarrow\infty,$

so that

$\displaystyle E\big{[}\big{|}\hat{\beta}^{u}-\beta^{u}\big{|}\big{]}$
$\displaystyle=E\Big{[}E\big{[}|\hat{\beta}^{u}-\beta^{u}|\mid\mathcal{C}_{M}\big{]}\Big{]}$
$\displaystyle=E\Big{[}\|\hat{\beta}^{u}-\beta^{u}\|_{\mathcal{C}_{M},1}\Big{]}$
$\displaystyle=E\Big{[}\Big{\|}\frac{1}{M}\sum_{m\in\mathcal{M}_{N}}\tilde{Y}_{M,m}^{u}\Big{\|}_{\mathcal{C}_{M},1}\Big{]}$
$\displaystyle\rightarrow{}0\qquad M\rightarrow\infty,$

where the last implication is a consequence of the Dominated Convergence
Theorem. In view of Assumption 2.1, this is true also with respect to $N$
going to infinity.

Since it holds by the Markov inequality that for any $c>0$

$\displaystyle\mathrm{Pr}\big{(}|\hat{\beta}^{u}-\beta^{u}|>c\big{)}\leq\frac{E\big{[}|\hat{\beta}^{u}-\beta^{u}|\big{]}}{c},$

it then follows that

$\displaystyle\mathrm{Pr}\big{(}|\hat{\beta}^{u}-\beta^{u}|>c\big{)}\rightarrow
0,$

as $N\rightarrow\infty$. Hence we have

$\displaystyle|\hat{\beta}^{u}-\beta^{u}|\stackrel{{\scriptstyle
p}}{{\rightarrow}}0\qquad as\quad N\rightarrow\infty.$

Finally, we can invoke the Cramér-Wold device to obtain

$\displaystyle\|\hat{\beta}-\beta\|_{2}\stackrel{{\scriptstyle
p}}{{\rightarrow}}0\qquad as\quad N\rightarrow\infty,$

as desired. $\square$

###  B.3 Lemma

Here we establish a lemma that is used repeatedly throughout the proofs in
this paper.

######  Lemma B.1.

Define
$A\coloneqq\lim_{N\rightarrow\infty}\frac{1}{M}\sum_{k\in\mathcal{M}}E\Big{[}x_{M,k}x_{M,k}^{\prime}\Big{]}$
and assume that Assumptions 2.1 and 3.3 hold.

  * (i)

$A^{-1}\coloneqq\lim_{N\rightarrow\infty}\Big{(}\frac{1}{M}\sum_{k\in\mathcal{M}}E\big{[}x_{M,k}x_{M,k}^{\prime}\big{]}\Big{)}^{-1}$
exists with finite elements and positive definite.

  * (ii)

Suppose, moreover, that Assumption 4.2 holds. Then,
$\Big{\|}\Big{(}\frac{1}{M}\sum_{k\in\mathcal{M}_{N}}x_{M,k}x_{M,k}^{\prime}\Big{)}^{-1}-\Big{(}\frac{1}{M}\sum_{k\in\mathcal{M}_{N}}E\big{[}x_{M,k}x_{M,k}^{\prime}\big{]}\Big{)}^{-1}\Big{\|}_{F}\stackrel{{\scriptstyle
p}}{{\rightarrow}}0$.

###### Proof.

$(i)$ The fact that it is positive definite follows from Assumption 3.3. The
fact that the elements are finite is proved by considering element-by-element
convergence. Let $x_{k,i}$ denote the $i$-th element of $x_{M,k}$. Then the
$(i,j)$ entry of
$\frac{1}{M}\sum_{k\in\mathcal{M}}E\big{[}x_{M,k}x_{M,k}^{\prime}\big{]}$ is
given by: $\frac{1}{M}\sum_{k\in\mathcal{M}}E\big{[}x_{k,i}x_{k,j}\big{]}$.

We write the $(i,j)$ entry of $A$ as $A_{i,j}$.

From Assumption 3.3, there exists a nonnegative finite constant $C_{0,1}$ such
that

$\displaystyle C_{0,1}=\sup_{N\geq
1}\max_{m\in\mathcal{M}_{N}}E\big{[}x_{m,i}x_{m,j}\big{]},$

so that

$\displaystyle A_{i,j}$
$\displaystyle=\lim_{N\rightarrow\infty}\frac{1}{M}\sum_{k\in\mathcal{M}_{N}}\underbrace{E\big{[}x_{k,i}x_{k,j}\big{]}}_{\leq
C_{0,1}}$
$\displaystyle\leq\lim_{N\rightarrow\infty}\frac{1}{M}\sum_{k\in\mathcal{M}_{N}}C_{0,1}$
$\displaystyle=C_{0,1}\lim_{N\rightarrow\infty}\frac{1}{M}\underbrace{\sum_{k\in\mathcal{M}_{N}}1}_{M}$
$\displaystyle=C_{0,1}.$

Hence $A_{i,j}$ exists with being finite. By repeating the same argument for
all $i,j=1,\ldots,K$, it holds that $A$ exists with finite elements.

$(ii)$ To begin with, observe that

$\displaystyle\Big{\|}\Big{(}\frac{1}{M}\sum_{k\in\mathcal{M}_{N}}x_{M,k}x_{M,k}^{\prime}\Big{)}^{-1}-\Big{(}\frac{1}{M}\sum_{k\in\mathcal{M}_{N}}E\big{[}x_{M,k}x_{M,k}^{\prime}\big{]}\Big{)}^{-1}\Big{\|}_{F}$
$\displaystyle=\Big{\|}\Big{(}\frac{1}{M}\sum_{k\in\mathcal{M}_{N}}x_{M,k}x_{M,k}^{\prime}\Big{)}^{-1}-A^{-1}+A^{-1}-\Big{(}\frac{1}{M}\sum_{k\in\mathcal{M}_{N}}E\big{[}x_{M,k}x_{M,k}^{\prime}\big{]}\Big{)}^{-1}\Big{\|}_{F}$
$\displaystyle\leq\Big{\|}\Big{(}\frac{1}{M}\sum_{k\in\mathcal{M}_{N}}x_{M,k}x_{M,k}^{\prime}\Big{)}^{-1}-A^{-1}\Big{\|}_{F}+\Big{\|}A^{-1}-\Big{(}\frac{1}{M}\sum_{k\in\mathcal{M}_{N}}E\big{[}x_{M,k}x_{M,k}^{\prime}\big{]}\Big{)}^{-1}\Big{\|}_{F}.$

Note that convergence of the second term follows from (i). Hence, we wish to
prove that:

$\displaystyle\Big{\|}\frac{1}{M}\sum_{k\in\mathcal{M}_{N}}x_{M,k}x_{M,k}^{\prime}-A\Big{\|}_{F}\stackrel{{\scriptstyle
p}}{{\rightarrow}}0$

To do so, we follow a strategy employed in Aronow et al. (2015) and Tabord-
Meehan (2019). In light of (i), it remains to show

$\displaystyle
Var\Big{(}\frac{1}{M}\sum_{k\in\mathcal{M}_{N}}x_{M,k}x_{M,k}^{\prime}\Big{)}\rightarrow
0.$

As in (i), we consider the element-by-element convergence, using the same
notation. The variance can be expressed as a sum of covariances:

$\displaystyle
Var\Big{(}\frac{1}{M}\sum_{k\in\mathcal{M}_{N}}x_{k,i}x_{k,j}\Big{)}$
$\displaystyle=\frac{1}{M^{2}}\sum_{m\in\mathcal{M}_{N}}\sum_{m^{\prime}\in\mathcal{M}_{N}}Cov\big{(}x_{m,i}x_{m,j},x_{m^{\prime},i}x_{m^{\prime},j}\big{)}$
$\displaystyle=\frac{1}{M^{2}}\sum_{s\geq
0}\sum_{m\in\mathcal{M}_{N}}\sum_{m^{\prime}\in\mathcal{M}_{N}^{\partial}(m;s)}Cov\big{(}x_{m,i}x_{m,j},x_{m^{\prime},i}x_{m^{\prime},j}\big{)}.$

Again from Assumption 3.3, there exists a nonnegative finite constant
$C_{0,2}$ such that

$\displaystyle C_{0,2}=\sup_{N\geq
1}\max_{m,m^{\prime}\in\mathcal{M}_{N}}Cov\left(x_{m,i}x_{m,j},x_{m^{\prime},i}x_{m^{\prime},j}\right).$

Hence,

$\displaystyle
Var\Big{(}\frac{1}{M}\sum_{k\in\mathcal{M}_{N}}x_{k,i}x_{k,j}\Big{)}$
$\displaystyle\leq\frac{1}{M^{2}}\sum_{s\geq
0}\sum_{m\in\mathcal{M}_{N}}\sum_{m^{\prime}\in\mathcal{M}_{N}^{\partial}(m;s)}C_{0,2}$
$\displaystyle=\frac{C_{0,2}}{M^{2}}\sum_{s\geq
0}\sum_{m\in\mathcal{M}_{N}}\sum_{m^{\prime}\in\mathcal{M}_{N}^{\partial}(m;s)}1$
$\displaystyle=\frac{C_{0,2}}{M^{2}}\sum_{s\geq 0}M\delta_{M}^{\partial}(s;1)$
$\displaystyle=C_{0,2}\underbrace{\frac{1}{M}\sum_{s\geq
0}\delta_{M}^{\partial}(s;1)}_{\rightarrow 0}$ $\displaystyle\rightarrow 0,$

where the last implication is due to Assumption 4.2. By repeating the same
argument for all $i,j=1,\ldots,K$, we obtain

$\displaystyle
Var\Big{(}\frac{1}{M}\sum_{k\in\mathcal{M}_{N}}x_{M,k}x_{M,k}^{\prime}\Big{)}\rightarrow
0.$

Now, by the Chebyshev’s inequality, we arrive at

$\displaystyle\Big{\|}\frac{1}{M}\sum_{k\in\mathcal{M}_{N}}x_{M,k}x_{M,k}^{\prime}-A\Big{\|}_{F}\stackrel{{\scriptstyle
p}}{{\rightarrow}}0.$

Furthermore, applying the Continuous Mapping Theorem yields

$\displaystyle\Big{\|}\Big{(}\frac{1}{M}\sum_{k\in\mathcal{M}_{N}}x_{M,k}x_{M,k}^{\prime}\Big{)}^{-1}-A^{-1}\Big{\|}_{F}\stackrel{{\scriptstyle
p}}{{\rightarrow}}0,$

obtaining the result. Therefore,

$\displaystyle\Big{\|}\Big{(}\frac{1}{M}\sum_{k\in\mathcal{M}_{N}}x_{M,k}x_{M,k}^{\prime}\Big{)}^{-1}-\Big{(}\frac{1}{M}\sum_{k\in\mathcal{M}_{N}}E\big{[}x_{M,k}x_{M,k}^{\prime}\big{]}\Big{)}^{-1}\Big{\|}_{F}\stackrel{{\scriptstyle
p}}{{\rightarrow}}0,$

as desired. ∎

###  B.4 Asymptotic Normality of $\hat{\beta}$

In this subsection, we prove Theorem 3.2 under a slightly milder condition
than Assumption A.3.

######  Assumption B.1 (Growth Rates of Variances).

There exists a sequence of (possibly random) positive numbers,
$\\{\pi_{N,M}\\}_{N>0}$, such that

$\displaystyle\frac{\sigma_{M}^{2}}{\pi_{N,M}\tau_{M}^{2}}\stackrel{{\scriptstyle
a.s.}}{{\rightarrow}}1\quad\quad\text{as}\quad N\rightarrow\infty.$

When $\pi_{N,M}=1$, this assumption simplifies to Assumption A.3, which is
used for the results in the main text.

When Assumption A.3 is replaced by Assumption B.1, Assumption 3.4 must also be
modified accordingly.

######  Assumption B.2.

$\lim_{N\rightarrow\infty}\frac{N\pi_{N,M}}{M^{2}}\sum_{m\in\mathcal{M}_{N}}\sum_{m^{\prime}\in\mathcal{M}_{N}}E\big{[}\varepsilon_{M,m}\varepsilon_{M,m^{\prime}}x_{M,m}x_{M,m^{\prime}}^{\prime}\big{]}$
exists with finite elements.

An important comparison of Assumption B.2 can be made to the variety of
assumptions used in the literature.

######  Remark B.1.

The requirement on the behavior of $AVar(\hat{\beta})$ mirrors Assumptions
2.4, 2.5 and 2.6 of Tabord-Meehan (2019): Assumption B.2 boils down to his
Assumption 2.4, if it is well-defined with $\pi_{N,M}=\frac{M}{N}$; it reduces
to his Assumption 2.5, if it is compatible with $\pi_{N,M}=\frac{M}{N^{2}}$;
and it coincides with Assumption 2.6, if it is maintained with
$\pi_{N,M}=\frac{M}{N^{r+1}}$ for $r\in[0,1]$. Moreover, if
$AVar(\hat{\beta})$ is well-defined for $\pi_{N,M}=1$, the expression (11)
simplifies to the assumption that appears in Lemma 1 of Aronow et al. (2015).

Proof of Theorem 3.2: From (10),

$\displaystyle\sqrt{N}(\hat{\beta}-\beta)$
$\displaystyle=\Big{(}\frac{1}{M}\sum_{j\in\mathcal{M}_{N}}x_{M,j}x_{M,j}^{\prime}\Big{)}^{-1}\frac{\sqrt{N}}{M}\sum_{m\in\mathcal{M}_{N}}Y_{M,m}$
$\displaystyle=\Big{(}\frac{1}{M}\sum_{j\in\mathcal{M}_{N}}x_{M,j}x_{M,j}^{\prime}\Big{)}^{-1}\frac{\sqrt{N}}{M}S_{M}.$

Since
$\big{(}\frac{1}{M}\sum_{j\in\mathcal{M}_{N}}x_{M,j}x_{M,j}^{\prime}\big{)}^{-1}$
converges to a well-defined limit (Lemma B.1), the asymptotic distribution of
$\sqrt{N}(\hat{\beta}-\beta)$ is dictated by that of
$\frac{\sqrt{N}}{M}S_{M}$.

First of all, we prove

$\displaystyle\frac{S_{M}^{u}}{\sigma_{M}}\stackrel{{\scriptstyle
d}}{{\rightarrow}}\mathcal{N}(0,1),$

as $N\rightarrow\infty$. Consider the scenario that $N\rightarrow\infty$, in
which Assumption 2.1 implies $M\rightarrow\infty$. Denote
$\tilde{S}_{M}^{u}\coloneqq\frac{S_{M}^{u}}{\sigma_{M}}$. Under Assumptions
A.2 and A.4, it holds by Theorem 3.2 of Kojevnikov et al. (2021) that for any
$\epsilon>0$, there exists $M_{0}>0$ such that for each $M>M_{0}$ and for each
$x\in\mathbbm{R}$,

$\displaystyle\sup_{t\in\mathbbm{R}}\big{|}\mathrm{Pr}(\tilde{S}_{M}^{u}\leq
t\mid X=x)-\Phi(t)\big{|}<\epsilon,$ (17)

where $\Phi(\cdot)$ is the CDF of a standard Normal distribution. Then, by the
law of total probability, we have

$\displaystyle\big{|}\mathrm{Pr}(\tilde{S}_{M}^{u}\leq t)-\Phi(t)\big{|}$
$\displaystyle=\Big{|}\int\mathrm{Pr}(\tilde{S}_{M}^{u}\leq t\mid
X=x)dF_{X}(x)-\Phi(t)\Big{|}$
$\displaystyle=\Big{|}\int\mathrm{Pr}(\tilde{S}_{M}^{u}\leq t\mid
X=x)-\Phi(t)dF_{X}(x)\Big{|}$
$\displaystyle\leq\int\big{|}\mathrm{Pr}(\tilde{S}_{M}^{u}\leq t\mid
X=x)-\Phi(t)\big{|}dF_{X}(x)$
$\displaystyle\leq\int\sup_{t\in\mathbbm{R}}\big{|}\mathrm{Pr}(\tilde{S}_{M}^{u}\leq
t\mid X=x)-\Phi(t)\big{|}dF_{X}(x),$ (18)

where $F_{X}(\cdot)$ denotes the probability distribution function of $X$. Now
pick arbitrarily $\epsilon>0$. Then there exists $M_{0}>0$ such that for each
$M>M_{0}$

$\displaystyle\big{|}\mathrm{Pr}(\tilde{S}_{M}^{u}\leq t)-\Phi(t)\big{|}$
$\displaystyle\leq\int\underbrace{\sup_{t\in\mathbbm{R}}\big{|}\mathrm{Pr}(\tilde{S}_{M}^{u}\leq
t\mid X=x)-\Phi(t)\big{|}}_{<\epsilon}dF_{X}(x)$
$\displaystyle\leq\int\epsilon dF_{X}(x)$ $\displaystyle\leq\epsilon,$ (19)

where the first and second inequalities come from (B.4) and (17),
respectively. Since the right hand side of (B.4) does not depend on $t$, we
then have that for each $M>M_{0}$,

$\displaystyle\sup_{t\in\mathbbm{R}}\big{|}\mathrm{Pr}(\tilde{S}_{M}^{u}\leq
t)-\Phi(t)\big{|}$ $\displaystyle\leq\epsilon,$

which implies

$\displaystyle\sup_{t\in\mathbbm{R}}\big{|}\mathrm{Pr}(\tilde{S}_{M}^{u}\leq
t)-\Phi(t)\big{|}$ $\displaystyle\rightarrow 0\qquad as\quad
M\rightarrow\infty,$

We have then shown that

$\displaystyle\sup_{t\in\mathbbm{R}}\big{|}\mathrm{Pr}(\tilde{S}_{M}^{u}\leq
t)-\Phi(t)\big{|}$ $\displaystyle\rightarrow 0\qquad as\quad
N\rightarrow\infty,$

from which we obtain

$\displaystyle\frac{S_{M}^{u}}{\sigma_{M}}\stackrel{{\scriptstyle
d}}{{\rightarrow}}\mathcal{N}(0,1)\qquad as\quad N\rightarrow\infty.$

Next this can be combined with Assumption B.1 by using the Slutsky’s Theorem,
yielding that

$\displaystyle\frac{S_{M}^{u}}{\tau_{M}\sqrt{\pi_{N,M}}}\stackrel{{\scriptstyle
d}}{{\rightarrow}}\mathcal{N}(0,1)\qquad as\quad N\rightarrow\infty.$

Moreover, applying the Cramér-Wold device gives

$\displaystyle\frac{\tau_{M}^{-1}}{\sqrt{\pi_{N,M}}}S_{M}\stackrel{{\scriptstyle
d}}{{\rightarrow}}\mathcal{N}(0,I_{K})\qquad as\quad N\rightarrow\infty,$

where $I_{K}$ is the $K\times K$ identity matrix and $\tau_{M}$ is understood
as the variance-covariance matrix.323232To save notation, we use the same
$\tau_{M}$ to denote the case of one-dimensional parameter and the case of
multiple-dimensional parameters.

Now notice that we have

$\displaystyle\sqrt{N}(\hat{\beta}-\beta)$
$\displaystyle=\Big{(}\frac{1}{M}\sum_{j\in\mathcal{M}_{N}}x_{M,j}x_{M,j}^{\prime}\Big{)}^{-1}\frac{\sqrt{N}}{M}\tau_{M}\sqrt{\pi_{N,M}}\underbrace{\frac{\tau_{M}^{-1}}{\sqrt{\pi_{N,M}}}S_{M,m}}_{\stackrel{{\scriptstyle
d}}{{\rightarrow}}\mathcal{N}(0,I_{K})}.$

Hence we obtain

$\displaystyle\sqrt{N}\big{(}\hat{\beta}-\beta\big{)}\stackrel{{\scriptstyle
d}}{{\rightarrow}}\mathcal{N}(0,AVar(\hat{\beta}))\qquad as\quad
N\rightarrow\infty,$

where

$\displaystyle
AVar(\hat{\beta})\coloneqq\lim_{N\rightarrow\infty}N\pi_{N,M}\Big{(}\sum_{k\in\mathcal{M}_{N}}E\big{[}x_{M,k}x_{M,k}^{\prime}\big{]}\Big{)}^{-1}\Big{(}\sum_{m\in\mathcal{M}_{N}}\sum_{m^{\prime}\in\mathcal{M}_{N}}E\big{[}\varepsilon_{M,m}\varepsilon_{M,m^{\prime}}x_{M,m}x_{M,m^{\prime}}^{\prime}\big{]}\Big{)}\Big{(}\sum_{k\in\mathcal{M}_{N}}E\big{[}x_{M,k}x_{M,k}^{\prime}\big{]}\Big{)}^{-1},$

which is well-defined due to Lemma B.1 (i) along with Assumption B.2. When
$\pi_{N,M}=1$, this is the result in the main text under Assumptions 3.3 and
3.4. $\square$

##  Appendix C Proofs of Main Theorems and Results

###  C.1 Lemma

In the proof of Theorem 4.1, we make use of the following lemma from
Kojevnikov et al. (2021), p.903:

######  Lemma C.1.

Define

$\displaystyle
H_{M}(s,r)\coloneqq\\{(m,j,k,l)\in\mathcal{M}_{N}^{4}:j\in\mathcal{M}_{N}(m;r),\
l\in\mathcal{M}_{N}(k;r),\ \rho_{M}(\\{m,j\\},\\{k,l\\})=s\\}.$

Then

$\displaystyle|H_{M}(s,r)|\leq 4Mc_{M}(s,r;2).$

###  C.2 Consistency of $\widehat{Var}(\hat{\beta})$

Proof of Theorem 4.1: Denote the variance of $\frac{S_{M}}{\sqrt{M}}$ as
$V_{N,M}\coloneqq Var\big{(}\frac{S_{M}}{\sqrt{M}}\big{)}$. It can readily be
shown that $V_{N,M}$ takes the form of $V_{N,M}=\sum_{s\geq
0}\Omega_{N,M}(s)$, where

$\displaystyle\Omega_{N,M}(s)\coloneqq\frac{1}{M}\sum_{m\in\mathcal{M}_{N}}\sum_{j\in\mathcal{M}_{N}^{\partial}(m;s)}E\big{[}Y_{M,m}Y_{M,j}^{\prime}\big{]}.$

Following Kojevnikov et al. (2021), we define the kernel heteroskedasticity
and autocorrelation consistent (HAC) estimator of $V_{N,M}$ as
$\hat{V}_{N,M}\coloneqq\sum_{s\geq 0}\omega_{M}(s)\hat{\Omega}_{N,M}(s)$,
where $\omega_{M}(s)\coloneqq\omega\big{(}\frac{s}{b_{M}}\big{)}$ and

$\displaystyle\hat{\Omega}_{N,M}(s)\coloneqq\frac{1}{M}\sum_{m\in\mathcal{M}_{N}}\sum_{j\in\mathcal{M}_{N}^{\partial}(m;s)}\hat{Y}_{M,m}\hat{Y}_{M,j}^{\prime}.$

Moreover, we define an empirical analogue of $V_{N,M}$, though infeasible, by
$\tilde{V}_{N,M}\coloneqq\sum_{s\geq 0}\omega_{M}(s)\tilde{\Omega}_{N,M}(s)$,
where

$\displaystyle\tilde{\Omega}_{N,M}(s)\coloneqq\frac{1}{M}\sum_{m\in\mathcal{M}_{N}}\sum_{j\in\mathcal{M}_{N}^{\partial}(m;s)}Y_{M,m}Y_{M,j}^{{}^{\prime}}.$

Additionally, we denote a conditional version of $V_{N,M}$ by
$V_{N,M}^{c}\coloneqq
Var\big{(}\frac{S_{M}}{\sqrt{M}}\mid\mathcal{C}_{M}\big{)}$, i.e.,
$V_{N,M}^{c}=\sum_{s\geq 0}\Omega_{N,M}^{c}(s)$, where

$\displaystyle\Omega_{N,M}^{c}(s)\coloneqq\frac{1}{M}\sum_{m\in\mathcal{M}_{N}}\sum_{j\in\mathcal{M}_{N}^{\partial}(m;s)}E\big{[}Y_{M,m}Y_{M,j}^{{}^{\prime}}\mid\mathcal{C}_{M}\big{]}.$

Notice that since $E\big{[}Y_{M,m}\mid\mathcal{C}_{M}\big{]}=0$ $a.s.$, it
follows from the law of total variance that
$V_{N,M}=E\big{[}V_{N,M}^{c}\big{]}$. Notice furthermore that it holds that

$\displaystyle
Var(\hat{\beta})=\frac{N}{M}\Big{(}\frac{1}{M}\sum_{k\in\mathcal{M}_{N}}E\big{[}x_{M,k}x_{M,k}^{\prime}\big{]}\Big{)}^{-1}V_{N,M}\Big{(}\frac{1}{M}\sum_{k\in\mathcal{M}_{N}}E\big{[}x_{M,k}x_{M,k}^{\prime}\big{]}\Big{)}^{-1},$

and

$\displaystyle\widehat{Var}(\hat{\beta})=\frac{1}{M}\Big{(}\frac{1}{M}\sum_{k\in\mathcal{M}_{N}}x_{M,k}x_{M,k}^{\prime}\Big{)}^{-1}\hat{V}_{N,M}\Big{(}\frac{1}{M}\sum_{k\in\mathcal{M}_{N}}x_{M,k}x_{M,k}^{\prime}\Big{)}^{-1}.$

Since

$\displaystyle\|N\widehat{Var}(\hat{\beta})-Var(\hat{\beta})\|_{F}$
$\displaystyle=\frac{N}{M}\Bigg{\|}\Big{(}\frac{1}{M}\sum_{k\in\mathcal{M}_{N}}x_{M,k}x_{M,k}^{\prime}\Big{)}^{-1}\hat{V}_{N,M}\Big{(}\frac{1}{M}\sum_{k\in\mathcal{M}_{N}}x_{M,k}x_{M,k}^{\prime}\Big{)}^{-1}$
$\displaystyle\phantom{=\frac{N}{M}\Bigg{\|}}-\Big{(}\frac{1}{M}\sum_{k\in\mathcal{M}_{N}}E\big{[}x_{M,k}x_{M,k}^{\prime}\big{]}\Big{)}^{-1}V_{N,M}\Big{(}\frac{1}{M}\sum_{k\in\mathcal{M}_{N}}E\big{[}x_{M,k}x_{M,k}^{\prime}\big{]}\Big{)}^{-1}\Bigg{\|}_{F},$

and $\frac{N}{M}$ is bounded due to Assumption 2.1, it thus suffices to show
that

  * (i)

$\Big{\|}\big{(}\frac{1}{M}\sum_{k\in\mathcal{M}_{N}}x_{M,k}x_{M,k}^{\prime}\big{)}^{-1}-\big{(}\frac{1}{M}\sum_{k\in\mathcal{M}_{N}}E\big{[}x_{M,k}x_{M,k}^{\prime}\big{]}\big{)}^{-1}\Big{\|}_{F}\stackrel{{\scriptstyle
p}}{{\rightarrow}}0$;

  * (ii)

$\|\hat{V}_{N,M}-V_{N,M}\|_{F}\stackrel{{\scriptstyle p}}{{\rightarrow}}0$.

Part (i) is already shown in Lemma B.1 (ii). Hence, it remains to prove Part
(ii).

To begin with, observe that by the technique of add and subtract as well as
the triangular inequality,

$\displaystyle\|\hat{V}_{N,M}-V_{N,M}\|_{F}$
$\displaystyle=\|\hat{V}_{N,M}-\tilde{V}_{N,M}+\tilde{V}_{N,M}-V^{c}_{N,M}+V^{c}_{N,M}-V_{N,M}\|_{F}$
$\displaystyle\leq\|\hat{V}_{N,M}-\tilde{V}_{N,M}\|_{F}+\|\tilde{V}_{N,M}-V^{c}_{N,M}\|_{F}+\|V^{c}_{N,M}-V_{N,M}\|_{F}.$

We thus aim to prove

  * (1)

$\|V^{c}_{N,M}-V_{N,M}\|_{F}\stackrel{{\scriptstyle p}}{{\rightarrow}}0$;

  * (2)

$\|\tilde{V}_{N,M}-V^{c}_{N,M}\|_{F}\stackrel{{\scriptstyle
p}}{{\rightarrow}}0$;

  * (3)

$\|\hat{V}_{N,M}-\tilde{V}_{N,M}\|_{F}\stackrel{{\scriptstyle
p}}{{\rightarrow}}0$.

  * (1)

$\|V^{c}_{N,M}-V_{N,M}\|_{F}\stackrel{{\scriptstyle p}}{{\rightarrow}}0$:
The proof proceeds multiple steps:

    * (a)

$E\big{[}\big{\|}V_{N,M}^{c}-V_{N,M}\big{\|}_{F}^{2}\big{]}\rightarrow 0$;

    * (b)

$\big{\|}V_{N,M}^{c}-V_{N,M}\big{\|}_{F}\stackrel{{\scriptstyle
p}}{{\rightarrow}}0$,

in the sequel.

    * (a)

$E\big{[}\big{\|}V_{N,M}^{c}-V_{N,M}\big{\|}_{F}^{2}\big{]}\rightarrow 0$:
We prove this by showing the element-wise convergence. With a slight abuse of
notation, we denote the $(a,b)$ entry of $V_{N,M}^{c}$ and $V_{N,M}$ as
$V_{a,b}^{c}$ and $V_{a,b}$, respectively. Then it is enough to verify that

$\displaystyle E\big{[}(V_{a,b}^{c}-V_{a,b})^{2}\big{]}\rightarrow 0.$

Notice that $V_{a,b}^{c}$ and $V_{a,b}$ are given by

$\displaystyle V_{a,b}^{c}=\sum_{s\geq
0}\frac{1}{M}\sum_{m\in\mathcal{M}_{N}}\sum_{j\in\mathcal{M}_{N}^{\partial}(m;s)}E\big{[}Y_{m,a}Y_{j,b}\mid\mathcal{C}_{M}\big{]}$

and

$\displaystyle V_{a,b}=\sum_{s\geq
0}\frac{1}{M}\sum_{m\in\mathcal{M}_{N}}\sum_{j\in\mathcal{M}_{N}^{\partial}(m;s)}E\big{[}Y_{m,a}Y_{j,b}\big{]},$

where $Y_{m,a}$ and $Y_{m,b}$ stand for the $a$-th and $b$-th element of
$Y_{M,m}$, respectively. Note moreover that $E[V_{a,b}^{c}]=V_{a,b}$. Hence,
we can write

$\displaystyle E\big{[}(V_{a,b}^{c}-V_{a,b})^{2}\big{]}$
$\displaystyle=Var(V_{a,b}^{c})$
$\displaystyle=E\big{[}(V_{a,b}^{c})^{2}\big{]}-\big{(}V_{a,b}\big{)}^{2}$
$\displaystyle\leq E\big{[}(V_{a,b}^{c})^{2}\big{]}.$

Observe that

$\displaystyle E\big{[}(V_{a,b}^{c})^{2}\big{]}$
$\displaystyle=E\Big{[}\Big{(}\sum_{s\geq
0}\frac{1}{M}\sum_{m\in\mathcal{M}_{N}}\sum_{j\in\mathcal{M}_{N}^{\partial}(m;s)}E\big{[}Y_{m,a}Y_{j,b}\mid\mathcal{C}_{M}\big{]}\Big{)}^{2}\Big{]}$
$\displaystyle=E\Big{[}\frac{1}{M^{2}}\sum_{s\geq
0}\sum_{m\in\mathcal{M}_{N}}\sum_{j\in\mathcal{M}_{N}^{\partial}(m;s)}\sum_{t\geq
0}\sum_{k\in\mathcal{M}_{N}}\sum_{l\in\mathcal{M}_{N}^{\partial}(k;t)}E\big{[}Y_{m,a}Y_{j,b}\mid\mathcal{C}_{M}\big{]}E\big{[}Y_{k,a}Y_{l,b}\mid\mathcal{C}_{M}\big{]}\Big{]}$
$\displaystyle=\frac{1}{M^{2}}\sum_{s\geq
0}\sum_{m\in\mathcal{M}_{N}}\sum_{j\in\mathcal{M}_{N}^{\partial}(m;s)}\sum_{t\geq
0}\sum_{k\in\mathcal{M}_{N}}\sum_{l\in\mathcal{M}_{N}^{\partial}(k;t)}E\Big{[}E\big{[}Y_{m,a}Y_{j,b}\mid\mathcal{C}_{M}\big{]}E\big{[}Y_{k,a}Y_{l,b}\mid\mathcal{C}_{M}\big{]}\Big{]}.$

By the Cauchy-Schwartz inequality,

$\displaystyle
E\big{[}\varepsilon_{m}\varepsilon_{j}\mid\mathcal{C}_{M}\big{]}\leq\big{(}E\big{[}\varepsilon_{m}^{2}\mid\mathcal{C}_{M}\big{]}\big{)}^{\frac{1}{2}}\big{(}E\big{[}\varepsilon_{m}^{2}\mid\mathcal{C}_{M}\big{]}\big{)}^{\frac{1}{2}},$

it then follows from from Assumption 4.1 (a) that there exists an a.s.-bounded
function $\bar{C}_{1}$ such that
$E\big{[}\varepsilon_{m}\varepsilon_{j}\mid\mathcal{C}_{M}\big{]}\leq\bar{C}_{1}\
a.s.$ Similarly, we have an a.s.-bounded function $\bar{C}_{2}$ such that
$E\big{[}\varepsilon_{k}\varepsilon_{l}\mid\mathcal{C}_{M}\big{]}\leq\bar{C}_{2}\
a.s.$ Then,

$\displaystyle E\big{[}Y_{m,a}Y_{j,b}\mid\mathcal{C}_{M}\big{]}$
$\displaystyle=E\big{[}\varepsilon_{m}\varepsilon_{j}x_{m,a}x_{j,b}\mid\mathcal{C}_{M}\big{]}$
$\displaystyle=x_{m,a}x_{j,b}\underbrace{E\big{[}\varepsilon_{m}\varepsilon_{j}\mid\mathcal{C}_{M}\big{]}}_{\leq\bar{C}_{1}}$
$\displaystyle\leq x_{m,a}x_{j,b}\bar{C}_{1}\quad a.s.,$

where $x_{m,a}$ represents the $a$-th element of $x_{M,m}$ and $x_{j,b}$ the
$b$-th element of $x_{M,j}$. Analogously, one obtains
$E\big{[}Y_{k,a}Y_{l,b}\mid\mathcal{C}_{M}\big{]}\leq
x_{k,a}x_{l,b}\bar{C}_{2}\ a.s.$ Once again, through the multiple application
of the Cauchy-Schwartz inequality, it follows that

$\displaystyle
E\big{[}E\big{[}Y_{m,a}Y_{j,b}\mid\mathcal{C}_{M}\big{]}E\big{[}Y_{k,a}Y_{l,b}\mid\mathcal{C}_{M}\big{]}\big{]}$
$\displaystyle\leq
E\big{[}\bar{C}_{1}\bar{C}_{2}x_{m,a}x_{j,b}x_{k,a}x_{l,b}\big{]}$
$\displaystyle\leq\big{(}E\big{[}(\bar{C}_{1}\bar{C}_{2})^{2}\big{]}\big{)}^{\frac{1}{2}}\big{(}E\big{[}(x_{m,a}x_{j,b}x_{k,a}x_{l,b})^{2}\big{]}\big{)}^{\frac{1}{2}}$
$\displaystyle\leq\big{(}E\big{[}(\bar{C}_{1})^{4}\big{]}\big{)}^{\frac{1}{4}}\big{(}E\big{[}(\bar{C}_{2})^{4}\big{]}\big{)}^{\frac{1}{4}}$
$\displaystyle\phantom{\leq\quad}\times\big{(}E\big{[}x_{m,a}^{8}\big{]}\big{)}^{\frac{1}{8}}\big{(}E\big{[}x_{j,b}^{8}\big{]}\big{)}^{\frac{1}{8}}\big{(}E\big{[}x_{l,a}^{8}\big{]}\big{)}^{\frac{1}{8}}\big{(}E\big{[}x_{l,b}^{8}\big{]}\big{)}^{\frac{1}{8}}.$

We note here that Assumption 3.3 ensures that there exists a nonnegative
finite constant $C_{m,a}$ such that $E\big{[}x_{l,a}^{8}\big{]}<C_{m,a}$, with
the same argument holding true for $x_{j,b}$, $x_{k,a}$ and $x_{l,b}$ as well.
Hence,

$\displaystyle
E\big{[}E\big{[}Y_{m,a}Y_{j,b}\mid\mathcal{C}_{M}\big{]}E\big{[}Y_{k,a}Y_{l,b}\mid\mathcal{C}_{M}\big{]}\big{]}$
$\displaystyle\leq\bar{C},$

where $\bar{C}$ is a nonnegative finite constant that is appropriately
defined.

Substituting this into the inequality above,

$\displaystyle E\big{[}(V_{a,b}^{c})^{2}\big{]}$
$\displaystyle\leq\frac{1}{M^{2}}\sum_{s\geq
0}\sum_{m\in\mathcal{M}_{N}}\sum_{j\in\mathcal{M}_{N}^{\partial}(m;s)}\sum_{t\geq
0}\sum_{k\in\mathcal{M}_{N}}\sum_{l\in\mathcal{M}_{N}^{\partial}(k;t)}\bar{C}$
$\displaystyle=\frac{\bar{C}}{M^{2}}\sum_{s\geq
0}\sum_{m\in\mathcal{M}_{N}}\sum_{j\in\mathcal{M}_{N}^{\partial}(m;s)}\sum_{t\geq
0}\sum_{k\in\mathcal{M}_{N}}\sum_{l\in\mathcal{M}_{N}^{\partial}(k;t)}1$
$\displaystyle=\frac{\bar{C}}{M^{2}}\sum_{s\geq 0}\sum_{(m,j,k,l)\in
H_{M}(s;b_{M})}1$ $\displaystyle=\frac{\bar{C}}{M^{2}}\sum_{s\geq
0}\underbrace{|H_{M}(s;b_{M})|}_{\leq 4Mc_{M}(s,b_{M};2)}$
$\displaystyle\leq\frac{\bar{C}}{M^{2}}\sum_{s\geq 0}4Mc_{M}(s,b_{M};2)$
$\displaystyle=4\bar{C}\underbrace{\frac{1}{M}\sum_{s\geq
0}c_{M}(s,b_{M};2)}_{\rightarrow 0}$ $\displaystyle\rightarrow 0,$

where the second inequality comes from Lemma C.1, and the last implication is
due to Assumption 4.2. Therefore we have shown that

$\displaystyle E\big{[}(V_{a,b}^{c}-V_{a,b})^{2}\big{]}\rightarrow 0.$

By repeating the same argument for each $a,b=1,\ldots,K$, it follows that

$\displaystyle E\big{[}\|V_{N,M}^{c}-V_{M,M}\|_{F}^{2}\big{]}\rightarrow 0.$

    * (b)

$\big{\|}V_{N,M}^{c}-V_{N,M}\big{\|}_{F}\stackrel{{\scriptstyle
p}}{{\rightarrow}}0$:
By the Chebyshev’s inequality and the result of part (a), we complete part (1)
as it follows that for any $\eta>0$,

$\displaystyle\mathrm{Pr}(\big{\|}V_{N,M}^{c}-V_{N,M}\big{\|}_{F}>\eta)$
$\displaystyle<\frac{1}{\eta^{2}}\underbrace{E\Big{[}\big{\|}V_{N,M}^{c}-V_{N,M}\big{\|}_{F}^{2}\Big{]}}_{\rightarrow
0}$ $\displaystyle\rightarrow 0.$

  * (2)

$\|\tilde{V}_{N,M}-V^{c}_{N,M}\|_{F}\stackrel{{\scriptstyle
p}}{{\rightarrow}}0$:
This immediately follows from applying Proposition 4.1 of Kojevnikov et al.
(2021)333333Notice that the definitions of $V_{N,M}$, $\hat{V}_{N,M}$,
$\tilde{V}_{N,M}$ and $V_{N,M}^{c}$ are slightly different from those used in
Proposition 4.1 of Kojevnikov et al. (2021). and the Dominated Convergence
Theorem in the Markov inequality: i.e.,

$\displaystyle\mathrm{Pr}\big{(}\|\tilde{V}_{N,M}-V^{c}_{N,M}\|_{F}\geq\eta\big{)}$
$\displaystyle\leq\frac{1}{\eta}E\big{[}\|\tilde{V}_{N,M}-V^{c}_{N,M}\|_{F}\big{]}$
$\displaystyle=\frac{1}{\eta}E\big{[}\underbrace{E\big{[}\|\tilde{V}_{N,M}-V^{c}_{N,M}\|_{F}\mid\mathcal{C}_{M}\big{]}}_{\stackrel{{\scriptstyle
a.s.}}{{\rightarrow}}0}\big{]}$ $\displaystyle\rightarrow 0,$

for any $\eta>0$.

  * (3)

$\|\hat{V}_{N,M}-\tilde{V}_{N,M}\|_{F}\stackrel{{\scriptstyle
p}}{{\rightarrow}}0$:
First, we have343434To lighten the notational burden, we drop the $M$
subscript from $\\{x_{M,m}\\}_{m\in\mathcal{M}_{N}}$ and
$\\{\varepsilon_{M,m}\\}_{m\in\mathcal{M}_{N}}$ in the rest of the proof.

$\displaystyle\big{\|}\hat{V}_{N,M}-\tilde{V}_{N,M}\big{\|}_{F}$
$\displaystyle=\Bigg{\|}\sum_{s\geq
0}\omega_{M}(s)\frac{1}{M}\sum_{m\in\mathcal{M}_{N}}\sum_{j\in\mathcal{M}_{N}^{\partial}(m;s)}\hat{\varepsilon}_{m}\hat{\varepsilon}_{j}x_{m}x_{j}^{\prime}-\sum_{s\geq
0}\omega_{M}(s)\frac{1}{M}\sum_{m\in\mathcal{M}_{N}}\sum_{j\in\mathcal{M}_{N}^{\partial}(m;s)}\varepsilon_{m}\varepsilon_{j}x_{m}x_{j}^{\prime}\Bigg{\|}_{F}$
$\displaystyle=\Bigg{\|}\sum_{s\geq
0}\omega_{M}(s)\frac{1}{M}\sum_{m\in\mathcal{M}_{N}}\sum_{j\in\mathcal{M}_{N}^{\partial}(m;s)}\big{(}\hat{\varepsilon}_{m}\hat{\varepsilon}_{j}-\varepsilon_{m}\varepsilon_{j}\big{)}x_{m}x_{j}^{\prime}\Bigg{\|}_{F}$
$\displaystyle\leq\Bigg{\|}\sum_{s\geq 0}\underbrace{|\omega_{M}(s)|}_{\leq
1}\frac{1}{M}\sum_{m\in\mathcal{M}_{N}}\sum_{j\in\mathcal{M}_{N}^{\partial}(m;s)}\big{(}\hat{\varepsilon}_{m}\hat{\varepsilon}_{j}-\varepsilon_{m}\varepsilon_{j}\big{)}x_{m}x_{j}^{\prime}\Bigg{\|}_{F}$
$\displaystyle\leq\Bigg{\|}\sum_{s\geq
0}\frac{1}{M}\sum_{m\in\mathcal{M}_{N}}\sum_{j\in\mathcal{M}_{N}^{\partial}(m;s)}\big{(}\hat{\varepsilon}_{m}\hat{\varepsilon}_{j}-\varepsilon_{m}\varepsilon_{j}\big{)}x_{m}x_{j}^{\prime}\Bigg{\|}_{F}$
$\displaystyle\leq\sum_{s\geq
0}\frac{1}{M}\sum_{m\in\mathcal{M}_{N}}\sum_{j\in\mathcal{M}_{N}^{\partial}(m;s)}\big{|}\hat{\varepsilon}_{m}\hat{\varepsilon}_{j}-\varepsilon_{m}\varepsilon_{j}\big{|}\big{\|}x_{m}x_{j}^{\prime}\big{\|}_{F}.$

Observe that, by definition, $\hat{\varepsilon}_{m}$ can be written as
$\hat{\varepsilon}_{m}=\varepsilon_{m}-x_{m}^{\prime}(\hat{\beta}-\beta)$.
Hence

$\displaystyle\hat{\varepsilon}_{m}\hat{\varepsilon}_{j}-\varepsilon_{m}\varepsilon_{j}=-\varepsilon_{m}(\hat{\beta}-\beta)^{\prime}x_{j}-x_{m}^{\prime}(\hat{\beta}-\beta)\varepsilon_{j}+x_{m}^{\prime}(\hat{\beta}-\beta)(\hat{\beta}-\beta)^{\prime}x_{j},$

so that by the triangular inequality,

$\displaystyle\big{|}\hat{\varepsilon}_{m}\hat{\varepsilon}_{j}-\varepsilon_{m}\varepsilon_{j}\big{|}$
$\displaystyle\leq\big{\|}\hat{\beta}-\beta\big{\|}_{2}\big{\|}x_{j}\big{\|}_{2}\big{|}\varepsilon_{m}\big{|}+\big{\|}\hat{\beta}-\beta\big{\|}_{2}\big{\|}x_{m}\big{\|}_{2}\big{|}\varepsilon_{j}\big{|}+\big{\|}\hat{\beta}-\beta\big{\|}_{2}^{2}\big{\|}x_{m}\big{\|}_{2}\big{\|}x_{j}\big{\|}_{2},$

for each $m,j\in\mathcal{M}_{N}$. Hence
$\|\hat{V}_{N,M}-\tilde{V}_{N,M}\|_{F}$ can be bounded as

$\displaystyle\big{\|}\hat{V}_{N,M}-\tilde{V}_{N,M}\big{\|}_{F}$
$\displaystyle\leq\big{\|}\hat{\beta}-\beta\big{\|}_{2}\frac{1}{M}\sum_{s\geq
0}\sum_{m\in\mathcal{M}_{N}}\sum_{j\in\mathcal{M}_{N}^{\partial}(m;s)}\|x_{m}\|_{2}\|x_{j}\|_{2}^{2}|\varepsilon_{m}|+\big{\|}\hat{\beta}-\beta\big{\|}_{2}\frac{1}{M}\sum_{s\geq
0}\sum_{m\in\mathcal{M}_{N}}\sum_{j\in\mathcal{M}_{N}^{\partial}(m;s)}\|x_{m}\|_{2}^{2}\|x_{j}\|_{2}|\varepsilon_{j}|$
$\displaystyle\qquad+\big{\|}\hat{\beta}-\beta\big{\|}_{2}^{2}\frac{1}{M}\sum_{s\geq
0}\sum_{m\in\mathcal{M}_{N}}\sum_{j\in\mathcal{M}_{N}^{\partial}(m;s)}\|x_{m}\|_{2}^{2}\|x_{j}\|_{2}^{2}.$

Denote

$\displaystyle R_{N,1}\coloneqq\frac{1}{M}\sum_{s\geq
0}\sum_{m\in\mathcal{M}_{N}}\sum_{j\in\mathcal{M}_{N}^{\partial}(m;s)}\|x_{m}\|_{2}\|x_{j}\|_{2}^{2}|\varepsilon_{m}|;$
$\displaystyle R_{N,2}\coloneqq\frac{1}{M}\sum_{s\geq
0}\sum_{m\in\mathcal{M}_{N}}\sum_{j\in\mathcal{M}_{N}^{\partial}(m;s)}\allowbreak\|x_{m}\|_{2}^{2}\|x_{j}\|_{2}|\varepsilon_{j}|;$
$\displaystyle R_{N,3}\coloneqq\frac{1}{M}\sum_{s\geq
0}\sum_{m\in\mathcal{M}_{N}}\sum_{j\in\mathcal{M}_{N}^{\partial}(m;s)}\|x_{m}\|_{2}^{2}\|x_{j}\|_{2}^{2}.$

Now, since by Theorem 3.1, $\|\hat{\beta}-\beta\|_{2}\stackrel{{\scriptstyle
p}}{{\rightarrow}}0$, and the application of the Continuous Mapping Theorem
yields $\|\hat{\beta}-\beta\|_{2}^{2}\stackrel{{\scriptstyle
p}}{{\rightarrow}}0$, it thus suffices to prove that each of $R_{N,1}$,
$R_{N,2}$ and $R_{N,3}$ converges in probability to a finite number. In
proving this, we follow a strategy employed in Aronow et al. (2015) and
Tabord-Meehan (2019).

First let us study the expectation of $R_{N,1}$. By applying the Cauchy-
Schwartz inequality repeatedly, we have that

$\displaystyle E\big{[}R_{N,1}\big{]}\leq\frac{1}{M}\sum_{s\geq
0}\sum_{m\in\mathcal{M}_{N}}\sum_{j\in\mathcal{M}_{N}^{\partial}(m;s)}\big{(}\big{(}E\big{[}\|x_{m}\|_{2}^{2}\big{]}\big{)}^{\frac{1}{2}}\big{(}E\big{[}\|x_{j}\|_{2}^{8}\big{]}\big{)}^{\frac{1}{2}}\big{)}^{\frac{1}{2}}\big{(}E\big{[}E\big{[}|\varepsilon_{m}|^{2}\mid\mathcal{C}_{M}\big{]}\big{]}\big{)}^{\frac{1}{2}}.$

Here, in light of Assumption 4.1, there exists an a.s.-bounded function
$C_{1}$ such that $C_{1}=\sup_{N\geq
1}\max_{m\in\mathcal{M}_{N}}E\big{[}|\varepsilon_{m}|^{2}\mid\mathcal{C}_{M}\big{]}$,
and moreover by Assumption 3.3, there exists a nonnegative finite number
$C_{2}>0$ such that $C_{2}=\sup_{N\geq
1}\max_{m\in\mathcal{M}_{N}}E\big{[}\|x_{m}\|_{2}^{8}\big{]}$. With a slight
abuse of notation, we have for every $N>0$

$\displaystyle E\left[R_{N,1}\right]$ $\displaystyle\leq\frac{1}{M}\sum_{s\geq
0}\sum_{m\in\mathcal{M}_{N}}\sum_{j\in\mathcal{M}_{N}^{\partial}(m;s)}C_{1}C_{2}$
$\displaystyle=\frac{C_{1}C_{2}}{M}\sum_{s\geq
0}\sum_{m\in\mathcal{M}_{N}}\sum_{j\in\mathcal{M}_{N}^{\partial}(m;s)}1$
$\displaystyle=C_{1}C_{2}\sum_{s\geq
0}\underbrace{\frac{1}{M}\sum_{m\in\mathcal{M}_{N}}|\mathcal{M}_{N}^{\partial}(m;s)|}_{\delta_{M}^{\partial}(s;1)}$
$\displaystyle=C_{1}C_{2}\underbrace{\sum_{s\geq
0}\delta_{M}^{\partial}(s;1)}_{<\infty}$ $\displaystyle<C,$

for some constant $C\in(0,\infty)$, where the last inequality is because of
Assumption 4.2.

Next let us study the variance of $R_{N,1}$. It suffices to show that
$E\big{[}R_{N,1}^{2}\big{]}\rightarrow 0$, By the Cauchy-Schwartz inequality,
it holds that

$\displaystyle E\big{[}R_{N,1}^{2}\big{]}$
$\displaystyle=E\Big{[}\frac{1}{M^{2}}\sum_{s\geq
0}\sum_{m\in\mathcal{M}_{N}}\sum_{j\in\mathcal{M}_{N}^{\partial}(m;s)}\sum_{t\geq
0}\sum_{k\in\mathcal{M}_{N}}\sum_{l\in\mathcal{M}_{N}^{\partial}(k;t)}\|x_{m}\|_{2}\|x_{j}\|_{2}^{2}\|x_{k}\|_{2}\|x_{l}\|_{2}^{2}|\varepsilon_{m}||\varepsilon_{k}|\Big{]}$
$\displaystyle\leq\frac{1}{M^{2}}\sum_{s\geq
0}\sum_{m\in\mathcal{M}_{N}}\sum_{j\in\mathcal{M}_{N}^{\partial}(m;s)}\sum_{t\geq
0}\sum_{k\in\mathcal{M}_{N}}\sum_{l\in\mathcal{M}_{N}^{\partial}(k;t)}\big{(}E\big{[}\|x_{m}\|_{2}^{2}\|x_{j}\|_{2}^{4}\|x_{k}\|_{2}^{2}\|x_{l}\|_{2}^{4}\big{]}\big{)}^{\frac{1}{2}}\big{(}E\big{[}|\varepsilon_{m}|^{2}|\varepsilon_{k}|^{2}\big{]}\big{)}^{\frac{1}{2}}.$

Here, by Assumption 3.3 and the Cauchy-Schwartz inequality, there exists a
nonnegative finite constant $C_{3}>0$ such that $C_{3}=\sup_{N\leq
1}\max_{m,j,k,l\in\mathcal{M}_{N}}E\big{[}\|x_{m}\|_{2}^{2}\|x_{j}\|_{2}^{4}\|x_{k}\|_{2}^{2}\|x_{l}\|_{2}^{4}\big{]}$.
Then, with a slight abuse of notation in writing $C_{3}^{\frac{1}{2}}$ as
$C_{3}$, we have

$\displaystyle E\big{[}R_{N,1}^{2}\big{]}$
$\displaystyle=\frac{C_{3}}{M^{2}}\sum_{s\geq
0}\sum_{m\in\mathcal{M}_{N}}\allowbreak\sum_{j\in\mathcal{M}_{N}^{\partial}(m;s)}\sum_{t\geq
0}\sum_{k\in\mathcal{M}_{N}}\sum_{l\in\mathcal{M}_{N}^{\partial}(k;t)}\big{(}E\big{[}|\varepsilon_{m}|^{2}|\varepsilon_{k}|^{2}\big{]}\big{)}^{\frac{1}{2}}$
$\displaystyle=\frac{C_{3}}{M^{2}}\sum_{s\geq 0}\sum_{(m,j,k,l)\in
H_{M}(s;b_{M})}\allowbreak\big{(}E\big{[}E\big{[}|\varepsilon_{m}|^{2}|\varepsilon_{k}|^{2}\mid\mathcal{C}_{M}\big{]}\big{]}\big{)}^{\frac{1}{2}}.$

Corollary A.2 of Kojevnikov et al. (2021) shows that there exists a
nonnegative finite constant $C_{4}$ such that
$E\big{[}|\varepsilon_{m}|^{2}|\varepsilon_{k}|^{2}\mid\mathcal{C}_{M}\big{]}\leq
C_{4}\bar{\theta}\theta_{M,s}^{1-\frac{4}{p}}$, where
$\bar{\theta}\coloneqq\sup_{M\geq 1}\max_{s\geq 1}\theta_{M,s}$. Upon applying
Lemma C.1 from the Appendix, we obtain

$\displaystyle
E\big{[}R_{N,1}^{2}\big{]}\leq\frac{C_{3}C_{4}^{\prime}}{M^{2}}\sum_{s\geq
0}\big{(}E\big{[}\theta_{M,s}^{1-\frac{4}{p}}\big{]}\big{)}^{\frac{1}{2}}4Mc_{M}(s,b_{M};2)=\frac{4C_{3}C_{4}^{\prime\prime}}{M}\sum_{s\geq
0}c_{M}(s,b_{M};2)\rightarrow 0,$

where we apply Assumption 4.2 for the last implication, and $C_{4}^{\prime}$
and $C_{4}^{\prime\prime}$ are nonnegative finite constants defined
appropriately. Hence we have shown that $R_{N,1}$ converges to a finite
constant.

The proof of $R_{N,2}$ is analogous.

It remains to show that $R_{N,3}$ converges in probability to a finite
constant. Let us first study the expectation of $R_{N,3}$. Observe that

$\displaystyle E\big{[}R_{N,3}\big{]}=\frac{1}{M}\sum_{s\geq
0}\sum_{m\in\mathcal{M}_{N}}\sum_{j\in\mathcal{M}_{N}^{\partial}(m;s)}E\big{[}\|x_{m}\|_{2}^{2}\|x_{j}\|_{2}^{2}\big{]}.$

By Assumption 3.3, there exists a nonnegative finite number $C_{5}>0$ such
that $C_{5}=\sup_{N\geq
1}\max_{m\in\mathcal{M}_{N}}E\big{[}\|x_{m}\|_{2}^{2}\|x_{j}\|_{2}^{2}\big{]}$.
Hence for every $N>0$,

$\displaystyle E\big{[}R_{N,3}\big{]}$ $\displaystyle=\frac{1}{M}\sum_{s\geq
0}\sum_{m\in\mathcal{M}_{N}}\sum_{j\in\mathcal{M}_{N}^{\partial}(m;s)}\underbrace{E\Big{[}\|x_{m}\|_{2}^{2}\|x_{j}\|_{2}^{2}\Big{]}}_{\leq
C_{5}}$ $\displaystyle\leq\frac{1}{M}\sum_{s\geq
0}\sum_{m\in\mathcal{M}_{N}}\sum_{j\in\mathcal{M}_{N}^{\partial}(m;s)}C_{5}$
$\displaystyle=\frac{C_{5}}{M}\sum_{s\geq
0}\sum_{m\in\mathcal{M}_{N}}\sum_{j\in\mathcal{M}_{N}^{\partial}(m;s)}1$
$\displaystyle=C_{5}\sum_{s\geq
0}\underbrace{\frac{1}{M}\sum_{m\in\mathcal{M}_{N}}|\mathcal{M}_{N}^{\partial}(m;s)|}_{\delta_{M}^{\partial}(s;1)}$
$\displaystyle=C_{5}\underbrace{\sum_{s\geq
0}\delta_{M}^{\partial}(s;1)}_{<\infty}$ $\displaystyle<C,$

where we apply Assumption 4.2 in the last implication and a constant
$C\in(0,\infty)$ is appropriately defined.

Next let us consider the variance of $R_{M,3}$:

$\displaystyle E\big{[}R_{N,3}^{2}\big{]}$
$\displaystyle=\frac{1}{M^{2}}\sum_{s\geq
0}\sum_{m\in\mathcal{M}_{N}}\sum_{j\in\mathcal{M}_{N}^{\partial}(m;s)}\sum_{t\geq
0}\sum_{k\in\mathcal{M}_{N}}\sum_{l\in\mathcal{M}_{N}^{\partial}(k;t)}E\big{[}\|x_{m}\|_{2}^{2}\|x_{j}\|_{2}^{2}\|x_{k}\|_{2}^{2}\|x_{l}\|_{2}^{2}\big{]}.$

Once again, Assumption 3.3 and the Cauchy-Schwartz inequality imply that there
exists a nonnegative finite number $C_{6}>0$ such that $C_{6}=\sup_{N\geq
1}\max_{m,j,k,l\in\mathcal{M}_{N}}\allowbreak
E\big{[}\|x_{m}\|_{2}^{2}\|x_{j}\|_{2}^{2}\|x_{k}\|_{2}^{2}\|x_{l}\|_{2}^{2}\big{]}$.
Then by Lemma C.1,

$\displaystyle E\big{[}R_{N,3}^{2}\big{]}\leq\frac{1}{M^{2}}\sum_{s\geq
0}\sum_{m\in\mathcal{M}_{N}}\sum_{j\in\mathcal{M}_{N}^{\partial}(m;s)}\sum_{t\geq
0}\sum_{k\in\mathcal{M}_{N}}\sum_{l\in\mathcal{M}_{N}^{\partial}(k;t)}C_{6}=\frac{4C_{6}}{M}\sum_{s\geq
0}c_{M}(s,b_{M};2)\rightarrow 0,$

where the last implication is a consequence of Assumption 4.2 (ii).

Therefore we have shown that

$\displaystyle\big{\|}\hat{V}_{N,M}-\tilde{V}_{N,M}\big{\|}_{F}$
$\displaystyle\leq\big{\|}\hat{\beta}-\beta\big{\|}_{2}\underbrace{\frac{1}{M}\sum_{s\geq
0}\sum_{m\in\mathcal{M}_{N}}\sum_{j\in\mathcal{M}_{N}^{\partial}(m;s)}\|x_{m}\|_{2}\|x_{j}\|_{2}^{2}|\varepsilon_{m}|}_{R_{M,1}}+\big{\|}\hat{\beta}-\beta\big{\|}_{2}\underbrace{\frac{1}{M}\sum_{s\geq
0}\sum_{m\in\mathcal{M}_{N}}\sum_{j\in\mathcal{M}_{N}^{\partial}(m;s)}\|x_{m}\|_{2}^{2}\|x_{j}\|_{2}|\varepsilon_{j}|}_{R_{M,2}}$
$\displaystyle\qquad+\big{\|}\hat{\beta}-\beta\big{\|}_{2}^{2}\underbrace{\frac{1}{M}\sum_{s\geq
0}\sum_{m\in\mathcal{M}_{N}}\sum_{j\in\mathcal{M}_{N}^{\partial}(m;s)}\|x_{m}\|_{2}^{2}\|x_{j}\|_{2}^{2}}_{R_{M,3}}$
$\displaystyle=\underbrace{\big{\|}\hat{\beta}-\beta\big{\|}_{2}}_{\stackrel{{\scriptstyle
p}}{{\rightarrow}}0}\underbrace{R_{M,1}}_{<\infty}+\underbrace{\big{\|}\hat{\beta}-\beta\big{\|}_{2}}_{\stackrel{{\scriptstyle