Identification and Estimation of Network Models with Nonparametric Unobserved Heterogeneity

Identification and Estimation of Network Models with Nonparametric Unobserved Heterogeneity

Introduction

Unobserved heterogeneity is pervasive in economics. The importance of accounting for unobserved heterogeneity is well recognized in microeconometrics, in general, as well as in the network context, in particular. For example, since \citet*{Abowd1999}, estimating a linear regression with additive (two-way) fixed effects has become a standard approach to analyzing interaction data. Originally employed to account for workers and firms fixed effects in the wage regression context, this technique has become a standard tool to control for two-sided unobserved heterogeneity and decompose it into agent specific effects.\footnote{For example, recent applications to employer-employee matched data feature \citet*{Card2013,Helpman2017,Song2019} among others. The numerous applications of this approach also include the analysis of students-teachers \citep*{Hanushek2003,Rivkin2005,Rothstein2010}, patients-hospitals \citep*{Finkelstein2016}, firms-banks \citep*{Amiti2018}, and residents-counties matched data \citep*{chetty2018impacts}.} Since the seminal work of \citet*{Anderson2003}, the importance of controlling for exporters and importers fixed effects has also been well acknowledged in the context of the international trade network, including nonlinear settings of \citet*{SantosSilva2006} and \citet*{Helpman2008}. \citet*{Graham2017} stresses the importance of accounting for agents' degree heterogeneity (captured by the additive fixed effects) in network formation models.

While the additive fixed effects framework is commonly employed to control for unobservables in networks, it is not flexible enough to capture more complicated forms of unobserved heterogeneity, which are likely to appear in many settings. This concern can be vividly illustrated in the context of estimation of homophily effects, one of the main focuses of the empirical network analysis.\footnote{The term homophily typically refers to the tendency of individuals to assortatively match based on their characteristics. For example, individuals tend to form social connections based on gender, race, age, education level, and other socioeconomic characteristics. Similarly, countries that share a border, have the same legal system, language or currency, are more likely to have higher trade volumes.} Since homophily (assortative matching) based on observables is widespread in networks (e.g., \citealp*{McPherson2001}), homophily based on unobservables (fixed effects) is also likely to be an important determinant of the interaction outcomes. Since observed and unobserved characteristics (i.e., covariates and fixed effects) are typically correlated, the presence of latent homophily significantly complicates identification of the homophily effects associated with the observables (e.g., \citealp*{Shalizi2011}). Failing to properly account for homophily based on unobservables (and other complex forms of unobserved heterogeneity, in general) is likely to result in inconsistent estimators and misleading policy implications.

To address the concern discussed above, we consider a dyadic network model with a flexible (nonparametric) form of unobserved heterogeneity, where the outcome of the interaction between agents $i$ and $j$ is given by

Here, $W_{ij}$ is a $p \times 1$ vector of pair-specific observed covariates, $\beta_0 \in \mathbb R^p$ is the parameter of interest, $\xi_i$ and $\xi_j$ are unobserved fixed effects, and $\varepsilon_{ij}$ is an idiosyncratic error. The fixed effects are allowed to interact via the coupling function $g(\cdot,\cdot)$, which is treated as $\emph{unknown}$. Importantly, we do not require $g(\cdot,\cdot)$ to have any particular structure and do not specify the dimension of $\xi$. Finally, $F(\cdot)$ is a known (up to location and scale normalizations) invertible link function. The presence of $F(\cdot)$ ensures that (ref) is flexible enough to cover a broad range of the previously studied dyadic network models with unobserved heterogeneity, including important nonlinear specifications such as network formation or Poisson regression models.\footnote{For example, with $F(\cdot)$ equal to the logistic CDF and $g(\xi_i,\xi_j) = \xi_i + \xi_j$, (ref) corresponds to the network formation model of \citet*{Graham2017}.}

Being agnostic about the dimensions of the fixed effects and the nature of their interactions, (ref) allows for a wide range of forms of unobserved heterogeneity, including homophily based on unobservables.

We study identification and estimation of (ref) under the assumption that we observe a single network of a growing size.\footnote{The large single network asymptotics is standard for the literature focusing on identification and estimation of network models with unobserved heterogeneity. See, for example, \citet*{Graham2017,Dzemski2018,Candelaria2016,Jochmans2018,Toth2017,gao2020nonparametric,gao2023logical}.} First, we focus on a simpler version of (ref)

We argue that the outcomes of the interactions can be used to identify agents with the same values of the unobserved fixed effects. Specifically, we introduce a certain pseudo-distance $d_{ij}$ measuring similarity between agents $i$ and $j$ in terms of their latent characteristics. We argue that (i) $d_{ij} = 0$ if and only if $\xi_i = \xi_j$, and (ii) $d_{ij}$ is identified and can be estimated from the data. Consequently, agents with the same values of $\xi$ can be identified based on the pseudo-distance $d_{ij}$. Then, the variation in the observed characteristics of such agents allows us to identify the parameter of interest $\beta_0$ while controlling for the impact of the fixed effects. Importantly, this result is not driven by the linearity or particular functional form of (ref): we show that a similar identification argument also applies when $W_{ij}'\beta_0$ is replaced by an unknown (nonparametric) function of observables.

Demonstrating that the introduced pseudo-distance $d_{ij}$ is identified plays a central role in the argument outlined above. First, we show that, if the idiosyncratic errors $\varepsilon_{ij}$ are homoskedastic, $d_{ij}$ can be readily estimated from the data using simple pairwise-difference regressions. Second, we extend this argument to models with general heteroskedasticity by leveraging and advancing recent developments in the matrix estimation/completion literature. Specifically, we demonstrate that the error free outcomes $Y_{ij}^* \coloneqq W_{ij}' \beta_0 + g(\xi_i,\xi_j)$ are identified and can be uniformly (across all pairs of agents) consistently estimated. This is a powerful result allowing us to treat $Y_{ij}^*$ as effectively observed and thus greatly simplifying the analysis. In particular, working with $Y_{ij}^*$ instead of $Y_{ij}$ effectively reduces (ref) to a model without the error term $\varepsilon_{ij}$, which can be interpreted as an extreme form of homoskedasticity. This, in turn, allows us to establish identification of $\beta_0$ by applying the same argument as in the homoskedastic model.

Building on these ideas, we construct an estimator of $\beta_0$ and characterize its rate of convergence. Following the identification argument, we first estimate the pseudo-distances $\hat d_{ij}$, which are used to find agents with similar latent characteristics. Then, we estimate $\beta_0$ by combining pairwise-difference regressions of $Y_{ik} - Y_{jk}$ on $W_{ik} - W_{jk}$ for all pairs of agents $i$ and $j$ sufficiently similar in terms of $\hat d_{ij}$. Consistency of this estimator crucially relies on the matched agents being different in terms of their observables, resulting in sufficient residual variation in $W_{ik} - W_{jk}$. We characterize the asymptotic behavior of the estimator's bias due to imperfect matching, and provide conditions under which it is consistent.

In the general heteroskedastic case, construction of $\hat d_{ij}$ also involves preliminary estimation of the error free outcomes $Y_{ij}^*$. To provide estimators of $Y_{ij}^*$ and $d_{ij}$ valid under heteroskedastic errors, we build on and extend the approach of \citet*{Zhang2017} originally employed in the context of nonparametric graphon estimation. Moreover, to formally establish identification of the error free outcomes, we also propose a modified version of Zhang2017's estimator and demonstrate its consistency in the max (matrix) norm, meaning that, in a large network, our estimator recovers $Y_{ij}^*$'s for all pairs of agents with a high precision at once. To the best of our knowledge, this result is also new to the statistics literature on graphon estimation, which has previously focused on constructing estimators $\hat Y_{ij}^*$ and deriving their rates of convergence in terms of the mean square error/Frobenius norm Chatterjee2015,Gao2015,Klopp2017,Zhang2017,li2019nearest.

Finally, we want to stress that identification of the error free outcomes is a powerful result, which applies to a general class of dyadic network models beyond (ref) and can be used as a foundation for establishing new identification results. To the best of our knowledge, this result has not been previously recognized and leveraged in the econometrics literature. In particular, building on it, we demonstrate how the proposed identification and estimation strategies can be naturally extended to cover model (ref), as well as its nonparametric version. We also argue that the pair-specific fixed effects $g_{ij} = g(\xi_i,\xi_j)$ are identified for all pairs of agents $i$ and $j$ and can be (uniformly) consistently estimated. Identification of $g_{ij}$ is an important result in itself since in many applications the fixed effects are the central objects of interest. Moreover, this result is also of special significance when $F$ is nonlinear because identification of $g_{ij}$ allows us to identify important policy relevant quantities such as pair-specific and average partial effects.

This paper contributes to the literature on econometrics of networks and, more generally, two-way models. The distinctive feature of our model is allowing for flexible nonparametric unobserved heterogeneity: the fixed effects can interact via the unknown coupling function $g(\cdot,\cdot)$. Importantly, we do not require $g(\cdot,\cdot)$ to have any particular structure (other than satisfying a weak smoothness requirement) and do not specify the dimensionality of the fixed effects. This is in contrast to most of the existing approaches, which either explicitly specify the form of $g(\cdot,\cdot)$ or impose additional restrictive assumptions on its shape and smoothness.

Among explicitly specified forms of $g(\cdot,\cdot)$, the additive fixed effects structure ${g(\xi_i,\xi_j) = \xi_i + \xi_j}$ is by far the most popular way of incorporating unobserved heterogeneity in dyadic network models, e.g., see Graham2017, Charbonneau2017, Jochmans2018, Dzemski2018, \citet*{Yan2019}, Candelaria2016, gao2020nonparametric and Toth2017 among others. While this specification provides a practical way of controlling for degree heterogeneity, it does not account for more complicated forms of unobserved heterogeneity including latent homophily. Recent semiparametric extensions of Graham2017 and related frameworks feature gao2023logical relaxing separability between $W_{ij}' \beta_0$, $\xi_i$, and $\xi_j$. However, gao2023logical still require $\xi$ to be scalar and assume that the linking probability is increasing in $\xi_i$ and $\xi_j$, thus ruling out homophily based on unobservables.

The linear factor specification $g(\xi_i,\xi_j) = \xi_i' \xi_j$ is widely used in both network and panel models as a generalization of the additive fixed effect framework.\footnote{Recent studies considering network models with unobserved effects having a linear factor structure include, among others, chen2021nonlinear,ma2022detecting,zeleneev2025tractable.} While, as argued in \citet*{chen2021nonlinear}, this specification allows for certain forms of latent homophily, approximating a general function $g(\cdot, \cdot)$ by the linear factor model requires a growing number of factors (e.g., fernandez2021low). Such low-rank approximations, for example, are utilized by freeman2023linear and beyhum2024inference who consider a panel variation of (ref). To control the accuracy of the proposed low-rank approximations, both of these papers require $g(\cdot,\cdot)$ to have sufficiently many continuous derivatives and, more importantly, to ensure consistency of their estimator of $\beta_0$, they require regressors $W$ to be “high-rank”. Both of these requirements are restrictive in the network setting;\footnote{Non-differentiable functions $g(\cdot,\cdot)$, such as the one provided in Example (ref) for $\zeta = 1$, are a common feature of popular latent homophily models (e.g., Hoff2002,handcock2007model).} see Section (ref) for a detailed comparison of the frameworks and approaches to identification. Allowing for both low-rank regressors $W$ and possibly non-differentiable $g(\cdot,\cdot)$ is an important and unique combination of features of our setting differentiating this paper from the rest of the literature.\footnote{Low-rank regressors are typically ruled out to ensure identification of $\beta_0$ even in linear factor models; see, for example, bai2009panel,moon2015linear,armstrong2022robust.}

Employing clustering methods to partition units into groups with similar unobserved characteristics is another method commonly employed to control for interactive (time-varying) unobserved heterogeneity. This approach was proposed and applied to likelihood panel models in the seminal work by bonhomme2022discretizing and then extended to semiparametric panel regressions by beyhum2024inference. While the specific clustering methods (involving the observed covariates as well) employed in these papers proved to be instrumental in panel settings, they would not allow one to identify and estimate $\beta_0$ in network models when the regressors take the typical form of $W_{ij} = w(X_i,X_j)$, where $X_i$ and $X_j$ are observed characteristics of agents $i$ and $j$. In this case, the previously employed methods would cluster units into groups with similar values of both $\xi$ and $X$ preventing one from disentangling the effects of observables and unobservables and thus from identifying $\beta_0$ as well.\footnote{Moreover, the clustering methods employed in bonhomme2022discretizing and beyhum2024inference are based on individual-specific moments. Such clustering methods might fail to meaningfully group units in network settings, in which informative individual-specific moments might not exist at all; see Section (ref) and specifically Footnote (ref) for details.} Importantly, unlike the previously employed clustering methods, our approach can be used to find units with similar values of $\xi$ and yet different values of $X$ allowing us to identify $\beta_0$. Thus, we also complement the methodology of bonhomme2022discretizing by providing a new grouping approach which can be instrumental in both network and panel models with low-rank regressors.

When $\xi$ is discrete, the considered model (ref) belongs to the class of stochastic block models (SBM) with covariates (e.g., mele2023spectral,ma2022detecting). While our framework and estimation approach are general enough to cover this special setup, in this paper, we focus on the case when unobserved heterogeneity is continuous.\footnote{In fact, when $\xi$ is discrete, agents can be correctly classified into groups with the same values of $\xi$ based on the already introduced pseudo-distance $\hat d_{ij}$, which greatly simplifies the asymptotic analysis; we will discuss this in more detail shortly after we introduce Assumption (ref) in Section (ref).} In a recent study, kitamura2024estimating study a nonseparable nonparametric variation of SBM with covariates. Since their approach crucially relies on the discreteness of $\xi$ whereas ours exploits separability between the observables and unobservables as in (ref), the frameworks considered in this paper and by kitamura2024estimating are non-nested and complementary.

Another strand of the literature emphasizes the importance of using network data to control for endogeneity in peer effects and other related models (e.g., \citealp*{goldsmith2013social,johnsson2021estimation,auerbach2022identification,starck2025improving}). In these papers, the agents' latent characteristics $\xi$ affect both the individual outcomes of interest as well as the network formation process. To tackle this problem, auerbach2022identification and johnsson2021estimation propose using certain network statistics to identify agents with the same values of $\xi$ allowing the authors to control for the unobservables in the other (cross-sectional) regression of interest (e.g., in a peer effects model). The important difference between their and our work is that we use the same network data both to control for unobserved heterogeneity and to identify the parameters of interest, allowing us to disentangle the effects of observables and unobservables within a single network model.\footnote{For example, similarly to the clustering methods of bonhomme2022discretizing and beyhum2024inference, a direct application of auerbach2022identification's approach to the observed network only allows one to identify agents with the same values of both $\xi$ and $X$ once again precluding identification of $\beta_0$.}

Finally, we emphasize that the considered model (ref) does not incorporate interaction externalities. Specifically, we assume that conditional on the agents' observed and unobserved characteristics, the interaction outcomes are independent. This assumption is plausible when the interactions are primarily bilateral. For excellent recent reviews of econometrics of networks, with and without strategic interactions, we refer the reader to graham2020econometric,GRAHAM2020111,de2020econometric.

{The rest of the paper is organized as follows. In Section (ref), we formally introduce the framework and provide (heuristic) identification arguments for the semiparametric regression model (ref). Section (ref) turns these ideas into estimators of the parameters of interest. In Section (ref), we establish consistency of the proposed estimators and derive their rates of convergence. In Section (ref), we generalize the proposed identification argument to cover more general settings including the nonlinear model (ref) as well as its nonparametric analogue. Section (ref) provides numerical and empirical illustrations. A supplementary appendix contains all proofs, and additional discussions and illustrations.

Identification of the Semiparametric Model

In this section, we introduce the framework and provide a conceptual discussion of identification of the semiparametric regression model. This discussion is supposed to illustrate the anatomy of the model and to highlight the main insights of our identification strategy, and is deliberately not formalized here. In Section (ref), we will turn these ideas into practical estimators. Identification then is formally demonstrated in Section (ref), where establish consistency of our estimators and provide their rates of convergence.

The model

We consider a network consisting of $n$ agents. Each agent $i$ is endowed with characteristics $Z_i = (X_i,\xi_i)$, where $X_i \in \mathcal X$ is observed by the econometrician, while $\xi_i \in \mathcal E$ is not. We consider the following semiparametric regression model, where the (scalar) outcome of the interaction between agents $i$ and $j$ is given by

Here, $w: \mathcal X \times \mathcal X \rightarrow \mathbb R^p$ is a known function, which transforms the observed characteristics of agents $i$ and $j$ into a pair-specific vector of covariates $W_{ij} \coloneqq w(X_i,X_j)$, $\beta_{0} \in \mathbb R^p$ is the parameter of interest, and $\varepsilon_{ij}$ is an unobserved idiosyncratic error. Note that unlike $w(\cdot,\cdot)$, the coupling function $g: \mathcal E \times \mathcal E \rightarrow \mathbb R$ is unknown, and the dimension of the fixed effect $\xi_i \in \mathcal E$ is not specified. For simplicity of exposition, we will suppose that $\xi_i \in \mathbb R^{d_\xi}$ even though, in principle, the same insights apply when $\mathcal E$ is a general metric space.

For concreteness, we will also focus on an undirected model with $Y_{i j} = Y_{j i}$, so $w(\cdot,\cdot)$ and $g(\cdot,\cdot)$ are symmetric functions, and $\varepsilon_{ij} = \varepsilon_{ji}$. The methodology presented in this paper straightforwardly extends to directed networks and general two-way settings including panel models; see Section (ref) for a more detailed discussion.

The following assumption formalizes the sampling process.

Assumption (ref) is standard for the networks literature.\footnote{See, e.g., Graham2017,gao2020nonparametric,gao2023logical.} The sampling process could be thought of as follows. First, the characteristics of agents $\{Z_i\}_{i=1}^n$ are independently drawn from some population distribution. Then, conditional on the drawn characteristics, the idiosyncratic errors $\{\varepsilon_{ij}\}_{i < j}$ are independently drawn from the conditional distributions, which potentially depend on the characteristics of the corresponding agents $Z_i$ and $Z_j$.

Identification of $\beta_0$: main insights

We study identification and estimation of $\beta_0$ under the large network asymptotics, which takes $n \rightarrow \infty$. The identification argument is based on the following observation. Suppose that we can identify two agents $i$ and $j$ with the same unobserved characteristics, i.e., with $\xi_i = \xi_j$. Then, for any third agent $k$, the difference between $Y_{ik}$ and $Y_{jk}$ is given by

The conditional mean independence of the regression errors now guarantees that $\beta_0$ can be identified from the regression of $Y_{ik} - Y_{jk}$ on $W_{ik} - W_{jk}$, provided that we have “enough” variation in $W_{ik} - W_{jk}$. Formally, we have

provided that $\mathbb{E}\left[(W_{ik} - W_{jk})(W_{ik} - W_{jk})'|Z_i,Z_j\right]$ is invertible. Since agents $i$ and $j$ are treated as fixed, the expectations are conditional on their characteristics $Z_i$ and $Z_j$. At the same time, $Z_k$, the characteristics of agent $k$, and the idiosyncratic errors $\varepsilon_{ik}$ and $\varepsilon_{jk}$ are treated as random and integrated over. Note that the invertibility requirement insists on $X_i$ and $X_j$, the observed characteristics of agents $i$ and $j$, to be “sufficiently different”. Indeed, if not only $\xi_i = \xi_j$ but also $X_i = X_j$, this condition is clearly violated since $W_{ik} - W_{jk} = 0$ for any agent $k$: in this case, $\beta_{0}$ cannot be identified from the regression (ref).

Hence, the problem of identification of $\beta_0$ can be reduced to the problem of identification of agents $i$ and $j$ with the same values of the unobserved fixed effects ($\xi_i = \xi_j$) but with “sufficiently different” values of $X_i$ and $X_j$.

Let $Y_{ij}^*$ be the error free part of $Y_{ij}$, i.e.,

Consider the following (squared) pseudo-distance between agents $i$ and $j$

where $\mathcal B \ni \beta_0$ is some parameter space. Here, the expectation is conditional on the characteristics of agents $i$ and $j$ and is taken over $Z_k$. Clearly, $d_{ij}^2 = 0$ when $\xi_i = \xi_j$: in this case, the minimum is achieved at $\beta = \beta_0$. Moreover, under a suitable (rank) condition (which we will formally discuss in Section (ref)), $d_{ij}^2 = 0$ also necessarily implies that $\xi_i = \xi_j$. Consequently, if $d_{ij}^2$ were available, agents with the same values of $\xi$ could be identified based on this pseudo-distance.

However, the expectation (ref) cannot be directly identified, since the error free outcomes $Y_{ij}^*$ are not observed. In the following Sections (ref) and (ref), we will argue that the pseudo-distances $d_{ij}^2$ (or their close analogues) are identified for all pairs of agents $i$ and $j$ and, hence, can be used to identify agents with the same values of $\xi$ (and different values of $X$).

Comparison with the existing approaches to identification

Before we proceed with identification of $d_{ij}^2$, we would like to compare our identification argument with the existing results for (interactive fixed effects) panel and network models and highlight some fundamental differences.

One of the most common ways to ensure identification of $\beta_0$ in panel regression models with interactive fixed effects is to assume that covariates $W_{it}$ have sufficient independent variation over $i$ and $t$. This assumption is also commonly referred to as the “high-rank” regressors condition. It plays a crucial role in establishing consistency and deriving asymptotic properties of various estimators of $\beta_0$ in models with unobserved heterogeneity following a linear factor form (e.g., bai2009panel,moon2015linear,armstrong2022robust) and in settings with more general nonparametric structures of unobserved heterogeneity similar to the one studied in this work bonhomme2022discretizing,freeman2023linear,beyhum2024inference.

In the context of this paper, the “high-rank” regressors condition translates into assuming that $W_{ij} = w(X_i,X_j) + \eta_{ij}$, where $\eta_{ij}$ exhibits sufficient independent variation across dyads, e.g., $\eta_{ij}$'s are (conditionally) independent. Identification of $\beta_0$ then can be ensured using the “exogenous” (“high-rank”) variation in $\eta_{ij}$ orthogonal to the unobserved (low-dimensional or even low-rank) component $g(\xi_i,\xi_j)$. While this strategy is commonly employed in panel models with high-rank regressors, it cannot be applied in network models with $W_{ij}$ typically taking the (low-dimensional) form of $W_{ij} = w(X_i,X_j)$ without the additional exogenous “high-rank” component $\eta_{ij}$.

The absence of the $\eta_{ij}$ component in the model substantially complicates both identification of $\beta_0$ as well as the asymptotic analysis of the subsequently constructed analogue estimators. To see this, notice that if $W_{ij} = w(X_i,X_j) + \eta_{ij}$, it would suffice to match agents with the same values of both $X$ and $\xi$, i.e., with $X_i = X_j$ and $\xi_i = \xi_j$, because in this case $\beta_0$ can still be identified as in (ref) thanks to the remaining variation in $W_{ik} - W_{jk} = \eta_{ik} - \eta_{jk}$. This matching strategy, for example, is employed by beyhum2024inference who generalized the clustering approach of bonhomme2022discretizing to semiparametric panel regressions. However, when $W_{ij} = w(X_i,X_j)$ and the $\eta_{ij}$-component is absent, such matching approaches fail to identify $\beta_0$ due to the lack of the remaining variation in $W_{ik} - W_{jk} = 0$.

In this paper, we overcome this difficulty and establish identification of $\beta_{0}$ by demonstrating how to find agents with the same values of $\xi$ but different values of $X$ allowing us to use the variation in $W_{ik} - W_{ij}$ without relying on the “high-rank” exogenous variation in $\eta_{ij}$ like the panel literature does. The latter is the most important and fundamental difference between our work and the recent panel papers by freeman2023linear and beyhum2024inference. Likewise, our identification approach is new and structurally different from the ones previously employed in the literature, and while we focus on network models in this paper, the proposed methodology can also be instrumental in establishing new identification results in panels with low-rank regressors.

Finally, we want to stress that the proposed approach to finding units with the same values of $\xi$ and different values of $X$ is new and structurally different from the previously employed approaches used to match agents with similar latent characteristics in panel and network models. In particular, our framework features two types of individual specific-characteristics, observed $X_i$ and unobserved $\xi_i$, which are needed to be treated differently: we want to identify the “causal” effect of $X_i$ while keeping $\xi_i$ fixed. At the same time, the existing matching methods, including k-means clustering employed in bonhomme2022discretizing and beyhum2024inference, and the similarity based approach pioneered by Zhang2017 and auerbach2022identification, either drop $X_i$ altogether or effectively match on $Z_i = (X_i,\xi_i)$, which, as explained above, does not allow one to identify $\beta_0$.\footnote{Also, notice that, even if one abstracts from having or explicitly controlling for observed $X_i$ in the studied network model, clustering based on individual-specific moments $h_i$ originally proposed by bonhomme2022discretizing and subsequently employed in beyhum2024inference might fail to match agents with similar values of $\xi$. Specifically, their approach requires availability of (consistently estimable) $h_i = \varphi(\xi_i)$ informative about $\xi_i$ but such individual-specific moments might not exist at all in the network setting. To see this, consider a simplified version of the studied model $Y_{ij} = g(\xi_i,\xi_j) + \varepsilon_{ij}$ without covariates, where (i) $\xi_i$ is uniformly distributed on some sphere in $\mathbb R^{d_\xi}$, (ii) $g(\xi_i,\xi_j) = \mathscr g (\left\Vert \xi_i - \xi_j\right\Vert)$ is effectively determined by $\left\Vert \xi_i - \xi_j\right\Vert$ exclusively and captures homophily based on $\xi$, (iii) $\varepsilon_{ij}$'s are iid and independent from all the $\xi$'s. It is clear that, given the spherical symmetry of this model, any individual-specific moments $h_i$ are completely uninformative about $\xi_i$, rendering the described clustering method inappropriate.}

Identification under conditional homoskedasticity

In this section, we consider the case when the regression errors are homoskedastic, i.e., when

For a pair of agents $i$ and $j$, consider the following conditional expectation

Essentially, $q_{ij}^2$ is a feasible analogue of $d_{ij}^2$ with $Y_{ik}$ and $Y_{jk}$ replacing $Y_{ik}^*$ and $Y_{jk}^*$. Importantly, unlike $d_{ij}^2$, $q_{ij}^2$ is immediately identified and can be estimated by

Notice that since $Y_{ik} = Y_{ik}^* + \varepsilon_{ik}$ and $Y_{jk} = Y_{jk}^* + \varepsilon_{jk}$,

where the second and the third equalities follow from Assumption (ref)(ref). Hence, when the errors are homoskedastic and (ref) holds, we have

Thus, for every pair of agents $i$ and $j$, $q_{ij}^2$ differs from $d_{ij}^2$ by a constant term $2 \sigma^2$.

Imagine that for a fixed agent $i$, we are looking for a match $j$ with the same value of $\xi$. As discussed in Section (ref), such an agent can be identified by minimizing $d_{ij}^2$ over potential matches. Then, (ref) ensures that, in the homoskedastic setting, such an agent can also be identified by minimizing $q_{ij}^2$. Hence, agents with the same values of $\xi$ (and different values of $X$) can be identified based on $q_{ij}^2$, which can be directly estimated.

Identification under general heteroskedasticity

Under general heteroskedasticity of the errors, the identification strategy based on $q_{ij}^2$ no longer guarantees finding agents with the same values of $\xi$. Consider the same process of finding an appropriate match $j$ for a fixed agent $i$. As shown in (ref), $q_{ij}^2$ can be represented as a sum of three components. The first term $d_{ij}^2$, which we will call the signal, identifies agents with the same values of $\xi$. The second term $\mathbb{E}\left[\varepsilon_{ik}^2|X_i,\xi_i\right]$ does not depend on $j$. However, under general heteroskedasticity, the third term $\mathbb{E}[\varepsilon_{jk}^2|X_j,\xi_j]$ depends on $\xi_j$ and distorts the signal. Hence, the identification argument provided in Section (ref) is no longer valid in this case.

In this section, we will address this issue and extend the arguments of Sections (ref) and (ref) to a model with general heteroskedasticity. Specifically, we will (heuristically) argue that the error free outcomes $Y_{ij}^*$ are identified for all pairs of agents $i$ and $j$. As a result, the pseudo-distance $d_{ij}^2$ introduced in (ref) is also identified and can be directly employed to find agents with the same values of $\xi$ (and different values of $X$).

Identification of $Y_{ij}^*$

With $Y_{ij}^*$ and $Y_{ij} = Y_{ij}^* + \varepsilon_{ij}$ collected as entries of $n \times n$ matrices $Y^*$ and $Y$ (with diagonal elements of $Y$ missing), the problem of identification and estimation of $Y^*$ based on its noisy proxy $Y$ can be interpreted as a particular variation of the classic matrix estimation/completion problem. Specifically, it turns out that the considered network model (ref) is an example of the latent space model (see, for example, \citet*{Chatterjee2015} and the references therein). In the standard formulation of the latent space model, the entries of (symmetric) matrix $Y$ have the form of

where $f$ is some (unknown) symmetric function, $Z_1, \ldots, Z_n$ are some latent variables associated with the corresponding rows and columns of $Y$, and the errors $\{\varepsilon_{ij}\}_{i < j}$ are assumed to be (conditionally) independent.\footnote{As noted, for example, in \citet*{bickel2009nonparametric} and \citet*{bickel2011method}, the latent space model is natural in exchangeable settings due to the Aldous-Hoover theorem \citep*{aldous1981representations,hoover1979relations}. For a detailed discussion of this result and other representation theorems for exchangeable random arrays, see, for example, \citet*{kallenberg2005probabilistic} and \citet*{orbanz2015bayesian}.}\footnote{The problem of estimation of $Y_{ij}^* = f (Z_i,Z_j)$ from $Y$ is also known as nonparametric regression without knowing the design \citep*{Gao2015} or blind regression \citep*{li2019nearest}. If $Y_{ij}$ is binary, $Y$ can be interpreted as the adjacency matrix of a random graph. In this case, the function $f(\cdot,\cdot)$ is called a graphon, and this problem is commonly referred to as graphon estimation (see, for example, \citealp*{Gao2015,Klopp2017,Zhang2017}).} Notice, that the studied model fits this general formulation, even though, in our setting, we observe $X_i$ (a subvector of $Z_i$).

It turns out that the particular structure of the latent space model (ref) allows one to construct a consistent estimator of $Y^*$ based on a single measurement $Y$. For example, \citet*{Chatterjee2015,Gao2015,Klopp2017,Zhang2017} construct such estimators and establish their consistency in terms of the mean square error (MSE).

In particular, we build on the estimation strategy of \citet*{Zhang2017} to argue that the error free outcomes $Y_{ij}^*$ are identified for all pairs of agents $i$ and $j$. The proposed identification strategy consists of two main steps. First, we argue that we can identify agents with the same values of (both) $X$ and $\xi$. Then, building on this result, we demonstrate how $Y_{ij}^*$ can be constructively identified.

Step 1: Identification of agents with the same values of $X$ and $\xi$\\ Consider a subpopulation of agents with a fixed value of $X = x$ exclusively. Let $g_x (\xi_i,\xi_j) \coloneqq w(x,x)' \beta_0 + g(\xi_i,\xi_j)$ and $P_{\xi|X}(\xi|x)$ denote the conditional distribution of $\xi$ given $X = x$. In this subpopulation, consider the following (squared) pseudo-distance between agents $i$ and $j$

where the second equality uses Assumption (ref)(ref).

The finite sample analogue of $d_\infty^2$ was originally introduced in \citet*{Zhang2017} in the context of nonparametric graphon estimation. It is also closely related to the so-called similarity distance inducing a weak topology on graphons (e.g., see, lovasz2012large and the references therein).

First, notice that (under weak smoothness conditions) $d_{\infty}^2(i,j;x)$ is directly identified and, if a sample of $n_x$ agents with $X = x$ is available, it can be estimated by

Second, note that $d_\infty^2 (i,j;x) = 0$ implies that

for almost all $\xi_k$. Evaluating (ref) at $\xi_k = \xi_i$ and $\xi_k = \xi_j$ and subtracting the latter from the former, we conclude

Thus, $d_{\infty}^2(i,j;x) = 0$ implies that $g(\xi_i,\cdot)$ and $g(\xi_j,\cdot)$ are the same (in terms of the $L^2$ distance associated with the conditional distribution of $\xi|X=x$).\footnote{Similar arguments are also provided in \citet*{lovasz2012large} and \citet*{auerbach2022identification}.} Finally, since the equivalence of $g(\xi_i, \cdot)$ and $g(\xi_j, \cdot)$ is understood as $\xi_i = \xi_j$, we can use $d_{\infty}^2(i,j;x)$ to identify agents with the same values of both $X$ and $\xi$.

Step 2: Identification of $Y_{ij}^*$\\ Now, being able to identify agents with the same values of $X$ and $\xi$, we can also identify the error free outcome $Y_{ij}^* = w(X_i,X_j)' \beta_0 + g(\xi_i,\xi_j)$ for any pair of agents $i$ and $j$. Specifically, for a fixed agent $i$, we can construct a collection of agents with $X = X_i$ and $\xi = \xi_i$, i.e., $\mathcal N_i \coloneqq \{i': X_{i'} = X_i, \xi_{i'} = \xi_i \}$. Similarly, we construct $\mathcal N_j \coloneqq \{j': X_{j'} = X_j, \xi_{j'} = \xi_j\}$. Then,

where $n_i$ and $n_j$ denote the number of elements in $\mathcal N_i$ and $\mathcal N_j$, respectively. Since in the population we can construct arbitrarily large $\mathcal N_i$ and $\mathcal N_j$, (ref) implies that $Y_{ij}^*$ is identified.

Identifiability of $Y_{ij}^*$ is a strong result, which, to the best of our knowledge, is new to the econometrics literature on identification of network and, more generally, two-way models. Importantly, it is not due to the specific parametric form or additive separability (in $X$ and $\xi$) of the model (ref). In fact, by essentially the same argument, the error free outcomes $Y_{ij}^* = f(X_i,\xi_i,X_j,\xi_j)$ are also identified in a fully non-separable nonparametric model

The established result implies that for studying identification of such models, the error free outcome $Y_{ij}^*$ can be treated as directly observed. Removing the error part greatly simplifies the analysis and provides a powerful foundation for establishing further identification results in more complicated settings; see Section (ref) for details and examples.

For example, in the particular context of the model (ref), identifiability of $Y_{ij}^*$ implies that the pseudo-distances $d_{ij}^2$ are also identified for all pairs of agents $i$ and $j$. Hence, as discussed in Section (ref), agents with the same values of $\xi$ (and different values of $X$) and, subsequently, $\beta_0$ can be identified based on $d_{ij}^2$.

Estimation of the Semiparametric Model

In this section, we turn the ideas of Section (ref) into an estimation procedure. First, we construct an estimator of $\beta_0$ assuming that some estimator of the pseudo-distances $\hat d_{ij}^2$ is already available for the researcher. Then, we discuss how to construct $\hat d_{ij}^2$ in the homoskedastic and general heteroskedastic settings.

Estimation of $\beta_0$

Suppose that we start with some (uniformly consistent) estimator of the pseudo-distances $d_{ij}^2$ denoted by $\hat d_{ij}^2$. Using $\hat d_{ij}^2$, we construct the following kernel based estimator of $\beta_0$

where $\Delta W_{ijk} \coloneqq W_{ik} - W_{jk}$ and $\Delta Y_{ijk} \coloneqq Y_{ik} - Y_{jk}$, $K: \mathbb R_+ \rightarrow \mathbb R$ is some kernel supported on $[0,1]$, and $h_n$ is a bandwidth, which needs to satisfy $h_n \rightarrow 0$ (and some additional requirements) as $n \rightarrow \infty$. Hereafter, we will also use the notations $\sum_{i < j} \coloneqq \sum_{i,j \in [n], i < j}$ and $\sum_{k \neq i,j} \coloneqq \sum_{k \in [n], k \neq i,j}$, where $[n] = \{1, \dots, n\}$.

As discussed previously, $\beta_0$ can be estimated by the regression of $\Delta Y_{ijk}$ on $\Delta W_{ijk}$ with fixed agents $i$ and $j$ satisfying $\xi_i = \xi_j$, and, consequently, $d_{ij}^2 = 0$; see (ref) and (ref). However, in a finite sample, we are never guaranteed to find a pair of agents with exactly the same values of unobserved characteristics. The proposed estimator $\hat \beta$ addresses this issue: it combines all of the pairwise-difference regressions weighted by $K({\hat d_{ij}^2}/{h_n^2})$. Typically, the smaller $\hat d_{ij}^2$ is, the closer agents $i$ and $j$ appear to be in terms of $\xi_i$ and $\xi_j$, and the higher weight is given to the corresponding pairwise-difference regression. Specifically, we show that, with probability approaching one, only the pairs that satisfy $\left\Vert \xi_i - \xi_j\right\Vert \leqslant \alpha h_n$ are given positive weights, where $\alpha$ is some positive constant. Since $h_n \rightarrow 0$, the quality of these matches increases and the bias introduced by the imperfect matching vanishes as the sample size grows. In Section (ref), we formalize this discussion by providing the necessary regularity conditions and establishing the rate of convergence for $\hat \beta$.

Estimation of $d_{ij}^2$

The kernel based estimator (ref) builds on the estimated pseudo-distances $\{\hat d_{ij}^2\}_{i \neq j}$. In this section, we construct particular estimators of $d_{ij}^2$ for both homoskedastic and general heteroskedastic settings. Their asymptotic properties will be established in Section (ref).

Estimation of $d_{ij}^2$ under conditional homoskedasticity of $\varepsilon_{ij}$

We start with considering the homoskedastic setting. Recall that in this case, the pseudo-distance of interest $d_{ij}^2$ is closely related to another quantity $q_{ij}^2$ defined in (ref). Specifically, according to (ref), $q_{ij}^2 = d_{ij}^2 + 2 \sigma^2$, where $\sigma^2$ stands for the conditional variance of $\varepsilon_{ij}$. Moreover, unlike $d_{ij}^2$, $q_{ij}^2$ can be directly estimated from the raw data as in (ref).

Then, a natural way to estimate $d_{ij}^2$ is to subtract $2 \hat \sigma^2$ from $\hat q_{ij}^2$, where $\hat \sigma^2$ is an estimator $\sigma^2$. One candidate estimator of $\sigma^2$ is given by

Indeed, in large samples, we expect $\min_{i, j \neq i} d_{ij}^2$ to be small since we are likely to find a pair of agents similar in terms of $\xi$. Hence, in large samples, $\min_{i, j \neq i} q_{ij}^2 = \min_{i, j \neq i} d_{ij}^2 + 2 \sigma^2$ is expected to be close to $2 \sigma^2$. Then, $d_{ij}^2$ can be estimated by

Estimation of $d_{ij}^2$ under general heteroskedasticity of $\varepsilon_{ij}$

As suggested by Section (ref), under general heteroskedasticity of the errors, the first step of estimation of $d_{ij}^2$ is to construct an estimator of $Y_{ij}^*$. Once $\hat Y_{ij}^*$ are constructed for all pairs of agents, $d_{ij}^2$ can be estimated by

As pointed out in Section (ref), the problem of estimation of $Y_{ij}^*$ is closely related to the classic matrix (graphon) estimation/completion problem, with a number of candidate estimators available in the literature (see, e.g., Chatterjee2015). The researcher can pick an appropriate estimator of $\hat Y_{ij}^*$ depending on the context.

In this paper, we propose estimating $Y_{ij}^*$ extending the approach of Zhang2017. For simplicity, first, we consider the case when $X$ is discrete and takes finitely many values. Formal statistical guarantees for the proposed estimator and a discussion of the general case are provided in Section (ref).

First, for all pairs of agents $i$ and $j$, we estimate

Then, for any agent $i$, we define its neighborhood $\hat{\mathcal N}_i (n_i)$ as a collection of $n_i$ agents closest to agent $i$ in terms of $\hat d_{\infty}^2$ among all agents with $X = X_i$, i.e.

Also notice that by construction, $i \in \hat{\mathcal N}_i (n_i)$, so agent $i$ is always included in its neighborhood. Essentially, for any agent $i$, its neighborhood $\hat{\mathcal N}_i (n_i)$ is a collection of agents with the same observed and similar unobserved characteristics. Note that since $X$ is discrete and takes finitely many values, we require $X_{i'} = X_i$. Also, note that the number of agents included in the neighborhoods should grow (at a certain rate) as the sample size increases.

Once the neighborhoods are constructed, we estimate $Y_{ij}^*$ by

where, for the ease of notation, we let $Y_{i' j} = 0$ whenever $i' = j$. Note that $\hat Y_{ij}^*$ is also defined for $i = j$: despite $Y_{ii}$ is not observed, we still can estimate the associated error free outcome $Y_{ii}^* \coloneqq w(X_i,X_i) + g(\xi_i,\xi_i)$.

Large Sample Theory

In this section, we formally study the asymptotic properties of the estimators we provided in Section (ref). The following set of basic regularity conditions will be used throughout the rest of the paper.

Conditions (ref) and (ref) are standard. Condition (ref) requires $g(\cdot,\cdot)$ to be (uniformly) Lipschitz continuous. Condition (ref) requires the conditional distribution of the error $\varepsilon_{ij}|X_i,\xi_i,X_j,\xi_j$ to be sub-Gaussian, uniformly over $(X_i,\xi_i,X_j,\xi_j)$. It allows us to invoke certain concentration inequalities and derive rates of uniform convergence.

Rate of convergence for $\hat \beta$

In this section, we provide necessary regularity conditions and establish the rate of convergence for the kernel based estimator $\hat \beta$ introduced in (ref). We will do that assuming that our candidate estimator $\hat d_{ij}^2$ converges to $d_{ij}^2$ (uniformly across all pairs) at a certain rate $R_n \rightarrow \infty$, i.e.,

We will first derive the result treating (ref) as a high-level assumption, and then verify that it holds and characterize $R_n$ for the candidate estimators (ref) and (ref).

We establish consistency and derive a guaranteed rate of convergence of $\hat \beta$ under the following regularity conditions. For simplicity of exposition, we state these conditions for the case when $\xi$ is scalar but the main result of this section still holds for multivariate $\xi$ under minimal appropriate modifications of the assumptions below.

Assumption (ref) describes the properties of the conditional distribution of $\xi|X$. Note that we focus on the case when $\xi|X=x$ is continuously distributed for all $x \in \text{supp}\left(X\right)$ even though our framework straightforwardly allows for the (conditional) distribution of $\xi$ to have point masses or to be discrete. In fact, the asymptotic analysis is substantially simpler in the latter case. Specifically, if $\xi$ is discrete (and takes finitely many values), the agents can be consistently clustered into groups with the same values of $\xi$ based on the same pseudo-distance $\hat d_{ij}^2$. In this case, $\hat \beta$ is asymptotically equivalent to the oracle pairwise-difference estimator using the knowledge of the true cluster membership, and, as a result, it is asymptotically normal and unbiased. Moreover, in this case, $\beta_0$ can also be estimated by the pooled linear regression, which includes additional interactions of the dummy variables for the estimated cluster membership.\footnote{Similar ideas are also explored in mugnier2025simple in the context of grouped panel models.}

Conditions (ref) and (ref) are weak smoothness requirements. The second part of Condition (ref) bounds the probability mass of $\xi|X=x$, for which $f_{\xi|X}(\xi|X)$ is not potentially continuous on a ball $B_\delta(\xi)$. It allows us to control the probability mass of $\xi$ close to the boundary of its support, where $f_{\xi|X}(\xi|X)$ is allowed to be discontinuous.

Assumption (ref) is a collection of identification conditions. Specifically, Condition (ref) is the identification condition for $\beta_0$. It ensures that in a growing sample, it is possible to find a pair of agents $i$ and $j$ such that (i) $X_i$ and $X_j$ are “sufficiently different”, so the minimal eigenvalue $\lambda_{min}(\mathcal C(X_i,X_j)) > \underline \lambda > 0$, (ii) and yet $\xi_i$ and $\xi_j$ are increasingly similar. The latter is guaranteed by $\int f_{\xi|X}(\xi|X_i) f_{\xi|X}(\xi|X_j) d \xi > \underline \delta$, which implies that the conditional supports of $\xi_i|X_i$ and $\xi_j|X_j$ have a non-trivial overlap. Condition (ref) is crucial for establishing consistency of $\hat \beta$.

Condition (ref) ensures that $d_{ij}^2$ is bounded away from zero whenever $\left\vert \xi_i - \xi_j\right\vert$ is. Notice that it also guarantees that agents that are close in terms of $d_{ij}^2$, must also be similar in terms of $\xi$. Hence, Condition (ref) justifies using the pseudo-distance $d_{ij}^2$ for finding agents with similar values of $\xi$ in finite samples. It also can be interpreted as a rank type condition: for fixed agents $i$ and $j$ with $\xi_i \neq \xi_j$, $g(\xi_i,\xi_k) - g(\xi_j,\xi_k)$ cannot be expressed as a linear combination of the components of $W_{ik} - W_{jk}$.

Condition (ref) is a local counterpart of Condition (ref). It says that, as a function of $\xi_j$, $d^2(X_i,X_j,\xi_i,\xi_j)$ can be locally quadratically approximated around $\xi_j = \xi_i$, and the approximation remainder can be uniformly bounded as $O(\left\vert \xi_j - \xi_i\right\vert^3)$. Condition (ref) also explains why we divide $\hat d_{ij}^2$ by $h_n^2$ for computing the kernel weights. Indeed, locally $d_{ij}^2 \propto (\xi_j - \xi_i)^2$, so the bandwidth $h_n$ effectively controls how large $\left\vert \xi_j - \xi_i\right\vert$ can be for the pair of agents $i$ and $j$ to get a positive weight $K(\frac{\hat d_{ij}^2}{h_n^2})$.

Assumption (ref) specifies the properties of the kernel $K$ and the bandwidth $h_n$. Condition (ref) imposes a number of fairly standard restrictions on $K$ including Lipschitz continuity. Condition (ref) restricts the rates at which the bandwidth is allowed to shrink towards zero. The requirement $n h_n / \ln n \rightarrow \infty$ ensures that we have a growing number of potential matches as the sample size increases. Additionally, to get the desired results we need $R_n h_n^2 \rightarrow \infty$: the bandwidth cannot go to zero faster than $R_n^{-1/2}$. This requirement allows us to bound the effect of the sampling variability coming from the first step (estimation of $\{d_{ij}^2\}_{i \neq j}$) on the second step (estimation of $\beta_0$).

Assumption (ref) is a weak smoothness requirement. It guarantees that as a function of $\xi_2$, the difference $g(\xi_1,\xi) - g(\xi_2,\xi)$ can be (locally) linearized around $\xi_2 = \xi_1$ provided that $\xi_2$ is close to $\xi_1$ relative to the distance between $\xi_1$ and $\xi$ (guaranteed by the restrictions $\left\vert \xi_1 - \xi\right\vert > \delta_n$ and $\left\vert \xi_2 - \xi_1\right\vert \leqslant \delta_n$). The goal of introducing these restrictions is to allow for a possibly non-differentiable $g$, e.g., $g(\xi_i,\xi_j) = \kappa \left\vert \xi_i - \xi_j\right\vert$, since such models of latent homophily are common in the literature (e.g., Hoff2002,handcock2007model). We provide an illustration of Assumption (ref) in Appendix (ref).

Theorem (ref) provides a guaranteed rate of convergence of $\hat \beta$ and, thus, formally establishes identification of $\beta_0$. The derived rate consists of four terms. The $h_n^2$ term accounts for the bias due to the imperfect $\xi$-matching. The terms involving $R_n$ are due to the fact that $\{d_{ij}^2\}_{i \neq j}$ are unknown and need to be estimated. Finally, $n^{-1}$ is the standard sampling variability term in a regression with $O(n^2)$ observations.

In Section (ref), we will derive $R_n$ for the considered estimators of $\{d_{ij}^2\}_{i \neq j}$ and complete the characterization of the rate of convergence for $\hat \beta$ provided in (ref). In particular, we will demonstrate that, under certain conditions, the candidate estimators (ref) and (ref) satisfy (ref) with $R_n = (\frac{n}{\ln n})^{1/2}$, or with an even slower-growing $R_n$ depending on $d_\xi$. For such values of $R_n$, (ref) effectively simplifies as

where we also used $\frac{R_n^{-1}}{h_n^2} = o(1)$ implied by Assumption (ref)(ref). In particular, under $h_n \propto R_n^{-1/3}$, $\hat \beta$ is guaranteed to achieve the following rate of convergence

The rate of convergence established by Theorem (ref) is not necessarily optimal and potentially can be improved, especially if additional smoothness conditions are imposed. The main goal of Theorem (ref) is to demonstrate consistency of $\hat \beta$ and thus complement Section (ref) by formally establishing identification of $\beta_0$ under minimal primitive conditions.

Finally, we want to highlight again two important features of the setting differentiating it from the previously studied frameworks and complicating the asymptotic analysis. First, as discussed in more detail in Section (ref), we consider low-dimensional (or even low-rank) regressors $W_{ij} = w(X_i,X_j)$ typical for network models instead of relying on “high-rank” regressors like most of the related literature does. Second, we only impose minimal smoothness assumptions on $g(\cdot, \cdot)$ and do not require it to be differentiable whereas the alternative approaches building low-rank approximations of interactive unobserved heterogeneity rely on the existence of multiple continuous derivatives of $g(\cdot,\cdot)$ (e.g., fernandez2021low,freeman2023linear,beyhum2024inference); see also Remark (ref) below for a clarification of the role of Assumption (ref).

As discussed earlier, these features are highly representative of network models. Moreover, their unique combination makes the asymptotic analysis highly nonstandard, potentially invalidating the previously proposed methods and established statistical guarantees. For this reason, in this paper, we do not attempt to construct an asymptotically normal estimator of $\beta_0$ or provide an inference method at the cost of imposing additional assumptions restrictive in network models. Instead, we focus on establishing identification of $\beta_0$ by providing a new estimator and showing its consistency in the general setting of the paper.

Rates of uniform convergence for $\hat d_{ij}^2$

In this section, we complement the result of Theorem (ref) by characterizing $R_n$, the rate of uniform convergence for $\hat d_{ij}^2$ defined in (ref), for the candidate estimators (ref) and (ref) in the homoskedastic and then in the general heteroskedastic settings.

Homoskedastic model

First, we consider the homoskedastic case, i.e., we assume that the idiosyncratic errors satisfy (ref). As discussed in Section (ref), in this case, the suggested estimator is given by $\hat d_{ij}^2 = \hat q_{ij}^2 - 2 \hat \sigma^2$, and $d_{ij}^2 = q_{ij}^2 - 2 \sigma^2$ (see (ref) and (ref), respectively). Thus,

This bound, together with the following lemma, allows us to characterize $R_n$.

The bounds (ref) and (ref) established by Lemma (ref) together with (ref) allow us to provide $R_n$ for the estimator (ref) in the homoskedastic setting. To make this characterization complete, we also need to provide a bound on $\left\Vert \xi_i - \xi_j\right\Vert^2$ in (ref).

In particular, since $\mathcal E$ is bounded (Assumption (ref)(ref)), we can guarantee that

for some $C > 0$. While this bound is conservative, it ensures that, for $d_\xi \leq 4$, the contribution of $\left\Vert \xi_i - \xi_j\right\Vert^2$ is negligible, resulting in the following corollary of Lemma (ref).

Corollary (ref) ensures that when the errors are homoskedastic and $d_\xi \leqslant 4$, $\hat d_{ij}^2$ given by (ref) satisfies (ref) with $R_n = \left(\frac{n}{\ln n}\right)^{1/2}$. Hence, when $h_n \propto R_n^{-1/3} = \left(\frac{\ln n}{n}\right)^{-1/6}$, (ref) implies that

so the guaranteed rate of convergence for $\hat \beta$ is $\left(\frac{n}{\ln n}\right)^{1/3}$.

Model with general heteroskedasticity

In this section, we establish the rate of uniform convergence for $\hat d_{ij}^2$ under general heteroskedasticity of the errors. First, we suppose that $X$ is discrete and derive $R_n$ for the estimator given by (ref)-(ref). Then, we discuss how this estimator can be modified to accommodate continuously distributed $X$. Finally, we also provide a general result establishing $R_n$ for $\hat d_{ij}^2$ based on a generic estimator of $Y_{ij}^*$, potentially other than (ref).

Estimation of $d_{ij}^2$ when $X$ is discrete\\ Now we formally derive the rate of uniform convergence for $\hat d_{ij}^2$ given by (ref)-(ref). This rate crucially depends on the asymptotic properties of the denoising estimator $\hat Y_{ij}^*$.

To formally state the asymptotic properties of $\hat Y_{ij}^*$, we introduce the following assumption.

As pointed out before, we suppose that $X$ is discrete and takes finitely many values. Condition (ref) is a weak condition imposed on the conditional distribution of $\xi|X$. If $\xi|X$ is continuously distributed, it is satisfied when the conditional density $f_{\xi|X}(\xi|x)$ is (uniformly) bounded away from zero and its support is not “too irregular”.

For any matrix $A \in \mathbb R^{n \times n}$, let $\left\Vert A\right\Vert_{2,\infty} \coloneqq \max_{i} \sqrt{\sum_{j=1}^n A_{ij}^2}$. Also let $\hat Y^*$ and $Y^*$ denote $n \times n$ matrices with entries given by $\hat Y_{ij}^*$ and $Y_{ij}^*$.

Theorem (ref) establishes two important asymptotic properties of $\hat Y_{ij}^*$. In fact, both results play key roles in bounding $R_n$, the rate of uniform convergence for $\hat d_{ij}^2$.

Part (ref) is analogous to the result of Zhang2017. While Zhang2017 only consider binary outcomes and do not allow for observed covariates, we extend their result by allows for (i) $d_\xi > 1$, (ii) possibly non-binary outcomes and unbounded idiosyncratic errors, (iii) observed (discrete) covariates $X$. Part (ref) is new. It allows us to substantially improve on $R_n$ compared to what Part (ref) can guarantee individually.\footnote{See Lemma (ref) for the comparison of the rates.}

Note that Theorem (ref) requires $n_i$, the number of agents included into $\hat{\mathcal N}_i (n_i)$, to grow at $(n \ln n)^{1/2}$ rate. As shown in Zhang2017, this rate is optimal.\footnote{The optimal choice of $n_i$ remains the same for $d_\xi > 1$.} In applications, the authors also recommend taking $n_i = C (n \ln n)^{1/2}$ with $C \simeq 1$ and document robustness of their numerical results to the choice of $C$.

Building on Theorem (ref), we now provide rate of uniform convergence for $\hat d_{ij}^2$.

Theorem (ref) establishes the rate of uniform convergence for the proposed estimator of $d_{ij}^2$. Specifically, it ensures that, for the considered $\hat d_{ij}^2$, (ref) holds with $R_n = \left(\frac{n}{\ln n}\right)^{\frac{1}{2 d_\xi}}$. Combined with (ref), this guarantees that under the appropriate choice of the bandwidth $h_n$, we have

Thus, when $X$ is discrete, $\beta_0$ can be estimated (at least) at $\left(\frac{n}{\ln n}\right)^{\frac{1}{3 d_\xi}}$ rate. If $\xi$ is scalar, the guaranteed rate of convergence for $\hat \beta$ is the same as in the homoskedastic case.

Estimation of $d_{ij}^2$ when $X$ is continuously distributed\\ The proposed estimator $\hat Y_{ij}^*$ given by (ref)-(ref) can be straightforwardly modified to allow for continuously distributed (components of) $X$. Since, in this case, the probability of finding two agents with exactly the same values of $X$ is zero, we have to modify the construction of $\hat{\mathcal N}_i (n_i)$ previously provided in (ref). One natural possibility is to consider

The parameter $\delta_n$ controls the quality of matching based on $X$. To guarantee consistency of $\hat Y_{ij}^*$, we will need $\delta_n$ to converge to zero but slowly enough to ensure that we can still find a growing number of good matches increasingly similar in terms of both $X$ and $\xi$. Once the constructed neighborhoods are appropriately modified, $\hat Y_{ij}^*$ can still be computed as in (ref), and, with some work, the result of Theorem (ref) can be generalized accordingly to allow for continuously distributed covariates.

Estimation of $d_{ij}^2$ using general matrix denoising techniques\\ Theorem (ref) establishes uniform consistency of $\hat d_{ij}^2$ leveraging the specific structure of $\hat Y_{ij}^*$ in (ref) and requires $X$ to be discrete. To complement Theorem (ref) and formally establish identification of $d_{ij}^2$ and $\beta_0$ in the general setting, we will now provide a generic consistency result, which does not require $\hat Y_{ij}^*$ to have any particular structure and holds regardless of whether $X$ is discrete or continuously distributed.

Lemma (ref) guarantees that $\hat d_{ij}^2$ in (ref) is uniformly consistent for $d_{ij}^2$ provided that $\hat Y_{ij}^*$ is consistent for $Y_{ij}^*$ in the $(2,\infty)$ norm, and thus it justifies using alternative estimators $\hat Y_{ij}^*$ available in the literature.

Together with the result of Theorem (ref), Lemma (ref) delivers consistency of $\hat \beta$ and thus formally establishes identification of $\beta_0$ in the general setting. In particular, as argued in Remark (ref), if $X$ is continuously distributed, we can still construct $\hat Y^*$ satisfying the requirement of the lemma with $\mathcal R_n = \left(\frac{n}{\ln n}\right)^{\frac{1}{2(d_\xi + d_X)}}$. Thus, Lemma (ref) guarantees that $\beta_0$ can be consistently estimated when $X$ is continuously distributed.

Uniformly consistent estimation of $Y_{ij}^*$

One of the contributions of this paper is establishing identification of the error free outcomes $Y_{ij}^*$'s. In Section (ref), we heuristically argued that $Y_{ij}^*$ is identified. In this section, we construct a uniformly consistent estimator of $Y_{ij}^*$ and, hence, formally prove its identification.

As before, first, we suppose that $X$ is discrete and takes finitely many values. The estimator we propose is an analogue of $\hat Y_{ij}^*$ given by (ref), which we used before to construct $\hat d_{ij}^2$. It utilizes exactly the same neighborhoods as in (ref) but, unlike $\hat Y_{ij}^*$, averages over all unique outcomes $Y_{i' j'}$ with $i' \in \hat{\mathcal N}_i(n_i)$ and $j' \in \hat{\mathcal N}_j(n_j)$.\footnote{Recall that in the studied undirected model, $Y_{ij} = Y_{ji}$, and $Y_{ij}$ is not observed for $i = j$.} For example, if $\hat{\mathcal N}_i(n_i)$ and $\hat{\mathcal N}_j(n_j)$ have no elements in common, then the proposed estimator takes a simple form as in (ref). More generally, for any $i$ and $j$, let

Essentially, $\hat{\mathcal M}_{ij}$ is a collection of unique unordered pairs of indices from the Cartesian product of $\hat{\mathcal N}_i (n_i)$ and $\hat{\mathcal N}_j (n_j)$. Then, $Y_{ij}^*$ is estimated by

where $m_{ij}$ denotes the number of elements in $\hat{\mathcal M}_{ij}$.

Theorem (ref) demonstrates that $\tilde Y_{ij}^*$ is uniformly consistent for $Y_{ij}^*$ and, consequently, it formally proves that $Y_{ij}^*$ is identified. We also stress that the previously employed estimator $\hat Y_{ij}^*$ is not necessarily uniformly consistent since it averages over only $n_i$ outcomes $Y_{i' j}$. At the same time, $\tilde Y_{ij}^*$ averages over $m_{ij} = O(n_i n_j)$ outcomes, which allows us to establish the desired result.

Another implication of Theorem (ref) is that the pair-specific fixed effects $g(\xi_i,\xi_j)$ can also be consistently estimated and, hence, are identified for all pair of agents $i$ and $j$. Consider

where $\tilde Y_{ij}^*$ is given by (ref). Since we have already demonstrated consistency of $\hat \beta$ and uniform consistency of $\tilde Y_{ij}^*$, $\hat g_{ij}$ is also uniformly consistent for $g_{ij} \coloneqq g(\xi_i,\xi_j)$.

Establishing nonparametric identification of the pair-specific fixed effects $g_{ij}$'s is another contribution of the paper. This result is also of high empirical importance since in certain applications, the fixed effects are the primary object of interest.

Extensions

Identification of single index and nonparametric models

In this section, we extend the identification arguments of Section (ref) to cover a wide range of network models, both semiparametric and nonparametric, beyond the model (ref).

First, recall that, as discussed in Section (ref), the error free outcomes are still identified in the most general analogue of (ref) given by

For example, we can formally establish identification of $Y_{ij}^* = f(X_i,\xi_i,X_j,\xi_j)$ by constructing a version of $\tilde Y_{ij}^*$ introduced in (ref), treating $X$ as unobserved, and and showing its uniform consistency using the result of Theorem (ref).

However, identification of $Y_{ij}^*$, the value of $f(X_i,\xi_i,X_j,\xi_j)$, for any pair of agents $i$ and $j$ is fundamentally different from identification of function $f$. Importantly, $Y_{ij}^*$ is not a causal object and cannot be directly employed in counterfactual analysis. Moreover, since $\xi$ is not observed, function $f$ or any features of it relevant for counterfactual analysis cannot be identified unless some additional structure is imposed on $f$.

For example, in Sections (ref)-(ref), we demonstrated that $\beta_0$ is identified and can be consistently estimated when $f(X_i,\xi_i,X_j,\xi_j) = W_{ij}'\beta_0 + g(\xi_i,\xi_j)$, which allows us to recover the ceteris paribus effect of $W_{ij}$ (or $X_{ij}$) on $Y_{ij}$ in this model. Below, we extend these identification results to more general forms of $f$ covering nonlinear and nonparametric models.

Identification of the semiparametric single index model

One empirically relevant generalization of (ref) is allowing $f$ to have a single index structure

where $F(\cdot)$ is a known invertible link function. Notice that the presence of the link function $F(\cdot)$ ensures that (ref) is flexible enough to cover a wide range of the previously studied nonlinear network models. For example, consider the following network formation model

where $Y_{ij}$ is a binary variable indicating whether agents $i$ and $j$ are connected or not, and $U_{ij}$'s are iid draws from some distribution, e.g., logistic or $N(0,1)$. This model is covered by (ref) with $F(\cdot)$ standing for the CDF of $U_{ij}$. If $F(\cdot) = \exp(\cdot)$, then (ref) generalizes the dyadic Poisson regression model commonly used to analyze trade networks.

The previously developed arguments can be immediately applied to establish identification of $\beta_0$ and to construct an analogue estimator. First, note that since $Y_{ij}^* = F(W_{ij}' \beta_0 + g(\xi_i,\xi_j))$ is identified and $F(\cdot)$ is invertible, we can also identify

Next, since $\mathcal Y_{ij}^*$ are effectively observed, we are back in the additively separable setting previously studied in Sections (ref)-(ref) with $\mathcal Y_{ij}^*$ replacing $Y_{ij}^*$, and all the relevant features of the model such as $\beta_0$ and the pair-specific fixed effects $g_{ij} = g(\xi_i,\xi_j)$ are identified for all $i$ and $j$.

In particular, we can still estimate $\beta_0$ as in (ref) with $\hat {\mathcal Y}_{ij}^*$ replacing $Y_{ij}$, using $\hat d_{ij}^2$ given by (ref) with $\hat {\mathcal Y}_{ij}^*$ replacing $\hat Y_{ij}^*$. In this case, the preliminary step involves denoising the observed outcomes by constructing $\hat Y_{ij}^*$ as before, and then computing $\hat{\mathcal Y}_{ij}^* = F^{-1} (\hat Y_{ij}^*)$.

Finally, we want to highlight the importance of identification of the pair-specific fixed effects $g_{ij} = g(\xi_i,\xi_j)$ in the model (ref). Since $F(\cdot)$ is potentially nonlinear, the knowledge of $\beta_0$ alone is not sufficient for identifying some policy relevant quantities such as partial effects. However, since the fixed effects $g_{ij}$'s are identified for all pairs of agents, we can also identify both pair-specific and average partial effects as well as other policy relevant counterfactuals.

Identification of the nonparametric model

In the previously considered models (ref) and (ref), the contribution of observables is parameterized by $\beta_0 \in \mathbb R^p$, which was the primary object of interest in our analysis. In this section, we consider a nonparametric version of (ref) with $f$ given by

where both $\mathscr h: \mathcal X \times \mathcal X \rightarrow \mathbb R$ and $g: \mathcal E \times \mathcal E \rightarrow \mathbb R$ are unknown symmetric functions.

In Sections (ref)-(ref), we considered a special case of (ref) with $\mathscr{h} (X_i,X_j) = w(X_i,X_j)' \beta_0$ and established identification of $\beta_0$. In this section, we will show that, in the studied setting, function $\mathscr{h}(\cdot, \cdot)$ is nonparametrically identified. Importantly, nonparametric identification of $\mathscr{h} (\cdot,\cdot)$ implies that the previously obtained identification results were not driven by the parametric restrictions or linearity imposed on $\mathscr{h} (\cdot, \cdot)$, justifying our focus on the semiparametric model (ref) as a practical approximation of the general nonparametrically identified model (ref).

We establish identification of $\mathscr{h} (\cdot, \cdot)$ and $g_{ij}$ in (ref) under the following assumption.

Discussion of Assumption (ref)(ref). The requirement $\mathscr h(x,x) = 0$ is a normalization. Indeed, since $g(\cdot,\cdot)$ and the dimension of $\xi$ are not specified, it is without loss of generality to let ${g(\xi_i,\xi_j) = \alpha_i + \alpha_j + \psi(\theta_i,\theta_j)}$, where $\psi(\cdot,\cdot)$ is symmetric, and $\xi = (\alpha, \theta')'$. Consider

where $\mathscr h(\cdot,\cdot)$ is symmetric. Then, we can construct an observationally equivalent model with

where $\tilde \xi_i = (\tilde \alpha_i, \theta_i')'$ with $\tilde \alpha_i = \alpha_i + \mathscr{h} (X_i,X_i)/2$. Note that $\tilde{\mathscr h}(\cdot,\cdot)$ is symmetric, continuous, and satisfies $\tilde{\mathscr h}(x,x) = 0$ for all $x \in \mathcal X$. Since $\tilde f(X_i,\tilde \xi_i, X_j, \tilde \xi_j) = f(X_i,\xi_i,X_j,\xi_j)$ for any pair of agents $i$ and $j$, the normalized model (ref) is equivalent to the original model (ref). \ensuremath{\blacksquare}

First, we argue that $\mathscr{h} (x, \tilde x)$ is identified for any fixed $x, \tilde x \in \text{supp}\left(X\right)$. Specifically, fix $X_i = x$ and $X_j = \tilde x$, and consider

Note that since $Y_{ik}^*$ and $Y_{jk}^*$ are identified in the general model (ref), $\mathscr{d}_{ij}^2 (x, \tilde x)$ and $\mu_{ij}^*(x, \tilde x)$ are also identified.

$\mathscr{d}_{ij}^2(x, \tilde x)$ is a nonparametric analogue of the previously considered pseudo-distance $d_{ij}^2$. Since we are interested in identifying $\mathscr{h} (x, \tilde x)$, we additionally condition on $(X_i,X_j) = (x, \tilde x)$ and only consider $X_k = x$ and $X_k = \tilde x$. Notice that when $\xi_i = \xi_j$, then we have

where the minimum is achieved at $\mu_{ij}^*(x,\tilde x) = \mathscr{h} (x, \tilde x)$. The following lemma establishes that the converse is also true: if $\mathscr{d}_{ij}(x, \tilde x)^2 = 0$, then we also necessarily have $\mu_{ij}^*(x, \tilde x) = \mathscr{h} (x, \tilde x)$.

Lemma (ref) establishes identification of $\mathscr{h} (x, \tilde x)$ by showing that $\mathscr{h} (x, \tilde x) = \mu^*_{ij}(x, \tilde x)$ for all agents $i$ and $j$, with $X_i = x$ and $X_j = \tilde x$, satisfying $\mathscr{d}_{ij} (x, \tilde x) = 0$; notice that Assumption (ref)(ref) guarantees that such agents exist. Moreover, since the value of $\mu_{ij}^*(x,\tilde x)$ should be the same for all such $i$ and $j$ (and this is true for any fixed $x, \tilde x \in \text{supp}\left(X\right)$), the nonparametric model (ref) as well as the previously considered models including (ref) are overidentified, and, in principle, can be falsified.

Next, notice that since $\mathscr{h} (X_i,X_j)$ is identified for any pair of agents $i$ and $j$, we can also identify the pair-specific fixed effect as $g_{ij} = Y_{ij}^* - \mathscr{h} (X_i,X_j)$. As a result, we conclude that both $\mathscr{h}$ and the pair-specific fixed effects are nonparametrically identified in model (ref).

Finally, consider a nonlinear single index version of (ref)

where $F(\cdot)$ is a known and invertible link function. This model is a nonparametric version of (ref) considered in Section (ref). By first constructing $\mathcal Y_{ij}^* = F^{-1}(Y_{ij}^*) = \mathscr{h} (X_i, X_j) + g(\xi_i,\xi_j)$ and then applying the results of this section, we conclude that $\mathscr{h} $, the fixed effects $g_{ij}$'s, as well as both pair-specific and average partial effects are also identified in (ref).

Incorporating missing outcomes

From the beginning of Section (ref), to simplify the exposition and facilitate the formal analysis, we assumed that $\{Y_{ij}\}_{i \neq j}$ are observed for all pairs of agents $i$ and $j$. While this assumption is standard in the network formation context, where the absence of an interaction between agents $i$ and $j$ is still recorded as observing $Y_{ij} = 0$, in many other applications, interaction outcomes are available only for a limited number of pairs of agents (e.g., in the matched employer-employee setting). Hence, it is important to discuss (i) how to properly adjust the constructed estimators to account for missing outcomes, and (ii) under which conditions the proposed method remains valid in this case. For simplicity, we will stick with considering an undirected network as before. We will discuss directed networks and general two-way settings covering the important matched employer-employee example in Section (ref) below.

Let $D_{ij}$ be a binary variable such that $D_{ij} = 1$ if $Y_{ij}$ is observed and $D_{ij} = 0$ otherwise, with $D$ denoting the resulting adjacency matrix (by construction, $D_{ii} = 0$). Also, let ${\mathcal O_{ij} \coloneqq \{k: D_{ik} = D_{jk} = 1 \}}$ denote a set of agents $k$ such that $Y_{ik}$ and $Y_{jk}$ are observed.

To fix ideas, consider the homoskedastic estimator first. We adjust (ref) and (ref) as

where $\hat d_{ij}^2$ is still computed as in (ref), and $\left\vert \mathcal O_{ij}\right\vert$ denotes the cardinality of $\mathcal O_{ij}$. Notice that to consistently estimate $q_{ij}^2$, we need $\left\vert \mathcal O_{ij}\right\vert \rightarrow \infty$. Hence, in practice one may want to limit their attention to pairs of agents for which $\left\vert \mathcal O_{ij}\right\vert$ is sufficiently large for $\hat q_{ij}^2$ to be informative.

Next, we want to discuss under which conditions the proposed method, with appropriate modifications as described above, remains valid when some interaction outcomes are missing. First, the selection mechanism needs to be exogenous conditional on the observed and unobserved characteristics of agents $\{(X_i,\xi_i)\}_{i=1}^n$, meaning that the $D$ should be (conditionally) independent of the errors $\{\varepsilon_{ij}\}$.\footnote{This assumption is satisfied if agents form connections based on $\{(X_i,\xi_i)\}_{i=1}^n$ but not on the idiosyncratic errors. For example, in a structural model, it can be rationalized if $\varepsilon_{ij}$'s are drawn after the network is formed.} This assumption is standard in the network regression literature including the seminal AKM model of Abowd1999 and subsequent work, assuming exogeneity of the employer-employee network conditional on the firms' and workers' additive fixed effect (and their observed characteristics). While this assumption is widely used in the AKM literature, it is also often criticized for severely restricting the network formation process. We want to stress that, since we allow for a much more general form of unobserved heterogeneity than the additive AKM model, the network exogeneity assumption is substantially less restrictive in our setting. Specifically, since we do not specify the dimensionality of the fixed effects and their role in the model (e.g., we allow for complementarity between the firm and worker fixed effects), our framework allows for a much broader class of selection mechanism compatible with the exogeneity assumption. Thus, rather than seeing network endogeneity as a potential threat to the validity of our approach, we consider it to be a new tool addressing this concern.

Second, for $\hat \beta$ in (ref) to be consistent, we need to have a growing number of pairs $i$ and $j$, for which $\left\vert \mathcal{O}_{ij}\right\vert \rightarrow \infty$ as $n \rightarrow \infty$. Specifically, the requirement $\left\vert \mathcal{O}_{ij}\right\vert \rightarrow \infty$ ensures that we can consistently estimate $\hat d_{ij}^2$ for a growing group of agents. At the same time, a growing pool of potential matches allows us to find increasingly similar agents controlling the bias of $\hat \beta$ and ensuring its consistency. Notice that this requirement still allows the network to be sparse, and that its adequacy can be evaluated in a given application. Moreover, in Section (ref), we will argue that this requirement becomes even less restrictive in general two-way settings and is plausible for many matched employer-employee data sets.

In the general heteroskedastic case, the first estimation step is to construct $\hat Y_{ij}^*$. In fact, recent developments in the matrix completion literature allow one to consistently estimate $Y^*$ (in terms of the MSE), even when the observed matrix $Y$ is sparse (e.g., \citealp*{Chatterjee2015,Klopp2017,li2019nearest}). Once $\hat Y_{ij}^*$ are constructed (for example, using one of the already developed matrix completion techniques), the rest of the estimation procedure remains the same. We provide an appropriate modification of the previously used estimator of $\hat Y_{ij}^*$ and discuss estimation of $\beta_0$ in more detail in Appendix (ref).

Extension to directed networks and two-way models

Finally, the proposed estimation procedure can also be generalized to cover directed networks and, more generally, two-way models. Specifically, consider a general interaction model

where $i \in \mathcal I$ and $k \in \mathcal K$ index senders and receivers, and $\xi_i$ and $\eta_k$ denote the sender and receiver fixed effects. As in Section (ref), we can also allow for missing interactions, with $Y_{ik}$ observed whenever $D_{ik} = 1$.

As before, our identification and estimation strategies are based on finding agents with similar values of the fixed effects. However, the considered interaction model consists of two types of agents, senders and receivers (e.g., firms and workers). As a result, we have the flexibility to decide whether we want to match senders or receivers depending on the context. For example, consider the sender-to-sender approach. To fix ideas, we also focus on the homoskedastic setting (the general heteroskedastic estimator can be constructed in a similar fashion). In this case, the sender-to-sender estimator of $\beta_0$ has exactly the same form as in (ref)-(ref) previously provided in Section (ref)

As discussed in the previous section, for the sender-to-sender estimator to be consistent, we need a growing pool of senders that we can potentially match satisfying the same requirement $\left\vert \mathcal O_{ij}\right\vert \rightarrow \infty$, so we can consistently determine if senders $i$ and $j$ are similar in terms of $\xi$ or not. Also notice that in this case receivers are allowed to participate only in a few interactions. For example, this suggests that in the matched employer-employee setting, the firm-to-firm approach can be appropriate even when the workers' mobility is limited, so long as we have a sufficient number of pairs of firms with a sufficient number of workers moving from one of them to another.

Numerical Evidence

Numerical experiment in a homoskedastic model

In this section, we illustrate the finite sample properties of the proposed estimators. Specifically, we consider the following homoskedastic variation of (ref):

where $\{\varepsilon_{ij}\}_{i < j}$ are independent draws from $N(0,1)$. The true value of the parameter of interest is $\beta_0 = -1$, so the considered model features homophily based on both $X$ and $\xi$.

We study the performance of the following estimators. The first estimator $\hat \beta_{\text{FE}}$ is produced by the standard linear regression with additive fixed effects. The second estimator $\hat \beta$ is the kernel based estimator (ref) with $\hat d_{ij}^2$ computed as in (ref) and using the Epanechnikov kernel.\footnote{The reported results are robust to the kernel choice.} We choose $h_n^2 = 0.9 \min \left\{\hat \sigma_{\hat d^2}, \text{IQR}_{\hat d^2}/1.349\right\} \binom{n}{2}^{-1/5}$ following the standard (kernel density estimation) rule of thumb. Here $\hat \sigma_{\hat d^2}$ and $\text{IQR}_{\hat d^2}$ stand for the standard deviation and the interquartile range of the estimated pseudo-distances $\{\hat d_{ij}^2\}_{i < j}$, and $\binom{n}{2}$ corresponds to the number of the estimated pseudo-distances.\footnote{Note that since the kernel weights in (ref) are $K(\frac{\hat d_{ij}^2}{h_n^2})$, the “effective” kernel density estimation bandwidth applied to $\{\hat d_{ij}^2\}_{i < j}$ is $h_n^2$, not $h_n$.} Finally, we also compute the 1 nearest neighbor pairwise-difference estimator $\hat \beta_{\text{NN1}}$, which, instead of using kernel weights as in (ref), matches every unit $i$ with exactly one unit $j$ closest to it in terms of $\hat d_{ij}^2$.

We simulate the model (ref) for $n \in \{30, 50, 100\}$ and $\rho \in \{0, 0.3 , 0.5, 0.7\}$. The simulated finite sample properties of the considered estimators are reported in Table (ref) below. The number of replications is 10,000. The naive estimator $\hat \beta_{\text{FE}}$ is biased whenever the observed and unobserved characteristics of agents are correlated. The magnitude of this bias increases rapidly as $\rho$ grows. The proposed estimators $\hat \beta$ and $\hat \beta_{\text{NN1}}$ effectively remove the bias even in networks of a moderate size with $n = 30$. Notice that the magnitudes of the bias for $\hat \beta$ and $\hat \beta_{\text{NN}1}$ are approximately the same but the kernel based estimator $\hat \beta$ is consistently less dispersed. This might suggest that in the studied setting, the 1 nearest neighbor estimator $\hat \beta_{\text{NN}1}$ tends to undersmooth. Finally, notice that the proposed estimators $\hat \beta$ and $\hat \beta_{\text{NN}1}$ dominate the naive estimator $\hat \beta_{\text{FE}}$ not only in terms of the bias but also in terms of the standard deviation/IQR even when $\rho = 0$, i.e., when $\hat \beta_{\text{FE}}$ is consistent. Indeed, when the fixed effects contribution to the variability in $Y_{ij}$ is large, controlling for the unobservables (as both $\hat \beta$ and $\hat \beta_{\text{NN1}}$ do by differencing them out) can substantially improve precision even at the cost of significantly reducing the effective sample size.

Empirical illustration: homophily in online social networks

In this section, we illustrate the usefulness of our method in the context of estimating homophily in online social networks using the Facebook 100 data set traud2012social. This dataset contains Facebook friendship network as well as nodal covariates (gender, major, graduation year, etc.) collected at 100 colleges and universities in the US in 2005.

Specifically, we consider the following logistic network formation model

where $Y_{ij}$ is a binary variable indicating whether students $i$ and $j$ are friend or not, $X_i$ is a gender dummy, and $\{U_{ij}\}_{i < j}$ are iid draws from the logistic distribution. In this model, $\beta_0$ captures gender homophily. Notice that (ref) can be represented in the regression form as

where $\Lambda(\cdot)$ stands for the logistic CDF. Note that in the studied network formation model the conditional mean is nonlinear and the errors $\varepsilon_{ij}$'s are heteroskedastic.

We estimate $\beta_0$ using the following three estimators. The first is $\hat \beta_{\text{MLE}}$, the naive MLE estimator that ignores unobserved heterogeneity. The second is $\hat \beta_{\text{TL}}$, the tetra-logit estimator introduced in Graham2017 allowing for additive fixed effects, i.e., assuming that ${g(\xi_i,\xi_j) = \xi_i + \xi_j}$. Finally, we also construct the kernel estimator $\hat \beta$ following the procedure described in Section (ref). Specifically, we use $n_i \approx 0.5 (n \ln n)^{1/2}$ for all $i$ for constructing $Y_{ij}^*$, and we choose $K$ and $h_n^2$ as described in Section (ref).

We report results for Princeton University, the class of 2004 (the results are qualitatively similar across the universities and cohorts). This network consists of 541 students with the mean degree of 37.6, i.e., on average the students in this sample are connected with about 7% of the whole class, so the studied network is characteristically sparse.

The results are reported in Table (ref) below. Both the MLE and tetra-logit homophily estimates are much higher compared to the one produced by $\hat \beta$. This is in line with one of the well known perils of estimating homophily effects: naive methods are likely to overestimate homophily associated with observables when agents also exhibit homophily based on their latent characteristics. For example, if students meet their friends in classrooms and gender is predictive of their choices of classes, naive methods might misattribute this effect to gender homophily.

However, since we do not have a confidence interval available for $\hat \beta$, it is not immediately clear if the difference in the produced estimates is systematic or simply attributed to the large sampling uncertainty. In order to address this concern and study the performance of our method in an empirically relevant setting, we perform the following numerical experiment designed to mimic the studied application. Specifically, we consider a version of (ref)

where $(X_i,\xi_i,V_i)$ are iid. We choose $(\alpha_0, \beta_0, \kappa_0, \pi_0) = (-1.5, 0.1, 2, 0.2)$ and $n=550$ to match the main features of the data including the average node degree, its standard deviation, as well as the estimates produced by the methods. Note that this model features homophily based on both $X$ and $\xi$, and ignoring the latter would result in overestimation of $\beta_0$.

The simulation results are also reported in Table (ref) below. First, as expected, $\hat \beta_{\text{MLE}}$ and $\hat \beta_{\text{TL}}$ indeed severely overestimate $\beta_0$, and their biases are substantially larger than their standard deviations. At the same time our estimator $\hat \beta$ does not suffer from this bias and correctly estimates the true effect. Moreover, we also find that, in the studied setting, its standard deviation is comparable to the standard deviations of the other estimators.

In summation, this numerical experiment demonstrates that our method can perform well in an empirically relevant setting featuring nonlinearity of the regression function, heteroskedasticity of the errors, and sparsity representative of real world social networks. Its results also support our empirical finding documenting that the standard approaches overestimate the gender homophily effects and suggest that the difference in the estimates is systematic and not likely to be (entirely) attributed to the sampling variability.