Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
95,767 characters · 14 sections · 114 citation commands
Empirical Challenges with Peers-of-Peers Instruments in the Linear-In-Means Model
Keywords: Social Networks, Weak Instruments, Peer Effects, Identification
\doublespacing
Humans are inherently social beings, frequently interacting in groups and affecting the behavior of their friends and neighbors. Thus, it comes as no surprise that the study of peer effects has become extremely popular in empirical research in Economics and Social Sciences more generally\footnote{Examples include peer effects in education, e.g., Sacerdote2001, calvo2009peer, worker productivity and labor markets (e.g., MasMoretti2009, caria2024village), Finance (e.g., LoughranSchultz2004 on the impact of IPOs on competitors); development and public goods (e.g., acemoglu2015), among many others. See bramoulle2020peer for a recent survey.}, especially with the emergence of high quality data on social interactions and advances in network statistics.
The most popular model of peer effects is arguably the linear-in-means model illustrated in equation (ref). In this model, one's outcome ($Y_i$) depends linearly on the mean outcome across $i$'s group, denoted $\bar{Y}_i$. The outcome may also depend on the exogenous characteristics of one's self $X_i$, the group itself ($P(i)$), or the average characteristics of the group with $\varepsilon_{i}$ capturing the unobserved error\footnote{For example, in calvo2009peer, the outcome (grades of student $i$) can depend on the average of its peers' grades ($\beta$), their own personal traits/parental education ($X_i$), the average group characteristics through $\delta$ and an unobserved error term ($\varepsilon_{i}$).}:
where $\bar{Y}_i = \frac{\sum_{j\in P_i}Y_{j}}{n_i}$ and $\bar{X}_i = \frac{\sum_{j\in P_i}X_{j}}{n_i}$, and $n_i$ is the number of $i$'s peers. Even if the error term, $\varepsilon_{i}$, is exogenous to $X_i$ and to the peer groups, endogeneity still arises in this model due to the simultaneous determination of behavior within the groups: the reflection problem (manski1993identification). After all, an increase in $\varepsilon_i$ affects one's $Y_i$, which then affects others' $Y_j$, leading to correlation between the average group outcome and the error.
Empirical papers typically solved this challenge by exploiting additional information such as external variables for instruments or randomization.\footnote{For example, Sacerdote2001 exploits randomization of individuals to groups, brock2001interactions exploits a specific block structure of groups, while other works use Instrumental Variables (IV) based on historical (e.g., acemoglu2015) or other external restrictions (ioannides2003neighbourhood and durlauf2008understanding)} Unfortunately, these solutions are unavailable for many settings. Yet, in a seminal contribution, bramoulle2009identification showed that $(\beta, \boldsymbol{\delta}', \boldsymbol{\gamma}')'$ could be identified using only the model above, the existing $X_i$ and the network structure itself, $\mathbf{G}$. They proposed using instruments based on the characteristics of friends-of-friends or higher-order connections (i.e. $\mathbf{G}^k \mathbf{X}$, where $\mathbf{X}$ is the matrix stacking $X_i'$).\footnote{Such instruments are valid because they are excluded from (ref) and $\mathbf{X}$ and $\mathbf{G}$ are exogenous. They are relevant because, in general, multiplying $\mathbf{G}_i$ on both sides of (ref) implies that $\bar{Y}_i$ is a function of $G_i^2 X_i$, where $\mathbf{G}_i, \mathbf{G}_i^2$ represent $i$'s friends and friends-of-friends, respectively.} This solution provided an easy to implement identification strategy, a natural estimator and straightforward inference using random variables readily available to the researcher.
In this paper, we first show that "friends-of-friends" instrumental variables are not a panacea for linear-in-means applications in Economics. In particular, we show that, in many empirical settings, these instruments are weak and/or lead to arbitrary estimates due to ill-defined first-stage estimands. In a nutshell, this happens when the underlying network is very sparse or very dense, so that the friends-of-friends network becomes arbitrarily similar to the original network. This results in instruments with low variance and potentially low covariance with the endogenous variable. Such sparse networks are also a prevalent issue in empirical work. Table (ref) shows the degree distribution statistics for each of the networks $(\mathbf{G})$ as well as the squared counterparts ($\mathbf{G}^2$) used to form the IV in salient examples in political economy (alumni networks in the U.S. Congress, battaglini2018) and development economics (network of allies/enemies in the Second Congo War, Konig2017). The networks are indeed very sparse: the modal degree for $\mathbf{G}$ is 0. But this is also true for $\mathbf{G}^2$! Such networks can then be visualized in Figures (ref)-(ref) and we revisit them further below.
Second, we characterize regimes where these challenges may arise focusing on Erdos-Renyi random graphs (erdds1959random), used extensively in both theoretical economics and econometrics (see jackson10, mele17, campbell24 for examples and discussions). Using tools from random graph theory we prove conditions on network sparsity/density, captured by the average degree $d_n$, that induce either problem as sample size grows. We formally show that the proposal of bramoulle2009identification works well when the networks are neither too sparse nor too dense. We analyze two scenarios: one in which the adjacency matrix is left unscaled and another in which it is scaled to control its spectral norm. In the unscaled model, when the average degree ($d_n$) decreases with $n$, the variance of the instrument collapses faster than its covariance with the endogenous regressor. Then, the first stage becomes asymptotically ill-defined. On the other extreme, when $d_n$ grows with $n$, $\mathbf G_n^{2}$ becomes asymptotically collinear with $\mathbf G_n$. Thus, the instrument adds little independent variation and first-stage relevance vanishes. Between these extremes, when the average degree is bounded, the strength of the first-stage depends on how quickly sampling noise dissipates which, in turn depends on network dependence and the assumed variance structure. Meanwhile, we find that appropriate scaling often stabilizes the spectrum and aligns the relevant growth rates, preventing purely mechanical divergence or collapse of the first-stage estimand. However, scaling is not a silver bullet: we find that weak identification can still arise in finite samples or near the stability boundary, particularly when spillover effects are large.
We then propose an adapted weak-IV robust inference, adapting an anderson1949estimation test with a network-spatial variance estimator of kojevnikov2021limit, and show that it provides asymptotically valid inference (under weak-IV asymptotics). This approach explicitly takes the network cross-sectional dependence into account and, in Monte Carlo simulations, is shown to perform well. However, we show that using a variance estimator under homoscedasticity is both simpler and performs very well in this setting because, for very sparse networks, network spillovers are limited.
We show that accounting for the very sparse nature of some economics networks (and their effects on estimation and inference) can lead to different conclusions in empirical examples. In particular, we revisit the setting of Konig2017 who studied the effects of allied (or enemy) networks across ethnicities in Africa and their effects on conflict. To account for endogeneity, they propose instruments that use the network structure of such linkages. Due to the sparsity of the network shown in Table (ref) and in Figures (ref) and (ref), the instrument was discussed as weak. Our proposed inference finds an associated confidence interval for the parameter of interest ($\beta$ above) that is much larger, and includes 0.
Finally, we conclude by discussing that our insights extend beyond the linear-in-means model, to other linear regression models with network-based instruments. Thus, close attention to network structure and its growth with sample size must be considered when implementing such instruments.
The rest of the paper is organized as follows: Section (ref) provides a brief review of literature, while Section (ref) contains the main theoretical results of our paper. Section (ref) provides results on our Monte Carlo simulations, while Section (ref) provides the empirical applications. We conclude in Section (ref).
There has been a steady growth in the literature dealing with econometric issues related to peer effects and social interactions. The issue of identification in a linear-in-means context was first studied by manski1993identification which spurred a large literature (see brock2001interactions, durlauf2004neighborhood, blume2005identifying for an extensive survey and bramoulle2020peer for a more recent one) looking specifically at the incidence of the "reflection problem". Initially, empirical research posed solutions to the problem of endogeneity that relied upon extra assumptions such as gaviria2001school assuming existence of only one type of social effect. Other solutions relied on an IV approach, such as ioannides2003neighbourhood and durlauf2008understanding which uses group analogues of individual characteristics satisfying an exclusion restriction. Some researchers used randomization as an appropriate identification strategy like field2016friendship. Meanwhile, papers like acemoglu2015 used historical variables that were exogenous to the network. All of these solutions, however, are case specific and not easily generalizable. In the absence of such identification strategies, bramoulle2009identification suggested using the network structure itself to generate valid instruments using friends-of-friends' characteristics or other higher powers of the adjacency matrix to construct instruments, as discussed above.
Our results speak to the literature on weak instruments (see andrews2019weak for a detailed survey). However, these papers - e.g., since dufour97 and staigerstock97 - typically consider the classical IV framework, rather than the linear-in-means model. One of our main theoretical contributions is to show how the latter can be transformed into a similar system, how they can be used for instruments, and the different assumptions for inference. This is where one of our main theoretical contributions lie. Then, we show how one can use inference that is robust to weak instruments (e.g., anderson1949estimation and moreira03) despite the different cross-sectional due to spillovers inherent to peer effects. This requires adapting results on consistent variance estimation with spillovers (e.g., kojevnikov2021limit or conley1999gmm).
A recent paper by Wang2025 also studies weak identification in peer-effect estimation, focusing on linear-in-means models. Their analysis is conducted under a near-degree-regularity assumption and characterizes the bias and convergence rates of Ordinary Least Squares (OLS) and Two Stage Least Squares (TSLS) estimators with i.i.d.\ covariates. In settings where neighborhoods become asymptotically identical as the network grows, they show that weak identification arises from asymptotic collinearity between network aggregates. In contrast, we study weak identification in network-based instrumental variables settings, where network sparsity and heterogeneous degree growth are key in determining strength of identification. To address collinearity under near-regularity, Wang2025 recommend working with a linear-in-sums specification in which the adjacency matrix is scaled to align the rates of different network regressors, rather than row-normalized. Our theoretical results reinforce and extend this insight by showing that such scaling is also essential in sparse and heterogeneous networks, where it stabilizes the first stage and prevents weak-instrument pathologies. Finally, unlike Wang2025, our analysis explicitly focuses on weak-IV–robust inference, which is central to our contribution.
We use the extended linear-in-means model of bramoulle2009identification where $N$-agents interact over an exogenously given network with the $N \times N$ adjacency matrix $\mathbf{G}$. Each element of this adjacency matrix is given by $g_{ij}$, where:
We consider the matrix version of the structural model in equation (ref) and, for simplicity, consider the case without group fixed-effects (called correlated effects), as they can be differenced out (see bramoulle2009identification).
where $Y$ and $\varepsilon$ is the $N \times 1$ vector of outcomes and errors, respectively, and $\mathbf{X}$ the $N \times d$ matrix of observable characteristics, and we have separated the constant from $\mathbf{X}$. $\boldsymbol{\gamma}$ and $\boldsymbol{\delta}$ are the $d \times 1$ vector of coefficients associated with $\mathbf{X}$ and $\mathbf{GX}$ respectively. We make the following standard assumptions:
Assumption (ref) (i) and (ii) are standard to guarantee the behavior of the model, including the invertibility and stability of the system. Condition (ii) ensures that ($\mathbf{I} - \beta \mathbf{G}$) is invertible and can be expanded into an infinite matrix (Neumann) series. Condition (iii) assumes that not all nodes are isolated, which is crucial for identification. Then, (iv) assumes strict exogeneity of individual characteristics and the network respectively. Now, we can write the reduced-form of ((ref)) as:
Under Assumptions (ref) (i)-(ii),we can expand $(\mathbf{I} - \beta \mathbf{G})^{-1}$ into a Neumann series:
Then, from the strict exogeneity assumption, we get:
The model is said to be identified if $\boldsymbol{\theta} = (\alpha, \beta, \boldsymbol{\gamma}', \boldsymbol{\delta}')'$ is identified.\footnote{We assume there is a super-population of exogenous networks from which the sample $\mathbf{G}$ is drawn, thereby defining identification relative to this super-population, and the DGP given in equation (ref).} The reflection problem discussed by manski1993identification is evident from equation (ref) where individual outcomes are affected by the respective expected group outcome which in itself is impacted by the former. bramoulle2009identification suggested using the exogenously given network structure to construct valid instruments. The main identification result of their paper states that, as long as attributes/characteristics of neighbors have some direct or indirect effect, i.e. $( \beta \boldsymbol{\gamma} + \boldsymbol{\delta}) \neq \mathbf{0}$ and $\mathbf{I}$, $\mathbf{G}$ and its higher powers are not linearly dependent, then $\mathbf{G}^k \mathbf{X}$ for $k \geq 2 $ can be used as valid instruments for $\mathbf{GY}$. This is summarized in their Proposition 1, rewritten for convenience below.
Hence, if the first stage is given by
the structural equation (ref) can be re-written as,
where $\eta = \epsilon + \beta \Tilde{\epsilon}$ and $\boldsymbol{\xi} = \beta \boldsymbol{\pi}$. Thus, the endogenous effect $\beta$ represents the proportionality constant linking the coefficient on $\mathbf{G}^{2} \mathbf{X}$ in the first-stage regression to that in the structural equation (ref) which can be estimated using an appropriate estimator $\hat{\boldsymbol{\xi}} = \hat{\beta} \hat{ \boldsymbol{\pi}}$.
Proposition (ref) shows that, for non-trivial combinations of parameters, identification hinges on the informational content of friends-of-friends' networks ($\mathbf{G}^2$) relative to the network $\mathbf{G}$ and to a constant. This suggests that, even when the exclusion restriction and rank conditions hold in principle, the effectiveness of $\mathbf{G}^2 \mathbf X$ as an instrument depends critically on whether it introduces sufficient independent variation beyond $\mathbf{GX}$. In the next section, we formalize the conditions for this to arise in a salient class of random graph models. Here, we outline the main ideas that apply more generally.
First, we demonstrate that the key assumptions in Proposition (ref), namely the linear independence of $\mathbf{G}$ and $\mathbf{G}^2$ as well as $\mathbf{G}^2$ being nonzero (see footnote 23 of bramoulle2009identification), are likely violated (or close to violated) in many empirical settings. Table (ref) already provides evidence of this problem arising in various empirically observed networks. This is also easily observed when the graphs themselves are plotted in Figures (ref), (ref) and (ref). These examples illustrate situations where $\mathbf{G}^2$ is extremely sparse. On the other hand, Figure (ref) showcases a much denser network where higher-order links become very close to a completed network. Our second main contribution is to characterize the behavior of the network-based instruments in different scenarios.
In particular, the two extremes mentioned above lead to distinct implications for identification and inference in these contexts. When the average degree tends to zero, the squared network $\mathbf{G}^2$ collapses to the zero matrix at a faster rate than $\mathbf{G}$, with linear independence being violated in the limit as $n \to \infty$. Hence, the variance of the instrument $\mathbf{G}^2 \mathbf{X}$ tends to zero faster than its covariance with the endogenous variable $\mathbf{GY}$. We show below that, rather than being weakly correlated with the endogenous variable, the limit estimand of the first stage regression becomes undefined as a very small covariance is normalized by an even smaller variance. Empirically this manifests as arbitrarily large and potentially unstable first-stage estimates of either sign. This is not identification failure in the usual sense, but rather asymptotic degeneracy of the instrument itself, where the limit estimand is ill-defined. At the opposite extreme, when the network becomes very dense, the second-order matrix satisfies $\mathbf{G}^2 \approx c\,\mathbf{G}$ for some scalar $c>0$. Thus, $\mathbf{G}^2\mathbf X$ introduces no new variation beyond $\mathbf{GX}$. Consequently, the population first stage coefficient tends to zero and the model becomes unidentified. This is precisely the identification failure highlighted in bramoulle2009identification: in dense networks, first- and second-order peer effects cannot be distinguished because the relevant regressors become collinear. This is further explored by Wang2025 in near-regular graphs where neighborhoods become asymptotically identical and hence higher order network structures do not provide additional independent variation beyond $\mathbf G$.
Between these two extremes lie two intermediate scenarios. First, when $\mathbf{G}^2$ is sparse but not degenerate and $\mathbf{G}^2 \mathbf X$ remains a valid instrument but may be weak. In the second, the network is moderately dense and $\mathbf{G}^2$ is not yet collinear with $\mathbf{G}$, so that $\mathbf{G}^2 X$ provides substantial independent variation and yields a strong first stage. Weakness in the sparse-but-non-degenerate case arises when $\mathrm{Cov}(\mathbf{G} \mathbf Y,\mathbf{G}^k \mathbf X)\quad (k\ge 2)$ is small relative to the variance of $\mathbf{G}^k \mathbf X$. This behavior is driven by sparsity: when the network has low average degree and high degree heterogeneity, $\mathbf{G}^k \mathbf X$ exhibits very low (but finite) variance, and its covariance with the endogenous variable $\mathbf{G} \mathbf Y$ decays to zero.
Following stock00 and andrews12, we formalize weak identification through the population first-stage coefficient $\boldsymbol{\pi}$ in (ref). As emphasized by andrews2019weak, non-standard asymptotic behavior of IV estimators arises when $\boldsymbol{\pi}$ is small relative to the sampling variability of $\hat{\boldsymbol{\pi}}$. In the standard i.i.d. case studied in staigerstock97, the sampling variability of $\hat{\boldsymbol{\pi}}$ is of order $n^{-1/2}$. Hence, instruments are weak when $\boldsymbol{\pi}$ is local to zero at the same rate, i.e.\ when $\boldsymbol{\pi} = O(n^{-1/2})$.\footnote{Formally, $\boldsymbol{\pi} = O(n^{-1/2})$ implies that there exists a fixed (non-random) matrix $\mathbf{C}$ such that $\sqrt{n}\,\boldsymbol{\pi} \leq \mathbf{C}$ as $n \to \infty$.} In models with network dependence, however, the sampling variability of $\hat{\boldsymbol{\pi}}$ depends on the strength and structure of cross-sectional dependence induced by the network. As a result, the rate at which the first-stage estimator concentrates around its population value need not be $\sqrt{n}$ and may vary with the size and connectivity of the network (see kojevnikov2021limit and Wang2025 among others). Thus, for our definition of weak instruments with network-based instruments, we introduce $k(n)$ below that captures the true rate of sampling variance decay with network-dependence.
To do so, we must introduce additional notation. Formally, let $Z_n=\mathbf{G}^{k}\mathbf{X}$ with $k\ge2$, where the subscript $n$ emphasizes the dependence of the instrument on the size and structure of the network. Let $W = (\iota,\mathbf{X},\mathbf{G}\mathbf{X})$ and define the residualized instrument $\tilde{Z}_n := M_W Z_n$, where $M_W=I-W(W'W)^{-1}W'$. Let $\widetilde{\mathbf{G} \mathbf Y}:=M_W(\mathbf{G} \mathbf Y)$ denote the corresponding residualized endogenous regressor. We let $\boldsymbol{\pi}_n$ denote the population first-stage coefficient, which is allowed to depend on $n$ and is given by the Frisch--Waugh--Lovell representation \[ \boldsymbol{\pi}_n := \left(Var(\tilde Z_n)\right)^{-1} Cov\!\left(\widetilde{\mathbf{G} \mathbf Y},\,\tilde Z_n\right), \] where $Var(\tilde Z_n)=\mathbb{E}[\tilde Z_n'\tilde Z_n]/n$ and $Cov(\widetilde{\mathbf{G}Y},\tilde Z_n) =\mathbb{E}[\tilde Z_n'\widetilde{\mathbf{G}Y}]/n$.\footnote{This dependence arises because both the strength of the instrument and the sampling variability of its estimator are functions of the network’s size and connectivity, so asymptotic behavior is governed by the sequence of networks $\{\mathbf{G}_n\}_{n\ge1}$ rather than by sample size alone.}
The sampling variability of \(\hat{\boldsymbol{\pi}}_n\) is driven by the first-stage error $\tilde\epsilon_n := \widetilde{\mathbf{G}Y}-\tilde Z_n\boldsymbol{\pi}_n$. Conditional on \(\mathbf{G}\), the variance of the first-stage estimator is therefore \[ Var(\hat{\boldsymbol{\pi}}_n\mid\mathbf{G}) = \left(Var(\tilde Z_n)\right)^{-1} \boldsymbol{\Omega}_{\pi,n} \left(Var(\tilde Z_n)\right)^{-1}, \] where \[ \boldsymbol{\Omega}_{\pi,n} := Var\!\left( \frac{1}{\sqrt{n}}\tilde Z_n'\tilde\epsilon_n \;\middle|\; \mathbf{G} \right). \] Accordingly, we scale the first-stage parameter by the network-dependent information index \[ k(n) \;\asymp\; \left\| \,Var(\hat{\boldsymbol{\pi}}_n\mid\mathbf{G})^{-1/2} \right\|, \] implying that $k(n)$ grows at the same rate as the conditional variance of the estimator.\footnote{This variance normalization coincides with that used by kojevnikov2021limit to establish a central limit theorem for network-dependent quadratic forms. Under conditional homoskedasticity of \(\tilde\epsilon_n\) given \(\mathbf{G}\), \(\boldsymbol{\Omega}_{\pi,n}\) is proportional to \(\mathbb{E}[\tilde Z_n'\tilde Z_n\mid\mathbf{G}]\), implying that the information content of the first stage is governed—up to constants—by the Frobenius norm \(\|\mathbf{G}^k\|_F^2\). This scaling is also adopted in Wang2025, who shows that it yields stable Gaussian limits under both sparse and dense network sequences.} Finally, network based instruments are considered weak if they satisfy the following definition.
Equivalently, this implies that instruments are weak when $k(n) \boldsymbol{\pi}_n=O(1)$ as $n\to\infty$.
When $\mathbf{X}$ is uni-dimensional, the population first-stage coefficient reduces to \[ \pi_n = \frac{\frac{1}{n}Cov(\widetilde{\mathbf{G}Y},\tilde Z_n)} {\frac{1}{n}Var(\tilde Z_n)}, \] with \[ \hat\pi_n-\pi_n = \left(\tilde Z_n'\tilde Z_n\right)^{-1}\tilde Z_n'\tilde\epsilon_n. \] Furthermore, define \[ \sigma_{\pi,n}^2 := Var(\hat\pi_n\mid\mathbf{G}) = \left(\frac{1}{n}\tilde Z_n'\tilde Z_n\right)^{-2} Var\!\left( \frac{1}{n}\tilde Z_n'\tilde\epsilon_n \;\middle|\; \mathbf{G} \right). \] The information index then simplifies to \[ k(n) \asymp \left(\frac{1}{\sigma_{\pi,n}}\right), \] with Definition (ref) reducing to
This definition mirrors the local-to-zero framework of staigerstock97. Sparsity therefore leads either to degeneracy (when $\mathbf{G}^2\to 0$ too quickly) or to weak identification (when $\mathrm{Cov}(\mathbf{G} \mathbf Y,\mathbf{G}^k \mathbf{X})$ decreases faster than $\mathrm{Var}(\mathbf{G}^k \mathbf X)$ relative to the sampling noise. Intuitively, in a sparse network the friends-of-friends matrix $\mathbf{G}^2$ (and other higher order terms) contains many zeros and a few nodes with disproportionately large reach. As a result, $\mathbf{G}^k \mathbf X$ tends to have low covariance with $\mathbf{G} \mathbf Y$ in the sparse case.
Figure (ref) reconciles the observations documented in Table (ref) with the weak-instrument definition stated in Definition (ref) which depends on the first-stage covariance and variance. The figure reports a heatmap of the variance-normalized partial covariance between the two endogenous variables in Konig2017, which capture total fighting incidents involving allies and enemies ($\mathbf{G}Y$), and the rainfall-based instruments constructed using higher-order powers of the ally and enemy networks ($\mathbf{G}^k \mathbf{X}$).\footnote{For expositional purposes, the partial covariances and variances are computed without conditioning on additional exogenous controls that may enter the first-stage regressions.} Across most instruments, the resulting ratios are of the order $10^{-3}$ to $10^{-1}$, indicating very small partial first-stage coefficients. Together with the sparsity of the underlying ally and enemy networks (Figures (ref) and (ref)), these findings are consistent with a weak-identification environment in the sense of Definition (ref).
In which settings are network-based IVs likely to be ill-defined, lead to weak identification or standard identification and inference? To answer this question and to explore the above characterization, we examine the behavior of peers-of-peers instruments within the widely studied Erdős–Rényi (ER) random graph model. In an ER graph $G(n,p)$, each potential link between two nodes is formed independently with probability $p$. Thus, sparseness or density of the network is therefore governed entirely by $p$, or equivalently by the average degree $d(n) = n p(n)$. Although conceptually simple, these graphs have been widely studied and can form a foundation for more complex models by providing a useful benchmark (jackson10). Within Economics, for example, some structural models of network formation are asymptotically indistinguishable from ER graphs (mele17), they are used to model such phenomena as market entry and diffusion (e.g., campbell24) and they are also a basis for many simulation designs and comparisons (see graham20, for example).
Our characterization and proofs rely on results linking the spectral norm of the adjacency matrix to the degree distribution and use it to establish bounds on the first-stage coefficient. More specifically, we require $(\mathbf{I} - \beta \mathbf{G})$ to be invertible for the first stage covariance to be finite and bounded. To this end, we make the following assumption.
This assumption holds automatically when the degree sequence is uniformly bounded, a condition imposed in several empirical and theoretical network models to control the size of peer effects and maintain stable influence (e.g., dePaula2018, Leung2020). This is also a common assumption made by discrete choice peer effect models such as lambotte2025peer, required for existence and uniqueness of equilibrium with peer effects.
However, this assumption fails for many Erdős–Rényi graphs because the largest eigenvalue $\lambda_1^A(n)$ grows proportionally to $\max\{d_n,\sqrt{\Delta_n}\}$, where $d_n$ denotes the average degree and $\Delta_n$ the maximum degree (see krivelevich2001). Row normalization, as in bramoulle2009identification, ensures that the spectral norm of $\mathbf G_n$ is bounded by one and allows $(\mathbf I-\beta \mathbf G_n)^{-1}$ to be expressed as a convergent Neumann series.
Alternatively, one can normalize the adjacency matrix by an explicit scaling factor that depends on the growth rates of the average and maximum degree, thereby tracking the spectral norm while preserving symmetry of the network operator. Formally, let $\mathbf A_n$ be an Erdős–Rényi adjacency matrix with average degree $d_n = n p_n$ and maximum degree $\Delta_n$. We work with the scaled matrix \[ \mathbf G_n = \frac{1}{w_n}\mathbf A_n, \qquad w_n := \max\{d_n,\,\sqrt{\Delta_n}\}, \] which ensures that the largest eigenvalue of $\mathbf G_n$ is stochastically bounded and is close to unity. Assumption (ref) can then be modified as follows.
Assumption (ref) holds for any fixed $|\beta|<1$ and reduces to Assumption (ref) when $w_n$ is chosen to equal $\lambda_1^A(n)$. Consequently, our main results below consider both unscaled and scaled adjacency matrices, with the corresponding proofs provided in the Appendix.\footnote{Relatedly, Wang2025 study the distinction between row-normalized and scaled adjacency matrices in near-degree-regular networks, where agents’ local neighborhoods become asymptotically identical. They show that while row normalization can induce weak identification in such settings, appropriate scaling can mitigate these issues by aligning the rates of network regressors. Our analysis is more general, applies to networks with heterogeneous degree distributions, and explicitly leverages results from random graph theory.}
In the following proposition, we establish an upper bound on the variance-normalized covariance between the endogenous variable and the network-based instrument on ER graphs. We find that contingent on the regime that we are in (defined by the average degree), the unscaled and scaled models may lead to very different conclusions about weak instruments and identification in general. Note that for tractability, we work with the non-residualized variance--normalized covariance since it is a convenient proxy for the partial ratio in Definition (ref) and remains informative about first-stage relevance.
Proposition (ref) establishes that the variance--normalized covariance between the friends-of-friends instrument and the endogenous regressor is asymptotically bounded above by a function of the expected degree. Moreover, under Erd\H{o}s--R\'enyi asymptotics and the assumptions of the Proposition, the variance of the friends-of-friends instrument admits a sharp rate: \[ \frac{1}{n}\operatorname{Var}(\mathbf G_n^{(2)}\mathbf X) = \frac{\sigma_x^2}{n}\,\mathbb E\|\mathbf G_n^{(2)}\|_F^2 = \sigma_x^2\left(d_n^2 + \frac{d_n^4}{n}\right) + O\!\left(\frac{d_n^2}{n}+\frac{d_n^4}{n^2}\right). \] In contrast, the corresponding covariance admits an exact decomposition\footnote{This decomposition follows from substituting the linear representation $\mathbf Y=(\mathbf I-\beta\mathbf G_n)^{-1}(\alpha\iota+\gamma\mathbf X+\delta\mathbf G_n\mathbf X+\boldsymbol\varepsilon)$ into $\operatorname{Cov}(\mathbf G_n^{(2)}\mathbf X,\mathbf G_n\mathbf Y)$ and using the Neumann-series expansion $(\mathbf I-\beta\mathbf G_n)^{-1}=\sum_{k\ge0}\beta^k\mathbf G_n^k$, which converges under Assumption (ref).} in which the leading term is proportional to the same quantity: \[ \frac{1}{n}\operatorname{Cov}(\mathbf G_n^{(2)}\mathbf X,\mathbf G_n\mathbf Y) = \frac{\sigma_x^2}{n}\Big( (\beta\gamma+\delta)\,\mathbb E\|\mathbf G_n^{(2)}\|_F^2 + \mathbb E[R_n] \Big), \] where
This representation implies that the covariance cannot decay faster than the instrument variance: the leading term is of the same order as $\|\mathbf G_n^{(2)}\|_F^2$, while higher-order contributions are controlled by powers of $\beta$ under Assumption (ref). Consequently, when $d_n=o(1)$ and the graph collapses asymptotically, the variance of $\mathbf G_n^{(2)}\mathbf X$ collapses rapidly while the covariance does not vanish faster, so the variance--normalized covariance diverges. In this extremely sparse regime, the population first-stage estimand is therefore ill-defined.
At the other extreme, when the graph becomes asymptotically dense as $d_n$ increases, higher-order neighborhoods become nearly deterministic and $\mathbf G_n^{(2)}$ becomes asymptotically proportional to $\mathbf G_n$. Consequently, $\mathbf G_n^{(2)}\mathbf X$ becomes asymptotically collinear with $\mathbf G_n\mathbf X$, and the friends-of-friends instrument adds little independent variation beyond first-order neighbors. This mirrors the dense-network identification failure documented in bramoulle2009identification and Wang2025. In this regime, the upper bound in Proposition (ref) vanishes, the population first-stage coefficient $\pi_n$ shrinks, and identification fails due to asymptotic collinearity.
Between these extremes lies the empirically relevant case in which the expected degree is asymptotically bounded, $d_n=O(1)$. In this regime, the upper bound derived in Proposition (ref) does not force the population first-stage coefficient $\pi_n$ to diverge or vanish, so population relevance is not mechanically ruled out. Identification strength is governed instead by the information index $k(n)$ in accordance with Definition (ref). Weak identification in this regime arises when $k(n)$ remains bounded or, in fact, falls, so that sampling uncertainty in the first stage does not vanish with $n$. A natural case when $d_n=O(1)$ is with the presence of a large fraction of isolated or weakly connected nodes.\footnote{In Erdős–Rényi graphs with bounded average degree, the probability that a node has degree zero or one does not vanish asymptotically. In fact, when $d_n<1$, the graph fails to form a giant component (see erdds1959random) and hence a non-negligible fraction of nodes are isolated.} Thus, depending on $d_n$, increasing $n$ may not eliminate the presence of many isolated or near-isolated nodes with empty or small higher-order neighborhoods. As these nodes contribute little variation to $\mathbf G_n^{(2)}\mathbf X$, sampling noise may not vanish at standard rates and may accumulate extremely slowly, giving rise to weak instruments.
Our results show that scaling the adjacency matrix mitigates these issues by stabilizing the population first-stage coefficient. Under Assumption (ref), this scaling controls the spectral norm of the network operator and permits the average degree to grow while ensuring that the model remains well defined. In sparse regimes, where $d_n=o(1)$ or $d_n=O(1)$, the maximum degree $\Delta_n$ still grows even when average degree is bounded and hence $w_n=\sqrt{\Delta_n}$ asymptotically. The variance and covariance of the scaled friends-of-friends instrument then satisfy \[ \operatorname{Var}(\mathbf G_n^{(2)}\mathbf X) = O\!\left(\frac{1}{w_n^4}\,\mathbb E\|\mathbf A_n^{(2)}\|_F^2\right), \qquad \operatorname{Cov}(\mathbf G_n^{(2)}\mathbf X,\mathbf G_n\mathbf Y) = O\!\left(\frac{1}{w_n^3}\,\mathbb E\|\mathbf A_n^{(2)}\|_F^2\right), \] so normalization aligns the rates at which the covariance and variance decay and prevents explosive behavior of the first-stage ratio.
Because scaling stabilizes the variance of the network-based instrument, weak identification in the scaled specification arises only if the effective information index $k(n)$ grows too slowly relative to $w_n$. In particular, when the upper bound on $|\pi_n|$ grows at most linearly in $w_n$, strong identification requires that $k(n)$ diverges faster than $w_n$ so that sampling uncertainty in the first stage vanishes. Conversely, if $k(n)=O(w_n)$ or grows more slowly, then the scaled local parameter remains bounded and weak identification persists despite normalization. Thus, scaling shifts the source of weakness from pathological variance collapse to a transparent comparison between information accumulation and network scaling. Even though the amount of independent information provided by the instrument is unchanged, scaling aligns the growth rates of these moments across regimes.\footnote{As shown in Wang2025, in near-degree-regular networks this adjustment restores standard $\sqrt{n}$ convergence rates.}
Assumptions (ref)--(ref) ensure that $(\mathbf I-\beta\mathbf G_n)$ is invertible and that the reduced-form is well defined. Lemma (ref) in Appendix Section (ref) shows, however, that invertibility alone does not guarantee stable or well-behaved first-stage relationships. In particular, conditional on $\mathbf G_n$, the population first-stage covariance admits the decomposition \[ \frac{1}{n}Cov(\mathbf G_n^{(2)}\mathbf X,\mathbf G_n\mathbf Y\mid \mathbf G_n) = \frac{\sigma_x^2}{n}\sum_{j=1}^n \frac{\lambda_j^3(n)\big(\gamma+\delta\lambda_j(n)\big)}{1-\beta\lambda_j(n)} \;+\; \frac{\sigma_x^2}{n}\,R_{n,\mathrm{diag}}, \] where the leading term aggregates the contributions of the spectral components of the network operator and the remainder term arises from the diagonal adjustment in $\mathbf G_n^{(2)}$. Moreover, the remainder satisfies the deterministic bound \[ |R_{n,\mathrm{diag}}| \;\le\; \|\mathbf D_n\|_F\, \|\mathbf G_n\|_2\, \|(\mathbf I_n-\beta\mathbf G_n)^{-1}\|_2\, \|\gamma\mathbf I_n+\delta\mathbf G_n\|_2, \] and therefore does not introduce additional amplification through factors of $(1-\beta\lambda_j(n))^{-1}$.\footnote{Note that the magnitude of $R_{n,\mathrm{diag}}$ depends on $\|\mathbf D_n\|_F^2$ as long as $(I_n - \beta \mathbf{G}_n)^{-1}$ is invertible and bounded. In the scaled version, taking expectations and using $\deg(i)\sim\mathrm{Bin}(n-1,p_n)$ gives \[ \mathbb E\|\mathbf D_n\|_F^2 =\frac{n}{w_n^4}\,\mathbb E[\deg(i)^2] =\frac{n}{w_n^4}\big(d_n(1-p_n)+d_n^2\big). \] Hence, by Markov's inequality, \[ \frac{\|\mathbf D_n\|_F}{n} =O_p\!\Big(\frac{\sqrt{d_n^2+d_n}}{w_n^2\sqrt n}\Big). \] If, in addition, $\|(\mathbf I_n-\beta\mathbf G_n)^{-1}\|_2=O_p(1)$ and $\|\mathbf G_n\|_2=O_p(1)$, then implies \[ \frac{1}{n}|R_{n,\mathrm{diag}}| \le \frac{\|\mathbf D_n\|_F}{n}\, \|\mathbf G_n\|_2\, \|(\mathbf I_n-\beta\mathbf G_n)^{-1}\|_2\, \|\gamma\mathbf I_n+\delta\mathbf G_n\|_2 = o_p(1), \] so the diagonal correction is asymptotically negligible and cannot affect the sign of the first-stage covariance provided the leading spectral term is bounded away from zero. Here $w_n=\max\{d_n,\sqrt{\Delta_n}\}$ with $\Delta_n$ the maximum degree. A similar conclusion holds in the unscaled case when the expected degree is uniformly bounded, $d_n=O(1)$.} Two important phenomena follow directly from this spectral decomposition.
To illustrate the previous discussions, in Figures (ref) and (ref), we report simulated upper bounds for the unscaled and scaled adjacency matrices described in Proposition (ref). For each $n\in\{200,400,800,1200,1600\}$ we generate $500$ Monte Carlo draws under several Erd\H{o}s--R\'enyi regimes.
From Figure (ref), we observe that in the unscaled model, the bound is approximately flat as n increases in the constant-degree and $\log\log n$ regimes. By contrast, in the extremely sparse regime ($d_n=o(1)$), the bound increases sharply as $n$ grows. This behavior is consistent with the analytical discussion above about the instability of the population first-stage estimand. In the dense regime where $d_n\to\infty$, the bound decreases toward zero, consistent with the fact that $\mathbf G_n^{(2)}$ becomes nearly collinear with $\mathbf G_n$, so the instrument adds little independent variation.
In the scaled specification, where $\mathbf G_n=\mathbf A_n/w_n$ and shown in Figure (ref), the simulated bound does not exhibit the explosive behavior seen in the extremely sparse regime. Rather, the average upper bound remains well behaved across all the regimes: it does not diverge, and it does not fall quickly. This is consistent with the role of scaling in Proposition (ref). Consequently, in the scaled model, weak identification is governed by the information index $k(n)$ in Definition (ref).
To connect these theoretical results to the linear-in-means model, we use the same four Erd\H{o}s--R\'enyi cases considered in the previous figures, but incorporate these networks within the linear-in-means model studied above. Following the set-up in bramoulle2009identification, we fix the network structure for each $n$ and, conditional on the graph, draw the error term $\varepsilon$ and uni-variate $\mathbf{X}$ independently across Monte Carlo repetitions. For each $n\in\{200,400,800,1200,1600\}$ we use $500$ Monte Carlo draws and compute the average first-stage $F$-statistic, the sample covariance between the endogenous regressor $\mathbf G\mathbf Y$ and the instrument $\mathbf G^{(2)}\mathbf X$, and the variance of the instrument. These are reported in Figures (ref), (ref), and (ref).
As we can see in Figure (ref), for the unscaled adjacency matrix, the first-stage $F$-statistic remains flat and close to zero across all regimes except the extremely sparse case. In the extremely sparse regime, the $F$-statistic increases with sample size, reflecting the instability of the first-stage ratio. This is because the instrument variance goes to 0, rather than genuine accumulation of identifying information. In contrast, for the scaled adjacency matrix, the first-stage $F$-statistic increases monotonically with $n$ across regimes, suggesting a strong and stable first stage.
Figures (ref) and (ref) further support this analysis. In the unscaled case, the covariance between the endogenous regressor and the instrument remains close to zero in most regimes, while the variance of the instrument is large except in the extremely sparse regime, where it collapses as higher-order neighborhoods vanish. By contrast, in the scaled specification, both the covariance and the variance of the instrument remain stable across sample sizes once the extremely sparse regime is excluded. Although normalization reduces the magnitude of these quantities in levels, it aligns their rates so that the variance--normalized covariance stabilizes as $n$ grows. This behavior is consistent with the theoretical upper bounds reported in Figure (ref) and illustrates how scaling by $w_n$ mitigates spurious first-stage behavior driven by sparsity or degree heterogeneity.
While the upper-bound plots and first-stage diagnostics (along with the theoretical results) presented above suggest that scaling the adjacency matrix mitigates many of the pathologies present in the unscaled model, scaling is not a panacea that can solve all possible perils with network based instruments. In particular, scaling stabilizes the population first stage across regimes, but it does not, by itself, rule out weak identification in some cases. As mentioned above, near-boundary instability can still arise in the scaled model and can lead to the problem of weak instruments. We explore this issue further using simulations in Section (ref).
Given that multiple regimes in Proposition (ref) have identification failures at the limit, the instruments based on peers-of-peers is likely to be weak empirically for those configurations, especially for the unscaled model. As a solution, we follow the literature on weak instruments such as moreira03 in considering inference that is robust to a weak first stage.
In the classical Instrumental Variables model, we could proceed by either implementing the Anderson-Rubin (AR) test from anderson1949estimation (known to be unbiased and asymptotically efficient in the just-identified case - see moreira2009tests) or the Conditional Likelihood Ratio test from moreira03. The AR-statistic is given by:
where $\sqrt{n}\hat{g} (\beta_0) = \sqrt{n}(\hat{\xi} - \beta_0 \hat{\pi}) \stackrel{H_0}{\to}_d N(0,\Omega(\beta_0))$ with an appropriate variance estimator for
Equation (ref) holds regardless of the strength of the instrument and, thus, the AR test is given as $\phi^{AR}_{n}(\alpha) = \mathbf{1}\{AR(\beta_0) > \chi^2_{k,1-\alpha} \}$ for the null hypothesis $H_0:\beta = \beta_0$. The confidence set of the test can take different forms including the extreme case where it is the entire real line when $\pi=0$. This is because $\beta$ is not identified in that case and any value of beta satisfies the restriction condition (the test will have zero power in this case).\footnote{In the over-identified case the AR test is still robust and unbiased, but it is no longer efficient under strong identification. The conditional likelihood test (CLR) of moreira03 can then be used in the homoskedastic case. We focus on the just-identified case and, thus, specialize to the AR test.}
However, there is one main distinction of our system (ref)-(ref) relative to the classical set-up: the heteroskedasticity induced by network-dependency of $\varepsilon$. This is the term $(\mathbf{I} - \beta \mathbf{G})^{-1} \mathbf{\varepsilon}$ in (ref), implying errors of the form,
causing them to be correlated across the connections even if $\varepsilon$ is assumed to be homoskedastic.
We provide two possible approaches to inference robust to weak instruments. First, to fully deal with the cross-sectional dependence arising from network spillovers, we propose the use of the variance estimator in kojevnikov2021limit (Proposition 4.3). As that variance estimator is consistent for $\Omega(\beta)$ given $\beta$, an application of Slutsky's Lemma guarantees the applicability of the feasible AR-test. This is summarized in the proposition below.\footnote{Other papers, such as acemoglu2015, model the cross-sectional dependence as spatially correlated data and use conley1999gmm as a consistent estimator for $\Omega(\beta)$. Alternatively, it is common to use clustered variance estimators. The latter requires a block structure (e.g., independence across villages, schools, families) and an asymptotic theory based on "many" networks.}
However, we note that even the homoscedastic implementation of the Anderson-Rubin test and, thus, the Conditional Likelihood Ratio test (CLR) perform well asymptotically. This is because the higher-order terms in (ref) are likely to be negligible when the network is sparse and converging to 0 and $\beta$ is small. Indeed, we have that
as $ \beta^k \mathbf{G}^k \to \mathbf{0}$.
Thus, asymptotic inference based on homoskedastic errors (e.g., moreira03) is likely to perform very well when network-based instruments are weak due to network sparsity, as the heteroskedasticity due to networks disappears when the network becomes very sparse. Similarly, when the network becomes close to complete, the variance converges to:
as $\mathbf{G}^k \to \mathbf{G}$, which is still homoskedastic.
We now provide Monte Carlo simulations to illustrate the finite-sample properties of our theoretical results. We base our data-generating process on those used in bramoulle2009identification, but with alternate network structures that showcase the issue of weak identification with peers-of-peers instruments for both the scaled and unscaled adjacency matrices.
We consider Erdös-Rényi random graphs (erdds1959random) with $\textit{d = np}$ being the average degree of the graph. We vary d across settings, taking values from the set $\{0.25,0.5,0.75,1,2,5\}$ with sample size $n$ varying in $\{250, 500, 1000, 2000\}$. Notice that the smaller the value of d, the sparser the network. For every setting, we fix d and n and set the Monte Carlo Simulation number to 1000, keeping the network structure fixed across all runs. Following bramoulle2009identification, we draw a uni-dimensional $X_i$ with approximately $5\%$ of values to be 0\footnote{Using a Bernoulli(0.9458333) as in bramoulle2009identification.} and the remaining $95\%$ follow an i.i.d. log-normal ($X_i \stackrel{i.i.d.}{\sim} LogNormal(1,3)$). These are fixed for a given data size and average degree. We draw the error terms $\varepsilon_i \stackrel{i.i.d.}{\sim} Normal(0,1)$, independently for each run of the simulation. The true coefficients are set at $\alpha = 0.7683, \gamma = 0.0834, \delta = 0.1507$, and $\beta \in \{0.466, 0.95\}$ to evaluate the impact of a change in intensity of the peer effect.
We estimate all parameters using Two-Stage Least Squares (TSLS), with the instrument $\mathbf{G}^{(2)} X$ for the endogenous variable $\mathbf{G}Y$, and the exogenous variables as instruments for themselves. We run the simulations for two models - the unscaled model where $\mathbf{G}_n = \mathbf{A}_n$ and the scaled model where $\mathbf{G}_n = \frac{1}{w_n}\mathbf{A}_n$ and $w_n = \max\{d_n,\sqrt{\Delta_n}\}$ with $\Delta_n$ being the maximum degree of the network. The results are shown in Tables (ref) - (ref) below. Appendix Section (ref) contains extended results from the simulation exercise.
Table (ref) reports the average TSLS estimates of $\beta$, together with the correlation between the endogenous regressor $\mathbf G\mathbf Y$ and the peers-of-peers instrument $\mathbf G^{(2)}\mathbf X$, and the corresponding first-stage $F$-statistic. Several features are worth highlighting. First, in the unscaled model, when the average degree is less than one, we observe very high correlations and large first-stage $F$-statistics. This behavior coincides with our theoretical results regarding extremely sparse regimes: in this extremely sparse regime the variance of the instrument collapses, so that the first-stage ratio becomes unstable and the population first stage is ill-defined asymptotically. As shown in Table (ref) in the Appendix, the covariance between $\mathbf G^{(2)}\mathbf X$ and $\mathbf G\mathbf Y$ is in fact extremely small in this regime, indicating that the large correlations are mechanically driven by vanishing instrument variance rather than genuine identifying power. Furthermore, as the average degree increases beyond one, the unscaled model exhibits very low correlations and weak first-stage $F$-statistics across all sample sizes. This pattern persists as $n$ grows, indicating that the resulting weak instruments problem is structural rather than a small-sample artifact.
Turning to the scaled model, we find that for $\beta=0.4666$ the instrument is strong: both the correlation and the first-stage $F$-statistic increase monotonically with $n$ for all values of average degree. However, by contrast, when $\beta=0.95$ the correlation between $\mathbf G\mathbf Y$ and $\mathbf G^{(2)}\mathbf X$ exhibits sign reversals and the first-stage $F$-statistic deteriorates before increasing again. This behavior reflects the near-boundary phenomenon discussed in Section (ref). Although scaling stabilizes the spectrum of the network operator and prevents explosive behavior, finite-sample realizations may still feature eigenvalues $\lambda_j$ such that $\beta\lambda_j$ is close to one. As the average degree increases further away from one, the correlation moves away from zero, consistent with the attenuation of this near-boundary effect. Hence, the simulations illustrate situations where even the scaled model might suffer from a weak first-stage.
To complement these findings, Table (ref) reports empirical coverage probabilities at the 95% nominal level for three inference procedures: the conventional $t$-test with homoscedastic standard errors, the $t$-test using Kojevnikov (kojevnikov2021limit) standard errors, and the Anderson--Rubin (AR) test with homoscedastic errors.\footnote{We also implement the AR test with Kojevnikov standard errors (reported in Table (ref) in Appendix). Given the homoscedastic data generating process, the resulting coverage rates are very similar—especially for larger samples.} Consistent with the discussion in Section (ref) and the first-stage diagnostics in Table (ref), we find that for the unscaled model the $t$-test fails to deliver correct coverage once the average degree is greater than or equal to one. In these regimes, coverage probabilities are often close to 100%, reflecting extremely wide confidence intervals—an outcome symptomatic of weak instruments and identification failure.
Even in the extremely sparse regime, where coverage improves toward the nominal 95% level as $n$ increases, we continue to observe systematic over-coverage. Using Kojevnikov standard errors mitigates some of these distortions, but the resulting confidence intervals remain conservative in many cases. This indicates that the problem cannot be resolved by variance correction alone and instead reflects a deeper identification issue.
To better understand the source of these distortions, we focus on two representative cases in the unscaled model: one in which the first-stage appears relatively strong ($d=0.5$) and coverage is close to nominal, and one in which the instrument is weak ($d=2$) and coverage is severely distorted. Figures (ref) and (ref) display the empirical distributions of the TSLS estimator $\hat{\beta}$ and the associated $t$-statistic computed using Kojevnikov standard errors, for $n=2000$ and $\beta_0=0.4666$ across 1,000 Monte Carlo replications. In both cases, the sampling distribution of $\hat{\beta}$ remains centered near the true parameter, but its shape varies sharply with instrument strength. Under a presumably strong first-stage, the distribution is approximately normal, in line with standard large-sample TSLS theory. Under weak instruments, however, $\hat{\beta}$ exhibits heavy tails and pronounced non-Gaussian features; in weakly identified settings, IV estimators are known to have Cauchy-like or multimodal finite-sample distributions, invalidating normal approximations and conventional inference (staigerstock97, andrews2019weak). This behavior is mirrored in the distribution of the corresponding $t$-statistics, which display clear departures from normality and pronounced multi-modality, explaining the observed over-coverage of confidence intervals.
A practical mechanism underlying these pathologies is instability in the first-stage projection. When the first stage is nearly singular, the matrix inversions that enter the TSLS variance formula become ill-conditioned, producing highly variable—and occasionally extremely large—standard errors. This suggests that TSLS standard errors may be poorly behaved in finite samples under weak identification, leading to non-standard distributions for the associated $t$-statistics. Figure (ref) illustrates this mechanism by plotting the first-stage $F$-statistic against the Kojevnikov-based TSLS standard errors. In the strong-first stage case, the relationship is approximately monotone and decreasing, with relatively few outliers. In contrast, under weak instruments the relationship becomes highly nonlinear: a nontrivial mass of realizations produces inflated standard errors even when the realized $F$-statistic is not particularly small. This further underscores that the observed distortions cannot be remedied by standard error corrections alone.
In contrast, AR test, even with homoskedastic standard errors, delivers coverage probabilities close to the nominal level across all degree regimes and sample sizes. This holds both in sparse and dense networks, and persists even in cases where the first stage is weak and TSLS-based inference severely over-covers. The robustness of the AR test reflects its invariance to weak identification and remains valid even when the reduced-form and first-stage relationships are unstable. Meanwhile, the scaled specification shows substantially more stable coverage behavior across degree regimes, although weak-instrument is likely when $\beta \to 1$.
In this section we look at two applications of network-based instruments which are related to our setting: an application in development (conflict across groups, in Konig2017) and one in trade (using gravity equations, from brancaccio2020).
Konig2017 investigate strategic complementarities in the use of violence within a network of armed groups during the Second Congo War and examine how networks of alliances and hostilities influence the intensity of conflict. Nodes correspond to armed actors, while links encode military relationships, distinguishing between alliances (groups fighting on the same side), enmities (groups that directly clash), and neutrality (groups that are neither allies nor enemies). They model conflict as a network game, and derive the corresponding Nash equilibrium with the optimal level of fighting intensity depending on the fighting of its allies and enemies. Endogeneity of these network-sum regressors is addressed using network-based instrumental variables constructed from exogenous weather shocks. In particular, rainfall in the homeland of linked groups is used as an excluded instrument, and—following bramoulle2009identification—the authors explicitly exploit second-degree instruments based on the rainfall of neighbors-of-neighbors.
We focus on the main empirical specification reported in Table 1 of Konig2017, which examines how a group’s own fighting effort responds to the fighting efforts of its network neighbors. The empirical model is over-identified and features three endogenous regressors: total fighting effort of allies (TFA), total fighting effort of enemies (TFE), and total fighting effort of neutral groups (TFN). Rainfall shocks defined on higher-order network neighborhoods, such as rainfall in the territories of allies’ allies and enemies’ enemies instruments correspond exactly to second-degree network instruments of the form $\mathbf G^{(2)}\mathbf X$ in our framework. The three most relevant specifications from Table 1 are reproduced in Tables (ref), (ref), and (ref), where we report point estimates, 95% non-robust confidence intervals based on clustered $t$-tests, and projected 95% Anderson–Rubin and conditional likelihood ratio confidence intervals.\footnote{The authors use a custom spatial TSLS estimator in Stata to account for spatial correlation. We instead use clustered standard errors via the ivreg2 command, which is also used in their replication files and is compatible with the weakiv package of weakiv2013. As a result, reported standard errors may differ slightly from those in Konig2017, while point estimates remain identical.}
The non-robust confidence intervals reported in Table 1 of Konig2017 would suggest statistically significant effects of allies’, enemies’, and neutral groups’ fighting efforts, thereby indicating strong first-stage relevance. However, once weak-IV–robust inference is applied, this conclusion is overturned. Across the same specifications, the Anderson–Rubin and conditional likelihood ratio confidence intervals not only include zero but—with the exception of a single case—are unbounded for all three endogenous variables. This pattern is characteristic of a weak-instrument environment in which the covariance between the instruments and the endogenous regressors is small relative to sampling variability.
The source of this weakness is the sparse structure of the underlying alliance–enmity networks. As shown in Table (ref), the networks have an average degree of approximately $2.4$--$3$ and a median degree of 1, implying that second-degree connections are limited. Consequently, the friends-of-friends instruments constructed from rainfall shocks provide little independent variation and are only weakly correlated with the corresponding network-sum regressors. This weak first-stage relationship is also evident in Figure (ref), which displays near-zero partial variance-normalized covariances between the endogenous variables and the network-based instruments. Together, these findings illustrate how standard inference can be misleading in network settings: non-robust confidence intervals mask weak identification due to sprasity and limited higher-order neighborhoods, while weak-IV–robust methods correct for weak identification. This empirical application therefore underscores the importance of weak-IV–robust inference when working with network-based instruments, in line with our theoretical and simulation results, and extends the simulation-based insights of andrews2019weak to a realistic empirical environment.
A second application of weak-IV–robust inference with network-based instruments is provided by brancaccio2020, who study how the global transportation network interacts with world trade and help shape trade flows and trade costs alongside traditional geographic determinants. Their key departure from the standard gravity-model framework is to treat transportation costs as endogenous, rather than relying on the iceberg-type trade costs of samuelson1954.\footnote{Transport costs are assumed to be linear in distance and are extracted in units of arriving volume, in contrast to the iceberg formulation where goods “melt away” proportionally during transit.} To this end, the authors develop a structural spatial model of the maritime transport network, using detailed data on ship movements—particularly in the dry bulk shipping sector—shipping prices, and port-to-port connections.
A central feature of their setting is that shipping prices exhibit strong network effects. Because ships frequently undertake ballast (empty) voyages, carriers choose routes by jointly considering spot prices at the destination port and at nearby ports in the broader region. As a result, shipping costs are shaped not only by bilateral conditions but also by congestion and demand elsewhere in the network. Within this framework, brancaccio2020 construct instrumental variables for shipping costs that exploit variation induced by the structure of the maritime network. These instruments therefore inherit the sparsity and interdependence of the underlying transportation network, raising concerns about weak identification—precisely the type of environment in which weak-IV–robust inference is essential.
The growth in measuring peer-effects in academia and policy should also bring renewed attention to the challenges in inference. We showed how undefined first-stage estimands and/or weak instruments can arise naturally in the linear-in-means model with network-based instruments. As discussed, the issue stems from the sparsity (or extreme density) of existing networks, and how these characteristics persist with higher order counterparts of linkages across peers. We characterized the conditions in terms of often used primitive networks, like Erdos-Renyi, and how they can affect inference. Yet, we provide solutions based on developments of weak-IV robust testing and consistent variance estimation for linear regression models with networked random variables.
Such concerns are relevant and, as we discussed, examples of sparse networks span many fields and salient examples in Economics (see Table (ref)). While our characterization results focus on the linear-in-means model, our main intuition extends to any linear regression framework with network-based instruments. Indeed, definition (ref) for weakness of network-based IVs do not rely on that specific model. Furthermore, issues with first-stage estimands and lack of information for higher-order networks would be prevalent in all such cases. Thus, we deem that such concerns are warranted in many more applications. Finally, it is likely that such issues extend to non-linear models of network interactions, such as discrete choice models with network-based instruments (e.g., volpe2025discrete), among others. Future work should investigate identification and inference in such settings.