The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
120,430 characters
Estimation of Peer Effects in Endogenous Social Networks: Control Function Approach
\begin{abstract}
We propose methods of estimating the linear-in-means model of peer effects in which the peer group, defined by a social network, is endogenous in the outcome equation for peer effects. Endogeneity is due to unobservable individual characteristics that influence both link formation in the network and the outcome of interest. We propose two estimators of the peer effect equation that control for the endogeneity of the social connections using a control function approach.
We leave the functional form of the control function unspecified and treat it as unknown. To estimate the model, we use a sieve semiparametric approach, and we establish asymptotics of the semiparametric estimator.
\\
{\sc Keywords: peer effects, endogenous network, sieve estimation, control function}
\\
{\sc JEL Classification: C14, C21}
\end{abstract}
\maketitle
\section{Introduction}
The ways in which interconnected individuals influence each other are usually referred to as peer effects. One of the first to formally model peer effects is \cite{Manski1993a}. He proposes the linear-in-means model, in which an individual's action depends on the average action of other individuals and possibly also on their average characteristics. \cite{Manski1993a} assumes that all individuals within a given group are connected. Later literature allows for more complex patterns of connections, in which an individual might be directly influenced by a subset of the group. Examples are \cite{Bramoulle2009}, \cite{Lee2010}, \cite{Lee2007} among others. Models of peer effects have been applied in various areas, such as education, health and development. Examples of applications are found in recent review papers such as \cite{Blume2010}, \cite{Manski2000}, \cite{Epple2011}, \cite{brock2001interactions} and \cite{Graham2011}.
Many models considered in earlier literature assume that connections between individuals are independent of unobserved individual characteristics that influence outcomes. However, assuming exogeneity of the network or peer group is restrictive in many applications. For example,
consider the following widely studied empirical application of peer effects: peer influence on scholarly achievement. The assumption that friendships are exogenous in the outcome equation for scholarly achievement means that there are no unobserved variables that influence both friendship formation and individual grades. However, even if a study controls for observable individual characteristics such as gender, age, race and parents' education, it is likely to omit factors that influence both students' choice of friends and their GPA; for example parental expectations, psychological disorders, or non-reported substance use. For more examples of endogenous peer groups see \cite{brock2001interactions}, \cite{Weinberg2007}, \cite{Shalizi2012} and \cite{hsieh2016social}, among others.
In this paper we propose a method for estimating a linear-in-means model of peer effects, where the peer group is defined by a network that is endogenous in the outcome equation. Our model allows for correlation between the unobserved individual heterogeneity that impacts network formation and the unobserved characteristics of the outcome. For this, we use a dyadic network formation model that allows the unobserved individual attributes of two different agents to influence link formation, and in which links are pairwise independent conditional on the observed and unobserved individual attributes. The network formation we consider in the paper is dense \label{dense1} and nonparametric.
The main contributions of the paper are methodological. First, given the endogenous peer group formation, we show that we can identify the peer effects by controlling the unobserved individual heterogeneity of the network formation equation. Second, we propose an empirically tractable implementation of the control function, whose functional form is not parametrically specified. For this, we propose two approaches, one based on an estimator of the unobserved individual heterogeneity and the other one based on the average node degrees of the network.\footnote{We acknowledge that this approach is developed based on an idea provided by one of the referees. We thank the referee.} Our estimation method is semiparametric because we do not restrict the functional form of the control function. Finally, we derive the limiting distributions of the estimators within a large single network. The main challenge of the asymptotics is handling the strong dependence of observables caused by the dense network.
Other peer effects papers that have considered endogenously formed peer groups and have controlled the endogeneity via various control functions include \cite{GoldsmithP2013}, \cite{hsieh2016social}, \cite{Qu2015}, \cite{Arduini2015a} and \cite{Auerbach2016}. We provide more detail on these papers in Section \ref{sec:related.literature}.
The remainder of the paper is organized as follows. In Section \ref{section: main.idea} we present a high level description of our approach and provide intuition as to its empirical applications.
In Section \ref{section: model of peer effects with endogenous network} we formally present our model. In Section \ref{section: identification} we show {\color{black} how to identify peer effects using control functions.} Estimation is discussed in Section \ref{section: estimation}, and in
Section \ref{subsection: estimation, limiting distributin of estimator} we discuss the limiting distribution of the estimator and propose standard errors.
In Section \ref{section: monte carlo} we present results of Monte Carlo simulations. There we compare the finite sample performance of our two semiparametric estimators against an estimator that assumes unobserved characteristics enter in a linear way, as well as an instrumental variables (IV) estimator that does not control for network endogeneity. {\color{black} We investigate both high degree and low degree networks.}
Section \ref{section: conclusions} concludes.
A word on notation: in what follows we denote scalars by lowercase letters, vectors by lowercase bold letters, and matrices by uppercase bold letters.
\section{Main Idea}\label{section: main.idea}
In this section we introduce a simple model in order to illustrate the main points of our approach. A more general model and detailed discussion of the model will follow later.
\subsection{Simple Model}
A simple peer effect model for the purpose of illustration of the main idea is
\begin{equation}\label{eq: outcome simplified}
y_i = \beta^0 \left( \frac{\sum_{j \neq i}d_{ij}x_j}{\sum_{j\neq i}d_{ij}} \right) + v_i, \quad i = 1,...,N,
\end{equation}
where $x_i$ is a measure of observable characteristics of individual $i$ and $d_{ij}$ is an indicator of individual $i$'s peer, so $d_{ij}=1$ if $i$ and $j$ are directly linked and $0$ otherwise. In (\ref{eq: outcome simplified}), the regressor of interest is the average of the characteristics of those individuals who are linked with $i$, $\frac{\sum_{j \neq i}d_{ij}x_j}{\sum_{j\neq i}d_{ij}}$. For simplicity, we assume that $x_i$ is exogenous with respect to all the unobserved components of the model; this will be relaxed later.
For the link formation, we consider the following dyadic network formation model,
\begin{equation}\label{eq: network simplified}
d_{ij} = \mathbb{I}(g(a_i,a_j)\geq u_{ij}) \mathbb{I}(i \neq j),
\end{equation}
where $a_i$ and $a_j$ are unobserved individual specific characteristics, $u_{ij}$ is a link specific component, and $g(\cdot,\cdot)$ is some function. It should be noted that this model of network formation does not allow for network effects in link formation, as a link between $i$ and $j$ only depends on the characteristics of $i$ and $j$.
The unobserved individual characteristic $a_i$ can be interpreted as social capital that increases the likelihood of forming a link. Depending on the context this could be factors like trustworthiness, socioeconomic status, or outspokenness.
For example, \cite{fafchamps2011} measure the risk sharing links between households in Tanzania and they construct links between households based on the question whom individuals could ``personally rely on for help.'' \cite{fafchamps2007risk} examine the formation of risk-sharing networks using data from the rural Philippines. \cite{Banerjee2013} examine how participation in micro-finance diffuses through a social network which they measure using lending and trust. In these settings, we can think of $a_i$ as a measure of individual trustworthiness and integrity in financial matters.
\cite{Ductor2011} analyze whether knowledge of a researcher's co-authorship network is helpful in predicting his or her productivity. In this setting $a_i$ can be interpreted as some unobserved productivity trait that induces the researcher to have more coauthors, and also to be more productive at writing papers.
\subsection{Control Function and Its Implementation}
The key feature of the peer effect model (\ref{eq: outcome simplified}) and (\ref{eq: network simplified}) is that individual $i$'s unobserved characteristic $a_i$, which impacts link formation, is correlated with $v_i$, $i$'s unobserved characteristic that affects the outcome $y_i$. For example, $a_i$ could be an unobserved component that affects a researcher's publication rate $y_i$, and also his or her co-authorship relationships, $d_{ij}$. Alternatively, we can think of a situation where there are two
types of agents: popular and unpopular. The popular agents are
more likely to be friends with other agents, and popular agents have better
outcomes even in the absence of a peer effect. Then the peer formation $d_{ij}$ becomes correlated with the unobserved component $v_i$ of the outcome, and, as a consequence, the regressor of the peer effect, $\frac{\sum_{j \neq i}d_{ij}x_j}{\sum_{j\neq i}d_{ij}}$, becomes endogenous.
In this paper we use a control function method to handle the endogenous peer group problem. Let $\mathbf{D}_N$ be the $N \times N$ adjacency matrix that describes the network links $d_{ij}$. Suppose that the unobserved characteristics $(a_i,v_i)$ and $u_{ij}$ are randomly drawn over $i$ and $(i,j)$, respectively. Also assume that $u_{ij}$ is independent of $(a_i,v_i)$. Then, for any $i \neq j$, the link $d_{ij} = \mathbb{I}(g(a_i,a_j)\geq u_{ij})$ and $v_i$ are dependent only through $a_i$. Therefore, controlling for $a_i$, the network $\mathbf{D}_N$ and $v_i$ become mean independent, that is,
\begin{align*}
\mathbb{E}(v_i \,|\, \mathbf{D}_N,a_i) = \mathbb{E}(v_i \,|\, a_i) =: h(a_i).
\end{align*}
Suppose that we observe $a_i$. Consider the outcome equation which controls for $a_i$ nonparametrically,
\[
y_i = \beta^0 \left( \frac{\sum_{j \neq i}d_{ij}x_j}{\sum_{j\neq i}d_{ij}} \right) + h(a_i) + \varepsilon_i,
\]
where $\varepsilon_i := v_i - h(a_i)$. Once we control the endogeneity of the network with $a_i$, then the regressor of the peer effect becomes exogenous, and we can estimate the peer effect coefficient $\beta^0$ using the conventional partially linear regression estimation method (e.g. \cite{Robinson1988}).
However, in most empirical applications, $a_i$ is not observed. Then the question becomes how to implement the control function. In this paper, as the main methodological contribution, we propose the following two procedures. Both procedures are implemented with a single snapshot of an observed network.\label{remark-snapshot}
\begin{enumerate}
\item[(i)] First, suppose that $a_i$ can be consistently estimated. An example can be found in \cite{Graham2017} with the specification $g(a_i,a_j) = a_i + a_j$. Then, we estimate $\beta^0$ by running the partially linear regression of $y_i$ on $ \frac{\sum_{j \neq i}d_{ij}x_j}{\sum_{j\neq i}d_{ij}}$ and $h(\widehat{a}_i)$ as in \cite{Robinson1988}.
\item[(ii)] The second method is to use an observed control function that asymptotically carries the same information as $a_i$.
For this, first notice by the WLLN,
\begin{align*}
{\rm deg}_i
&:= \frac{1}{N} \sum_{j \neq i} d_{ij} = \frac{1}{N} \sum_{j \neq i} \mathbb{I}(g(a_i,a_j) \geq u_{ij}) \rightarrow_p \mathbb{P}(d_{ij} = 1 \,|\, a_i).
\end{align*}
Suppose that the network formation probability conditional on $a_i$, $\mathbb{P}(d_{ij} = 1 \,|\, a_i)$, is a monotonic function of $a_i$. A sufficient condition for this is that $g(\cdot, a_j)$ is monotonic in the same direction for all $a_j$, for example\label{remark-g-conditions}
\begin{equation}
g(a_i,a_j) = a_i + a_j - \tau |a_i-a_j| \label{ex.simple.one-to-one}
\end{equation}
with $ 0 \leq \tau < 1$.
In this case, the limit of the average node degree, $ \lim_{N \rightarrow \infty} \frac{1}{N} \sum_{j \neq i} d_{ij}$, carries the same information as the control function $a_i$, which justifies ${\rm deg}_i$ as a proxy of the control function $a_i$, that is, $\mathbb{E}(v_i \,|\, a_i) \simeq \mathbb{E}(v_i \,|\, {\rm deg_i}) =: h_{*}({\rm deg}_i)$. The peer effect coefficient $\beta^0$ can be estimated by using ${\rm deg}_i$ as a control function. More specifically, we estimate $\beta^0$ by running the partially linear regression of $y_i$ on $\frac{\sum_{j \neq i}d_{ij}x_j}{\sum_{j\neq i}d_{ij}}$ and $h_{*}({\rm deg}_i)$.\\
Intuitively, unobserved characteristics $a_i$ drive heterogeneous degree sequences. We can therefore control for degree when estimating peer effects, ignoring the specific choice of a structural model explaining heterogeneous degrees.
\end{enumerate}
The use of degree as a control function requires much fewer restrictions on the specification of the network. Intuitively, the unobserved node (or individual) fixed effects $a_i$ control for heterogeneous degree sequences. Therefore, from an economic point of view, what needs to be controlled is the agent's degree, which validates the control function approach that uses $\rm{deg_i}$. This approach does not require a specification of the specific structural model explaining heterogeneous degree sequences. \label{remark-control-structural-interpretation}
Consistent estimation of $a_i$ usually requires a specific functional form. For example, \cite{Graham2017} assumed an additive model and \cite{ChenFernandez-ValWeidner2018} require an interactive form.
However, there is a disadvantage in the degree approach. The degree approach cannot identify the coefficient of the observed exogenous regressor if the same regressor also impacts the network formation. \label{remark-comparison-approaches-intro}
In Section \ref{section: model of peer effects with endogenous network}, we generalize the simple model (\ref{eq: outcome simplified}) by allowing for an additional peer effect, $\frac{\sum_{j \neq i}d_{ij}y_j}{\sum_{j\neq i}d_{ij}}$, known as the endogenous peer effect, which measures the effects of the outcomes of the peer group on an individual outcome. In this case we have to deal with two kinds of endogeneity in the peer effect regressors: one from the endogenous regressors $y_{j}$ and the other one from the endogenous peers $d_{ij}$.
In Section \ref{section: model of peer effects with endogenous network}, we also generalize the dyadic network formation model by introducing a dyadic component based on observed individual characteristics. We provide application examples of the general model and discuss its features there.
The identification of the peer effects in the general model will be discussed in Section \ref{section: identification}. In Section \ref{section: estimation} we shows how to implement the two aforementioned estimation methods in the general framework.
In the appendix we provide the regularity conditions that are required for the asymptotic results of the paper. All the technical proofs and comprehensive Monte Carlo simulation results are found in the Online Supplement material which is available in \cite{JohnssonMoon2019}.
\subsection{Related Literature} \label{sec:related.literature}
Closely related papers that adopt a control function approach include \cite{GoldsmithP2013}, \cite{hsieh2016social}, \cite{Qu2015}, \cite{Arduini2015a} and \cite{Auerbach2016}. Our paper adopts a frequentist approach based on a nonparametric specification of the network formation, while \cite{GoldsmithP2013} and \cite{hsieh2016social} use the Bayesian method based on a full parametric specification of the network formation and the outcome equation. Like our paper, \cite{Qu2015} assume the network (spatial weights in their model) to be endogenous through unobserved individual heterogeneity. However, our paper is different from \cite{Qu2015} in many aspects. They consider sparse network formation models while we consider a dense network. They restrict the functional form of the control function to be linear, while we impose no restriction on the functional form. The two papers propose different implementations of the control function. Also, in \cite{GoldsmithP2013}, unobserved components account for homophily in link formation, whereas in our setup they mainly drive degree heterogeneity but are allowed to account for homophily as well, as in the example (\ref{ex.simple.one-to-one}).
Our paper is different from \cite{Arduini2015a} regarding the main source of the endogeneity of the network and the form of the control function. \cite{Arduini2015a} assume that the endogeneity of the network is allowed through dependence between the outcome equation error and the idiosyncratic network formation error, like the conventional sample selection model. This model can be interpreted as meeting opportunities being correlated with unobserved ability of the agent that affects the outcome.
\cite{Arduini2015a} consider control functions (both parametric and semiparametric) to deal with the selection bias problem and propose a semiparametric estimator that uses a power series to approximate selectivity bias terms. Regarding asymptotics, in both \cite{Qu2015} and \cite{Arduini2015a}, the asymptotics are derived using near-epoch dependence and are based on the assumption that the number of connections does not increase at the same rate as the square of the network size.
\label{reference-Auerbach}Among the aforementioned related papers, probably the one most closely related to ours is \cite{Auerbach2016}. As a result, we would like to discuss the differences between the two papers in more detail.
The outcome model of \cite{Auerbach2016} is a partially linear regression model where the nonparametric component is an unknown function of the unobserved network heterogeneity,
\begin{align*}
y_i &= \beta^0 x_i + h(a_i) + \varepsilon_i, \\
d_{ij} &= \mathbb{I}(g(a_i,a_j)\geq u_{ij}) \mathbb{I}(i \neq j).
\end{align*}
In the simple peer effect example, the exogenous peer effect corresponds to the regressor $x_i$ above. The network formation is the same as (\ref{eq: network simplified}).
To compare the identification ideas, let's assume that $a_i \sim U[-1/2,1/2]$ and $u_{ij} \sim U[0,1]$. In this case, $d_i := (d_{i1},...,d_{in})'$ and the distribution of $d_i$ of node $i$, whose characteristic is $a_i$, is fully characterized by the link formation probability profile $g(a_i, \bullet)$.
The key condition of \cite{Auerbach2016} is that $h(a_i)$ and the the link formation distribution profile $g_i(\bullet) := g(a_i,\bullet)$ be one-to-one a.s., that is, $g(a,\bullet) \neq g(a^*,\bullet)$ a.s. if and only if $h(a) \neq h(a^*)$. Then, for any distance measure between the two profiles $g_i$ and $g_j$, $d(g_i,g_j)$, it follows that $d(g_i,g_j) = 0$ if and only if $h(a_i) = h(a_j)$.
Based on this, \cite{Auerbach2016} finds that one can control the network endogeneity by pair-wise differencing\footnote{This resembles \cite{Powell1987}, \cite{Heckmanetal1998}, and \cite{AbadieImbens2006}.} of the observations of the two individuals, $i$ and $j$, whose network formation distributions are the same, $d(g_i,g_j) = 0$, and proposes a semiparametric estimator based on matching pairs of agents
with similar columns of the squared adjacency matrix.
Notice that the identification condition of \cite{Auerbach2016} is satisfied if $g(a_i,\bullet)$ and $a_i$ have a one-to-one relation. However, our second identification is based on the condition that $a_i$ and the marginal network probability, $\int g(a_i,\tau) d \tau$, have a one-to-one relation. We admit that this condition is more restrictive than the identification condition of \cite{Auerbach2016}, because our restriction is a special case of his restriction. However, as mentioned in the introduction, our identification under the stronger condition allows for the omitted variable in the peer effects equation to be nonparametrically directly
estimated, which results in the peer effect estimator having the parametric convergence rate ($\sqrt{N}$). This feature is not necessarily guaranteed in the framework of \cite{Auerbach2016}.\footnote{We thank one of the referees for suggesting the comparisons.}
\section{General Model of Peer Effects with an Endogenous Network}\label{section: model of peer effects with endogenous network}
In this section, we introduce a general linear-in-means peer effect model that extends the simple illustrative outcome model with a peer effect in (\ref{eq: outcome simplified}) and the simple dyadic network formation model in (\ref{eq: network simplified}).
\subsection{General Linear-In-Means Peer Effects Model}
As in Section \ref{section: main.idea}, $d_{ij}$ are the observed binary variables that measure undirected links among individuals $i\in \{1,2,\ldots,N\}$.
We assume that individual outcomes are given by the linear-in-means model of peer effects
\begin{equation}\label{model:outcome}
y_i = \left( \sum_{j=1 \atop j\neq i}^Ng_{ij}y_j \right) \beta_1^0 + \mathbf{x}'_{1i}\beta_2^0 + \left(\sum_{j=1\atop j\neq i}^Ng_{ij}\mathbf{x}_{1j}\right)^\prime\beta_3^0+\upsilon_i,
\end{equation}
where $\mathbf{x}_{1i}$ are observed individual characteristics that affect the outcome $y_i$, $v_i$ are unobserved individual characteristics, and
\[
g_{ij} =
\left\{
\begin{array}{cc}
0 & \quad {\rm if } \quad i=j \\
\frac{d_{ij}} {\sum_{j \neq i} d_{ij}} & {\rm otherwise}
\end{array}
\right.
\]
is the weight of the peer effects.
Using the terminology of \cite{Manski1993a}, $\beta_1^0$ captures the endogenous social effect, and $\beta_3^0$ measures the exogenous social effect.
We let $\beta^0 := (\beta^0_1, \beta_2^{0'}, \beta_3^{0'})'$ and denote $\beta = (\beta_1, \beta_2^{'}, \beta_3^{'})'$.
We let $\mathbf{D}_N$ be the $ (N \times N)$ adjacency matrix of the network whose $(i,j)^{th}$ element is $d_{ij}$. We let $d_{ii}=0$ for all $i$, following convention. Let $\mathbf{G}_N$ be the matrix whose $(i,j)^{th}$ element is $g_{ij}$. Recall that $\mathbf{G}_N$ is obtained by row-normalizing $\mathbf{D}_N$.
Denote $\mathbf{X}_{1N}=(\mathbf{x}_{11}',\ldots,\mathbf{x}_{1N}')'$, $\mathbf{y}_N=(y_1,\ldots,y_N)'$ and $\bm{\upsilon}_N=(\upsilon_1,\ldots,\upsilon_N)'$.\label{remark-equilibirum}
{\color{black} Using this notation, we can express the linear-in-means peer effects model (\ref{model:outcome}) as
\begin{equation}
\mathbf{y}_N = \mathbf{G}_N \mathbf{y}_N \beta_1^0 + \mathbf{X}_{1N} \beta_2^0 + \mathbf{G}_N \mathbf{X}_{1N} \beta_3^0 + \bm{\upsilon}_N. \label{model.outcome.matrix}
\end{equation}
Throughout the paper, we assume that $| \beta_1^0 | < 1$.
It is known that when $\mathbf{G}_N$ is row normalized (i.e., $\sum_{j \neq i}g_{ij} = 1$) and $ | \beta_1^0 | < 1$, the (equilibrium) solution of the peer effect model uniquely exists (e.g., see \cite{Bramoulle2009}) as
\begin{align}
\mathbf{y}_N
&= (\mathbf{I}_N-\beta_1^0 \mathbf{G}_N)^{-1}(\mathbf{X}_{1N}\beta_2^0 + \mathbf{G}_N \mathbf{X}_{1N}\beta_3^0 + \bm{\upsilon}_N) \nonumber
\\
&= \sum_{k=0}^{\infty} \left( \beta_1^0 \mathbf{G}_N \right)^k(\mathbf{X}_{1N} \beta_2^0 + \mathbf{G}_N\mathbf{X}_{1N} \beta_3^0 + \bm{\upsilon}_N). \label{model.outcome.reduced.form}
\end{align}
}
In the standard linear-in-means model of peer effects, the main focus has been identification and estimation of peer effects, assuming that the peer group (or the network) is exogenous, that is, $\mathbb{E}[\upsilon_i|\mathbf{X}_{1N},\mathbf{G}_N]=0$. For example, see \cite{Manski1993a} and \cite{Bramoulle2009}, \cite{Lee2007}, and \cite{blume2015linear}.
To identify and estimate the linear-in-means model of peer effects when the peer group is exogenous, it is necessary to take into account the fact that the regressor $\sum_{i=1}^N g_{ij}y_{j}$ is correlated with the error term $\upsilon_i$. For example, if $\upsilon_i\sim\ i.i.d. (0,\sigma^2)$,
it is true that
\begin{equation}
\begin{split}
\mathbb{E}[(\mathbf{G}_N\mathbf{y}_N)'\bm{\upsilon}_N]&=[(\mathbf{G}_N(\mathbf{I}_N-\beta_1^0 \mathbf{G}_N)^{-1}(\mathbf{X}_{1N}\beta_2^0+\mathbf{G}_N\mathbf{X}_{1N}\beta_3^0+\bm{\upsilon}_N))'\bm{\upsilon}_N]\\
&=\mathbb{E}[(\mathbf{G}_N(\mathbf{I}_N-\beta_1^0 \mathbf{G}_N)^{-1}\bm{\upsilon}_N)'\bm{\upsilon}_N]=\sigma_0 tr(\mathbf{G}_N(\mathbf{I}_N-\beta_1^0 \mathbf{G}_N)^{-1})\neq 0.
\end{split}
\end{equation}
To solve this endogeneity problem different estimators have been proposed in the literature, see for example \citet{Kelejian1998}, \cite{Lee2003} and \cite{Lee2007a}. One of the widely used estimation methods is the Instrumental Variables (IV) approach.
\label{remark-IV}
{\color{black} In view of the expression of (\ref{model.outcome.reduced.form}), when $\beta_2^0 \neq 0$, we can use $\mathbf{G}^2_N\mathbf{X}_{1N}$ as the IV of the endogenous regressor $\mathbf{G}_N\mathbf{y}_N$ because $\mathbf{G}^2_N\mathbf{X}_{1N}$ is uncorrelated with $\bm{\upsilon}_N$ while it is correlated with the endogenous regressor $\mathbf{G}_N \mathbf{y}_N$ (see for example \cite{Kelejian1998}, \cite{Lee2003}, and \cite{Bramoulle2009})\footnote{ {\color{black} If $\beta_2^0 = 0$, $\mathbf{y}_N$ does not depend on $\mathbf{X}_{1N}$ and $\mathbf{G}_N^2\mathbf{X}_{1N}$ is not a relevant instrument for $\mathbf{G}_N \mathbf{y}_N$. }}. }Then, the natural estimator is the Two-Stage Least Squares (2SLS) estimator,
\begin{equation}
\widehat{\beta}_N^{2SLS}=(\mathbf{W}_N'\mathbf{Z}_N(\mathbf{Z}_N'\mathbf{Z}_N)^
{-1}\mathbf{Z}_N\mathbf{W}_N)^{-1}\mathbf{W}_N'
\mathbf{Z}_N(\mathbf{Z}_N'\mathbf{Z}_N)^{-1}\mathbf{Z}_N'\mathbf{y}_N,
\end{equation}
where $\mathbf{W}_N=[\mathbf{G}_N\mathbf{y}_N,\ \mathbf{X}_{1N},\ \mathbf{G}_N\mathbf{X}_{1N}]$ and $\mathbf{Z}_N=[\mathbf{X}_{1N},\ \mathbf{G}_N\mathbf{X}_{1N},\ \mathbf{G}^2_N\mathbf{X}_{1N}]$ is the matrix of instruments. For the IVs $\mathbf{Z}_N$ to be strong, we assume that $\beta_2^0 \neq 0$.
When the network matrix is endogenous, $\mathbb{E}[\mathbf{G}_N\bm{\upsilon}_N]\neq 0$, and the procedure used by \cite{Kelejian1998}, \cite{Lee2003}, \cite{Bramoulle2009} and others is no longer valid since the IV matrix $\mathbf{Z}_N=[\mathbf{X}_{1N},\ \mathbf{G}_N\mathbf{X}_{1N},\ \mathbf{G}^2_N\mathbf{X}_{1N}]$ is correlated with the error term $\bm{\upsilon}_N$.
Specifically, the validity of the 2SLS estimator depends on the orthogonality condition $\mathbb{E}[\bm{\upsilon}_N|\mathbf{Z}_N] = 0$, which is implied if $\mathbb{E}[\bm{\upsilon}_N|\mathbf{X}_{1N},\mathbf{G}_N] = 0$. However, it does not hold if the (row normalized) network $\mathbf{G}_N$ is correlated with $\bm{\upsilon}_N$, which is true if unobserved individual characteristics of $\mathbf{G}_N$ directly influence both link formation and individual outcomes.
In this paper, we consider the case where it may be that $\mathbb{E}[\bm{\upsilon}_N|\mathbf{X}_{1N},\mathbf{G}_N]\neq 0$, so that unobserved characteristics that influence link formation can also have a direct effect on individual outcomes. This is an important consideration in many common applications, like the impact of school friendships on scholarly achievement or substance use. Imagine kids from homes where parents help with homework who only form friendships with kids from similar homes. If this unobserved characteristic of parental behavior is not taken into account, and if this is what really determines grades, this effect might falsely be classified as a peer effect.\label{remark-scholarly} {\color{black} A more elaborate discussion of our framework and its empirical applications can be found in Section \ref{section: main.idea}.}
\subsection{Model of Network Formation} \label{sec: model of network formation}
Let $\mathbf{x}_{2i}$ be a vector of observable characteristics of individual $i$, and let $\mathbf{x}_i=\mathbf{x}_{1i}\cup \mathbf{x}_{2i}$.
Define $\mathbf{X}_{2N}$ analogously to $\mathbf{X}_{1N}$ and let $\mathbf{X}_N=\mathbf{X}_{1N}\cup\mathbf{X}_{2N}$.
We introduce $a_i$, a scalar unobserved characteristic of individual $i$, which is treated as an individual fixed effect, and hence, might be correlated with $\mathbf{x}_i$. We denote the vector of individual unobserved characteristics by $\mathbf{a}_N=(a_1,a_2,\ldots,a_N)'$. Individuals are connected by an undirected network $\mathbf{D}_N$, with the $(i,j)^{th}$ element $d_{ij} = 1$ if $i$ and $j$ are directly connected and $0$ otherwise. We assume the network to be undirected\footnote{{\color{black}Our analysis can be extended to the directed network case, but we do not pursue it in this paper.}}, $d_{ij} = d_{ji}$, and assume $d_{ii}=0$ for all $i$, following the convention. In this case, there are $n=\binom{N}{2}$ dyads. Let $\mathbf{t}_{ij}$ denote an $l_T\times 1$ vector of dyad-specific characteristics of dyad $ij$, and we assume that $\mathbf{t}_{ij}=t(\mathbf{x}_{2i},\mathbf{x}_{2j})$.
Agents form links according to
\begin{equation}
d_{ij}=\mathbb{I}( g( t(\mathbf{x}_{2i},\mathbf{x}_{2j}), a_i, a_j ) - u_{ij} \geq 0), \label{model.network.formation}
\end{equation}
where $\mathbb{I}( \bullet)$ is an indicator function. In this setup, link surplus is transferable across directly linked agents and consists of three components: $\mathbf{t}_{ij} := t(\mathbf{x}_{2i},\mathbf{x}_{2j}) $ is a systematic component that varies with observed dyad attributes and accounts for homophily, $a_i$ and $a_j$ account for unobserved dyad attributes (degree heterogeneity), and $u_{ij}$ is an idiosyncratic shock that is i.i.d. across dyads and independent of $\mathbf{t}_{ij}$ and $a_i$ for all $i,j$.
Since links are undirected, the surplus of link $d_{ij}$ must be the same for individual $i$ and $j$. Hence, we assume that the function $t_{ij}$ is symmetric in $i$ and $j$, and the function $g$ is symmetric in $a_i$ and $a_j$.
In the literature, various parametric versions of the network formation in (\ref{model.network.formation}) are used, ({\color{black} see for example \cite{jackson2005survey}, \cite{Graham2017})}). An important example of a parametric specification is the one in \cite{Graham2017},
\begin{equation}
d_{ij}=\mathbb{I}(t(\mathbf{x}_{2i},\mathbf{x}_{2j})'\lambda+a_i+a_j - u_{ij} > 0). \label{model.network.formation.parametric}
\end{equation}
For the purpose of the paper, particularly in constructing the estimators that we introduce in Section \ref{section: estimation}, we do not need a parametric specification.
Regarding the network formation (\ref{model.network.formation}), we impose restrictions (Assumption \ref{assumption:limit.dist} (iii) - (vi) in the Appendix) that imply the following two features. The first feature is that the link formation probability of individual $i$ with characteristics $(\mathbf{x}_{2i},a_i)$ is one-to-one with respect to the unobserved characteristic $a_i$, that is, for all $x_{2i}$,
\begin{equation}
a_i \neq a_i^* \text{ if and only if } \mathbb{P} \left( d_{ij} = 1 \,|\, \mathbf{x}_{2i},a_i \right) \neq \mathbb{P} \left( d_{ij} = 1 \,|\, \mathbf{x}_{2i},a_i^* \right). \label{eq.monotone.link.formation}
\end{equation}
Obviously, this condition is satisfied in the parametric model (\ref{model.network.formation.parametric}). This monotonic condition justifies the use of the average node degree in implementing the control function as introduced in Section \ref{section: main.idea} and will be discussed in Section \ref{sec: estimation with x,a as control}. The second feature is that the network formed by (\ref{model.network.formation}) is dense \label{remark-dense} in the sense that the expected number of connections is proportional to the square of the network size. This is satisfied if the error $u_{ij}$ is drawn randomly from a distribution with full support, while $g( \mathbf{t}_{ij}, a_i, a_j )$ is bounded (see Assumption \ref{assumption:limit.dist} (iii),(iv), and (v) in the Appendix). In this case, the probability of any two individuals forming a link is bounded away from zero and strictly less than one. The dense network model is appropriate for scenarios where any two individuals can plausibly form a link.
Notice that the dense network assumption and the sharing restriction on the net surplus function $g$ are necessary for implementing the control function in Section \ref{section: estimation} and establishing the asymptotic theory of the control function based estimators in Section \ref{subsection: estimation, limiting distributin of estimator}.
If $a_i$ is observed, we can identify and estimate peer effects without these assumptions (see Section \ref{section: identification}).
Regarding the network formation model (\ref{model.network.formation}), it is important to note that the network formation model (\ref{model.network.formation}) rules out interdependent link preferences, and it assumes that links are formed independently conditional on observed individual characteristics and unobserved fixed effects. \label{no network externalities}\label{remark-conditional-independence} As discussed in \cite{Graham2017}, this assumption is appropriate for settings where link formation is driven predominantly by bilateral concerns, such as certain types of friendship networks, trade networks and some models of conflict between nation-states. The model in (\ref{model.network.formation}) is not a good choice when important strategic aspects influence link formation, like when the identity of the nodes to which $j$ is linked influences $i$'s return from forming a link with $j$. A discussion of networks with interdependent links can be found in \cite{Graham2017} and \cite{dePaula2016}. Also, when network externalities are present, the additional complication of multiple equilibria has to be considered, see for example \cite{Sheng2012} for more details.
\section{Identification of peer effects using a control function approach}\label{section: identification}
In this section we provide an identification argument for the peer effect equation based on a control function when the network is endogenous.
\subsection{Control Function of Network Endogeneity}
In this subsection we discuss how to control the endogeneity of the peer group defined by the network formed in equation (\ref{model.network.formation}).
First we introduce a basic assumption that we will maintain throughout the paper.
\begin{assumption}[]\label{as:basic}
(i) $(\mathbf{x}_i,a_i,\upsilon_i)$ are i.i.d. for all $i$, $i=1,\ldots,N$,
(ii) $\{u_{ij}\}_{i,j=1,\ldots,N}$ are independent of $(\mathbf{X}_{N},\mathbf{a}_N, \bm{\upsilon}_N )$ and i.i.d. across $(i,j)$ with cdf $\Phi(\cdot)$, and
(iii) $\mathbb{E}(v_i|\mathbf{x}_i,a_i) = \mathbb{E}(v_i|a_i).$
\end{assumption}
Assumption \ref{as:basic}(i) implies that the observables $\mathbf{x}_i$ and the unobservable characteristics $(a_i,\upsilon_i)$ are randomly drawn. This is a standard assumption in the peer effects literature. Assumption \ref{as:basic}(ii) assumes that the link formation error $u_{ij}$ is orthogonal to all other observables and unobservables in the model. This means that the dyad-specific unobservable shock $u_{ij}$ from the link formation process does not influence outcomes $(y_1,\ldots,y_N)'$. However, we allow for endogeneity of the social interaction group through dependence between the two unobserved components $a_i$ and $\upsilon_i$. This means that the unobserved error $\upsilon_i$ in the outcome equation can be correlated with unobserved individual characteristics $a_i$ that are determinants of link formation. We also allow the observed characteristics $\mathbf{x}_i$ of the outcome equation and the network formation to be correlated with the unobserved components $(\upsilon_i,a_i)$, so that the regressor $\mathbf{x}_{1i}$ can be endogenous in the outcome equation, and the network formation observables $\mathbf{x}_{2i}$ can be arbitrarily correlated with the unobserved individual characteristic $a_i$. In Assumption \ref{as:basic}(iii), we assume that the dependence between $\mathbf{x}_i$ and $\upsilon_i$ exists only through $a_i$. That is, $a_i$ is the fixed effect of individual $i$ and controls the endogeneity of $\mathbf{x}_i$ with respect to $\upsilon_i$.
Notice that the network $\mathbf{D}_N$ defined in (\ref{model.network.formation}) and the (row normalized) network $\mathbf{G}_N$ are measurable functions of $ (\mathbf{x}_{2i},\mathbf{x}_{2,-i},a_i,\mathbf{a}_{-i},\{u_{ij}\}_{i,j=1,\ldots,N}),$ where $\mathbf{x}_{2,-i}=(\mathbf{x}_{2,1},\ldots,\mathbf{x}_{2,i-1},\mathbf{x}_{2,i+1},\ldots,\mathbf{x}_{2,N})$ and $\mathbf{a}_{-i}$ is defined analogously.
Under Assumption \ref{as:basic} we have
\begin{eqnarray}
\mathbb{E}[ \upsilon_i|\mathbf{X}_{N},\mathbf{G}_N,a_i ] &=& \mathbb{E}[\upsilon_i|\mathbf{x}_{-i},\mathbf{G}_N(\mathbf{x}_{2,-i},\mathbf{a}_{-i},\{u_{ij}\}_{i,j=1,\ldots,N},\mathbf{x}_{2i},a_i),\mathbf{x}_{i},a_i] \nonumber \\
&=&\mathbb{E}[\upsilon_i|\mathbf{x}_{i},a_i] \nonumber
= \mathbb{E}[\upsilon_i| a_i],
\label{eq.control.function}
\end{eqnarray}
where the second equality holds because $(\mathbf{x}_{-i},\mathbf{a}_{-i},\{u_{ij}\}_{i,j=1,\ldots,N})$ and $ (\mathbf{x}_{i},a_i,\upsilon_i)$ are independent under Assumptions \ref{as:basic} (i) and (ii). This shows $v_i$ and $(\mathbf{x}_{-i}, \mathbf{G}_N(\mathbf{x}_{2,-i},\mathbf{a}_{-i},\{u_{ij}\}_{i,j=1,\ldots,N},\mathbf{x}_{2i},a_i))$ are mean-independent conditioning on $(\mathbf{x}_{i},a_i)$. The last line follows by the fixed effect assumption, Assumption \ref{as:basic} (iii).
Result (\ref{eq.control.function}) shows that conditional on the unobserved heterogeneity $a_i$ in the network formation (and any subcomponents of $\mathbf{x}_i$), the unobserved characteristic $\upsilon_i$ that affects the outcome $y_i$ becomes uncorrelated with the (row normalized) network $\mathbf{G}_N$ (and the observables $\mathbf{X}_N$). This implies that the network endogeneity can be controlled by $a_i$ (or together with any subcomponents of $\mathbf{x}_i$).
We summarize the discussion above in the following lemma:
\begin{lemma}[Control Function of Peer Group Endogeneity]
\label{lemma:control.function}
Suppose that Assumption \ref{as:basic} holds. Then, $\mathbb{E}[\upsilon_i|\mathbf{X}_{N},\mathbf{G}_N,a_i]=\mathbb{E}[\upsilon_i| \mathbf{x}_{i}, a_i].$
\end{lemma}
\subsection{Identification of Peer Effects with $a_i$ as Control Function}
In this section we show how to identify the peer effects in the outcome question when the endogenous network is formed by (\ref{model.network.formation}). We provide two identification methods depending on whether we control the network (peer group) endogeneity with $a_i$ or $a_i$ together with $\mathbf{x}_{2i}$, in the case when $\mathbf{x}_{2i}$ and $\mathbf{x}_{1i}$ do not overlap.
First notice that regardless of the possible endogeneity of the (row normalized) network $\mathbf{G}_N$, we need to control for the endogeneity of the term $\sum_{j \neq i} g_{ij}y_j$ that represents the so-called endogenous peer effects. When the peer group $\mathbf{G}_N$ is exogenous and uncorrelated with $\upsilon_N$, $\mathbf{G}^2_N \mathbf{X}_{1N}$ is often used as an IV for the endogenous peer effects term $\mathbf{G}_N \mathbf{y}_N$ (See, for example, \cite{Kelejian1998}, \cite{Lee2003}, \cite{Bramoulle2009}.).
Let $\mathbf{Z}_N=[\mathbf{X}_{1N}, \mathbf{G}_N \mathbf{X}_{1N}, \mathbf{G}^2_N\mathbf{X}_{1N} ]$ be the usual IV matrix used in 2SLS estimation of the peer effects equation.
Note that $\mathbf{Z}_N$ is not a valid IV matrix anymore in our framework because the peer group defined by the network $\mathbf{G}_N$ is correlated with $\upsilon_N$ due to potential correlation between the unobserved $\upsilon_i$ and $a_i$.
Let $\mathbf{W}_N=[\mathbf{G}_N \mathbf{y}_N, \mathbf{X}_{1N}, \mathbf{G}_N \mathbf{X}_{1N}]$. Further, denote the transpose of the $i$th row of $\mathbf{Z}_N$ and $\mathbf{W}_N$ by $\mathbf{z}_i$ and $\mathbf{w}_i$, respectively.
Suppose that Assumption \ref{as:basic} holds and so $a_i$ controls the network endogeneity. Then,
\begin{eqnarray}
\mathbb{E}\left[\: \left( \mathbf{z}_i -\mathbb{E}[\mathbf{z}_i|a_i]\right)(\upsilon_i - \mathbb{E}(\upsilon_i|a_i))\: |\:a_i \right]
&=& \mathbb{E}[\mathbf{z}_i\upsilon_i \: | \: a_i] - \mathbb{E}[\mathbf{z}_i \: | \: a_i]\mathbb{E}[\upsilon_i \: | \: a_i] \nonumber \\
&=&\mathbb{E}\left[\mathbb{E}[\mathbf{z}_i\upsilon_i \: | \: a_i,\mathbf{X}_{1N},\mathbf{G}_N] \: | \: a_i\right] - \mathbb{E}[\mathbf{z}_i|a_i]\mathbb{E}[\upsilon_i \: | \: a_i] \nonumber \\
&=&\mathbb{E}\left[\mathbf{z}_i\mathbb{E}[\upsilon_i \: | \: a_i,\mathbf{X}_{1N},\mathbf{G}_N] \: | \: a_i\right] - \mathbb{E}[\mathbf{z}_i \:| \: a_i]\mathbb{E}[\upsilon_i \:| \: a_i] \nonumber \\
&\stackrel{(1)}{=}& \mathbb{E}\left[\mathbf{z}_i\mathbb{E}[\upsilon_i \:| \: a_i] \:|\: a_i\right] - \mathbb{E}[\mathbf{z}_i \:|\: a_i]\mathbb{E}[\upsilon_i \:|\: a_i] \nonumber \\
&=& 0,\label{eq.orthogonality}
\end{eqnarray}
where equality $(1)$ holds by Lemma \ref{lemma:control.function}(a).
This shows that the instrumental variables $\mathbf{z}_{i}$ or $\mathbf{z}_{i} -\mathbb{E}[\mathbf{z}_{i}|a_i]$ become orthogonal to $ \upsilon_i-\mathbb{E}[\upsilon_i|a_i],$ the residual of $\upsilon_i$ after projecting out $a_i$.
Furthermore, if $\mathbb{E}\left[ \left( \mathbf{z}_i -\mathbb{E}[\mathbf{z}_i|a_i] \right)
\left( \mathbf{w}_i - \mathbb{E}[\mathbf{w}_i|a_i]\right)^{\prime}\right]$ has full rank, then we can identify the peer effect coefficients $\beta^0$ as
\begin{eqnarray*}
0 &=& \mathbb{E}\left[ \left( \mathbf{z}_i -\mathbb{E}[\mathbf{z}_i|a_i]\right)
\left(y_i-\mathbf{w}'_i\beta-\mathbb{E}[y_i - \mathbf{w}'_i \beta |a_i] \right) \right] \\
&=& \mathbb{E}[\left( \mathbf{z}_i -\mathbb{E}[\mathbf{z}_i|a_i]\right)(\mathbf{w}_i-\mathbb{E}[\mathbf{w}_i|a_i])'](\beta-\beta^0)+\mathbb{E}[\left( \mathbf{z}_i -\mathbb{E}[\mathbf{z}_i|a_i]\right)(\upsilon_i-\mathbb{E}[\upsilon_i|a_i])]\\
&\stackrel{(1)}{=}&\mathbb{E}[\left( \mathbf{z}_i -\mathbb{E}[\mathbf{z}_i|a_i]\right)(\mathbf{w}_i-\mathbb{E}[\mathbf{w}_i|a_i])'](\beta-\beta^0) \\
& \stackrel{(2)} {\Leftrightarrow} & \beta=\beta^0,
\end{eqnarray*}
where equality $(1)$ follows by the orthogonality result in (\ref{eq.orthogonality}) and equality $(2)$ follows from the full rank condition.
\begin{assumption}[Rank condition]\label{assumption: rank}
$\mathbb{E}\left[ \left( \mathbf{z}_i -\mathbb{E}[\mathbf{z}_i|a_i] \right)
\left( \mathbf{w}_i - \mathbb{E}[\mathbf{w}_i|a_i]\right)^{\prime}\right]$ has full rank.
\end{assumption}
For the full rank condition in Assumption \ref{assumption: rank}, it is necessary that the IVs $\mathbf{z}_i$ and the regressors $\mathbf{w}_i$ have additional variation after projecting out the control function $a_i$. As shown in the Supplementary Appendix \ref{appendix: distribution of est}, when $N$ is large, both $\mathbf{z}_i$ and $\mathbf{w}_i$ become close to functions that depend only on $(\mathbf{x}_i,a_i)$. In this case, for the full rank condition to be satisfied, it is necessary that there be additional random components in $\mathbf{x}_{i}$ that are different from $a_i$, so that the limits of $\mathbf{z}_i$ and $\mathbf{w}_i$ are not linearly dependent. \label{rank - discussion}
As a summary, we have the following first identification theorem.
\begin{theorem}[Identification]Under Assumptions \ref{as:basic} and \ref{assumption: rank}, the parameter $\beta^0$ is identified by the moment condition $\mathbb{E}[\left( \mathbf{z}_i -\mathbb{E} (\mathbf{z}_i|a_i) \right)(y_i - \mathbb{E} (y_i|a_i) - (\mathbf{w}_i-\mathbb{E}(\mathbf{w}_i|a_i))'\beta^0)]=0$:
\[
\mathbb{E}[\left( \mathbf{z}_i -\mathbb{E} (\mathbf{z}_i|a_i) \right)(y_i - \mathbb{E} (y_i|a_i) - (\mathbf{w}_i-\mathbb{E}(\mathbf{w}_i|a_i))'\beta)]=0\ \iff\ \beta=\beta^0.
\] \label{theorem:identification}
\end{theorem}
Theorem \ref{theorem:identification} shows that we can identify the parameter $\beta^0$ by controlling the unobserved network heterogeneity $a_i$ in the outcome equation and taking the residuals $y_i - \mathbb{E} (y_i|a_i) - (\mathbf{w}_i-\mathbb{E}(\mathbf{w}_i|a_i))'\beta$ and using the instrumental variables $\mathbf{z}_i - \mathbb{E}[\mathbf{z}_i|a_i]$.
\subsection{Identification of Peer Effects using $(\mathbf{x}_{2i},a_i)$ as Control Function} \label{section: alternative identification}
In view of the derivation of the control function in (\ref{eq.control.function}) under Assumption \ref{as:basic}, it is possible to use any regressors in $\mathbf{x}_i$ in addition to the unobserved heterogeneity $a_i$.
In this section, we discuss identification of the peer effects using $(\mathbf{x}_{2i},a_i)$ as control function. The reason to consider this particular control function is that we can {\color{black} implement it in the absence of a consistent estimator of $a_i$}, which will be discussed in detail in Section \ref{section: estimation}.
First, suppose that there is no overlap between the regressors in the outcome equation $\mathbf{x}_{1i}$ and the regressors in the network formation equation $\mathbf{x}_{2i}$ and assume the conditions in Assumption \ref{as:basic}.\footnote{Later in this section, we will discuss a more general case where $\mathbf{x}_{1i}$ and $\mathbf{x}_{2i}$ intersect.}
\begin{assumption}[]\label{as:basic.alternative}
Assume that the conditions (i),(ii), and (iii) of Assumption \ref{as:basic} hold. Also, assume that (iv) the explanatory variables in $\mathbf{x}_{1i}$ and $\mathbf{x}_{2i}$ do not overlap (i.e., $\mathbf{x}_{1i} \, \cap \, \mathbf{x}_{2i} = \emptyset$).
\end{assumption}
Then, under Assumption \ref{as:basic} and by (\ref{eq.control.function}), it follows that
\begin{equation}
\mathbb{E}[ \upsilon_i|\mathbf{X}_{N},\mathbf{G}_N,a_i ]
= \mathbb{E}[\upsilon_i|a_i]
= \mathbb{E}[\upsilon_i|\mathbf{x}_{2i},a_i], \label{eq.control.function.alternative}
\end{equation}
where the last line holds by Assumption \ref{as:basic}(iii). Then, similar to (\ref{eq.orthogonality}), we can show that
\begin{eqnarray}
&& \mathbb{E}\left[\: \left( \mathbf{z}_i -\mathbb{E}[\mathbf{z}_i|\mathbf{x}_{2i},a_i]\right)(\upsilon_i - \mathbb{E}(\upsilon_i|\mathbf{x}_{2i},a_i))\: |\: \mathbf{x}_{2i}, a_i \right]
= 0. \label{eq.orthogonality.alternative}
\end{eqnarray}
Furthermore, suppose that the following full rank assumption is satisfied:
\begin{assumption}[Rank condition]\label{assumption:rank.alternative}
$\mathbb{E}\left[ \left( \mathbf{z}_i -\mathbb{E}[\mathbf{z}_i|\mathbf{x}_{2i},a_i] \right)
\left( \mathbf{w}_i - \mathbb{E}[\mathbf{w}_i|\mathbf{x}_{2i},a_i]\right)^{\prime}\right]$ has full rank.
\end{assumption}
Notice that if $\mathbf{x}_{1i}$ and $\mathbf{x}_{2i}$ are overlapped, then the full rank condition in Assumption \ref{assumption:rank.alternative} does not hold.
Using similar arguments that lead to Theorem \ref{theorem:identification}, we can identify the peer effect coefficients $\beta^0$ as
\begin{eqnarray}
0 &=& \mathbb{E}\left[ \left( \mathbf{z}_i -\mathbb{E}[\mathbf{z}_i|\mathbf{x}_{2i},a_i]\right)
\left(y_i-\mathbf{w}'_i\beta-\mathbb{E}[y_i - \mathbf{w}'_i \beta |\mathbf{x}_{2i},a_i] \right) \right]
\Leftrightarrow \beta=\beta^0, \label{eq.identification.alternative}
\end{eqnarray}
This is summarized in the following theorem.
\begin{theorem}[Alternative Identification]Under Assumptions \ref{as:basic}, \ref{as:basic.alternative}, and \ref{assumption:rank.alternative}, the parameter $\beta^0$ is identified by the moment condition
\[
\mathbb{E}[\left( \mathbf{z}_i -\mathbb{E}( \mathbf{z}_i|\mathbf{x}_{2i},a_i )\right)(( y_i - \mathbb{E}(y_i|\mathbf{x}_{2i},a_i) - ( \mathbf{w}'_i -\mathbb{E} (\mathbf{w}_i |\mathbf{x}_{2i},a_i))'\beta]=0\ \iff\ \beta=\beta^0.
\] \label{theorem:identification.alternative}
\end{theorem}
So far, we have considered the case where the regressors $\mathbf{x}_{i1}$ and $\mathbf{x}_{2i}$ do not intersect.
A more general case is when the regressors $\mathbf{x}_{1i}$ consist of two components, where one component is different from the observed control function $\mathbf{x}_{2i}$ and the other is part of $\mathbf{x}_{2i}$. That is, $\mathbf{x}_{1i} = (\mathbf{x}_{11i},\mathbf{x}_{12i})$, where $\mathbf{x}_{11i}$ does not share any elements with $\mathbf{x}_{2i}$ and $\mathbf{x}_{11i}$ is nonempty, and $\mathbf{x}_{12i} \subset \mathbf{x}_{2i}$. Let $\beta^0_2 = (\beta^0_{21},\beta^0_{22}), \beta^0_3 = (\beta^0_{31},\beta^0_{32})$ conformable to the dimensions of $(\mathbf{x}_{11i},\mathbf{x}_{12i})$. Similarly let $\beta_2 = (\beta_{21},\beta_{22}), \beta_3 = (\beta_{31},\beta_{32}).$
In this case, with a properly modified rank condition of $\mathbf{z}_{(2),i}$ and $\mathbf{w}_{(2),i}$ which excludes the variables associated with $\mathbf{x}_{12,i}$ and $\sum_{j=1,\neq i}^N g_{ij} \mathbf{x}_{12,j}$, we can identify the coefficients $\beta^0_{(2)} := (\beta^0_1,\beta^0_{21},\beta^0_{31})$ using the same argument that leads to the identification in (\ref{eq.identification.alternative}).
However, we cannot identify the coefficients that correspond to the variable $\mathbf{x}_{12,i}$ and $\sum_{j=1,\neq i}^N g_{ij} \mathbf{x}_{12,j}$. The reason is that controlling the network endogeneity with the control variable $(\mathbf{x}_{2i},a_i)$ wipes out the information in $(\mathbf{x}_{12,i},\sum_{j=1,\neq i}^N g_{ij} \mathbf{x}_{12,j})$:
\begin{align*}
\mathbf{x}_{12,i} - \mathbb{E}[ \mathbf{x}_{12,i} | \mathbf{x}_{2i},a_i] &= 0 \\
\sum_{j=1,\neq i}^N g_{ij} \mathbf{x}_{12,j} - \mathbb{E}\left[\sum_{j=1,\neq i}^N g_{ij} \mathbf{x}_{12,j} \left| \mathbf{x}_{2i},a_i \right. \right] &\rightarrow_p 0,
\end{align*}
where the second convergence holds because $\sum_{j=1,\neq i}^N g_{ij} \mathbf{x}_{12,j}$ converges to a function that depends only on $(\mathbf{x}_{2i}, a_i)$ (see Section \ref{appendix: distribution of est} in the Supplementary Appendix.).
Throughout the rest of the paper, when we consider $(\mathbf{x}_{2i},a_i)$ as control function, we will without loss of generality apply the restriction in Assumption \ref{as:basic.alternative} that $\mathbf{x}_{1i}$ and $\mathbf{x}_{2i}$ do not overlap.
\section{Estimation}\label{section: estimation}
In this section we present two estimation methods. In subsections \ref{subsection: estimation, network formation} and \ref{sec: estimation with x,a as control} we discuss estimation using $a_i$ and $(\mathbf{x}_{2i},a_i)$ as control functions, respectively.
\subsection{With $a_i$ as Control Function}\label{subsection: estimation, network formation}
The identification scheme of Theorem \ref{theorem:identification} identifies the parameter of interest $\beta^0$ with the two step procedure: (i) control $a_i$ in the outcome equation and yield $y_i - \mathbb{E}(y_i|a_i) = (\mathbf{w}_i - \mathbb{E}(\mathbf{w}_i|a_i))'\beta^0 + \upsilon_i - \mathbb{E}(\upsilon_i)$, and then (ii) use $\mathbf{z}_i - \mathbb{E}(\mathbf{z}_i|a_i)$ as IVs for $\mathbf{w}_i - \mathbb{E}(\mathbf{w}_i|a_i)$. If we observe $a_i$ and know the conditional mean functions $\mathbf{h}(a_i)=(h^y(a_i),\mathbf{h}^{w}(a_i),\mathbf{h}^z(a_i))
:=(\mathbb{E}[y_i|a_i], \mathbb{E}[\mathbf{w}_i|a_i],\mathbb{E}[\mathbf{z}_i|a_i])$, then $\beta^0$ can be estimated using 2SLS as
\begin{align}\label{def.2sls.inf.ai}
&\widehat{\beta}_{2SLS}^{\text{inf}} \nonumber \\
&= \left[ \sum_{i=1}^N (\mathbf{w}_i - \mathbf{h}^{w}(a_i))(\mathbf{z}_i - \mathbf{h}^{z}(a_i))' \left(\sum_{i=1}^N (\mathbf{z}_i - \mathbf{h}^{z}(a_i)) (\mathbf{z}_i - \mathbf{h}^{z}(a_i))'\right)^{-1}
\sum_{i=1}^N (\mathbf{z}_i - \mathbf{h}^{z}(a_i)) (\mathbf{w}_i - \mathbf{h}^{w}(a_i))' \right]^{-1} \nonumber \\
&\times
\left[ \sum_{i=1}^N (\mathbf{w}_i - \mathbf{h}^{w}(a_i))(\mathbf{z}_i - \mathbf{h}^{z}(a_i))' \left(\sum_{i=1}^N (\mathbf{z}_i - \mathbf{h}^{z}(a_i)) (\mathbf{z}_i - \mathbf{h}^{z}(a_i))'\right)^{-1}
\sum_{i=1}^N (\mathbf{z}_i - \mathbf{h}^{z}(a_i)) (y_i - h^{y}(a_i))' \right].
\end{align}
However, since the individual heterogeneity $a_i$ is not observed and the conditional mean functions $\mathbf{h}(a_i) = (\mathbb{E}(y_i|a_i), \mathbb{E}(\mathbf{w}_i|a_i),\mathbb{E}(\mathbf{z}_i|a_i))$ are not known either, the estimator $\widehat{\beta}_{2SLS}^{\text{inf}}$ is not feasible.
A natural implementation of the infeasible estimator $\widehat{\beta}_{2SLS}^{\text{inf}}$ is to replace the conditional mean function $\mathbf{h}(a_i)$ with its estimate. Suppose that $\widehat{a}_i$ is an estimator of $a_i$ and $\widehat{\mathbf{h}}(\widehat{a}_i)$ is a nonparametric estimator of $\mathbf{h}(a_i)$. Then we can implement the infeasible estimator $\widehat{\beta}_{2SLS}^{\text{inf}}$ with
\begin{align}
&\widehat{\beta}_{2SLS} \\
&:=
\left[ \sum_{i=1}^N (\mathbf{w}_i - \widehat{\mathbf{h}}^w(\widehat{a}_i))(\mathbf{z}_i - \widehat{\mathbf{h}}^z(\widehat{a}_i))'
\left(\sum_{i=1}^N (\mathbf{z}_i - \widehat{\mathbf{h}}^z(\widehat{a}_i)) (\mathbf{z}_i - \widehat{\mathbf{h}}^z(\widehat{a}_i))'\right)^{-1}
\sum_{i=1}^N (\mathbf{z}_i - \widehat{\mathbf{h}}^z(\widehat{a}_i)) (\mathbf{w}_i - \widehat{\mathbf{h}}^w(\widehat{a}_i))' \right]^{-1} \nonumber \\
&\times
\left[ \sum_{i=1}^N (\mathbf{w}_i - \widehat{\mathbf{h}}^w(\widehat{a}_i))(\mathbf{z}_i - \widehat{\mathbf{h}}^z(\widehat{a}_i))'
\left(\sum_{i=1}^N (\mathbf{z}_i - \widehat{\mathbf{h}}^z(\widehat{a}_i)) (\mathbf{z}_i - \widehat{\mathbf{h}}^z(\widehat{a}_i))'\right)^{-1}
\sum_{i=1}^N (\mathbf{z}_i - \widehat{\mathbf{h}}^z(\widehat{a}_i)) (y_i - \widehat{\mathbf{h}}^y(\widehat{a}_i))' \right].
\label{def.betahat}
\end{align}
See Section \ref{appendix:beta_hat} in the Appendix for more details on the estimator $\widehat{\beta}_{2SLS}$.
{\noindent \bf Estimation of $\mathbf{h}(\cdot)$:} We can estimate $\mathbf{h}(\cdot)$ using various standard nonparametric methods. In this paper we consider a (linear) sieve estimation method.\footnote{In principle we can use other nonparametric estimation methods such as kernel smoothing or local polynomial methods.\label{kernel}}
Suppose that $h^l(a)$ is the $l^{th}$ element in $\mathbf{h}(a)$ for $l=1,...,L$, where $L$ is the dimension of $(y_i,\mathbf{w}_i',\mathbf{z}_i')'$.
The sieve estimation method assumes that each function $h^l(a)$, $l=1,...,L$ is well approximated by a linear combination of base functions $(q_1(a),...,q_{K_N}(a))$:
\begin{equation}
h^l(a) \cong \sum_{k=1}^{K_N}q_k(a)\alpha_k^l, \label{eq.approximation}
\end{equation}
as the truncation parameter $K_N \rightarrow \infty$.
A linear sieve (or series) estimator of a function, for example $\widehat{h}^y(\widehat{a}_i)$, is the OLS projection of $y_i$ on the sieve basis $\mathbf{q}^K(\cdot) = (q_1(\cdot),...,q_K(\cdot))'$ with $\widehat{a}_i$ plugged in,
\[
\widehat{h}^y(\widehat{a}_i) := \mathbf{q}^K(\widehat{a}_i)' \left( \sum_{i=1}^N \mathbf{q}^K(\widehat{a}_i)\mathbf{q}^K(\widehat{a}_i)' \right)^{-1}
\sum_{i=1}^N \mathbf{q}^K(\widehat{a}_i)y_i.
\]
For the regularity conditions of the sieve basis $\mathbf{q}^K(a_i)$, we impose standard conditions such as those proposed by \cite{Newey1997} and \citet{Li2008a}. These assumptions ensure that $\sum_{i=1}^N \mathbf{q}^K(a_i) \mathbf{q}^K(a_i)'$ is asymptotically non-singular and control the rate of approximation of the sieve estimator. These assumptions are formally stated in Assumptions \ref{assumption:sieve basis} and \ref{assumption:sieve basis.alternative} of the Appendix.
Additionally, we require that the sieve basis satisfy a Lipschitz condition, which allows us to control for the error introduced by the estimation of $a_i$ with $\widehat{a}_i$ in the estimation of $\widehat{\beta}_{2SLS}$\footnote{This issue is similar to the two step series estimation problem in \cite{Newey2009}. \label{newey2009reference}Other papers that investigated the problem of nonparametric or semiparametric analysis with generated regressors include \cite{Ahn1993}, \cite{Mammenetal2012}, \cite{HahnRidder2013}, and \cite{Escancianoetal2014}, for example.} (see Assumptions \ref{assumption: Lipschitz condition} and \ref{assumtion:sieve with x2, Lipschitz}).
As an example, define the polynomial sieve as follows.
Let $Pol(K_N)$ denote the space of polynomials on $[-1,1]$ of degree $K_N$,
\[
Pol(K_N)=\left\{
\nu_0 + \sum_{k=1}^{K_N}\nu_k a^k, \ a \in [-1,1], \nu_k \in \mathbb{R}
\right\}.
\]
For any $k$ we have
\[
\big| a_1^k - a_2^k \big|
= k | \tilde{a}^k | | a_1- a_2 | \leq M k | a_1- a_2 |
,
\]
where $\tilde{a} \in [-1,1]$ and $M$ is a finite constant.
In sieve estimations an important issue is choosing the truncation parameter $K_N$. Well-known procedures for selecting $K_N$ are Mallows' $C_P$, generalized cross-validation and leave-one-out cross-validation. For more on these methods see Chapter 15.2 in \cite{Li2008a}, \cite{li1987asymptotic}, \cite{wahba1985comparison}, \cite{andrews1991asymptotic} and \cite{hansen2014nonparametric}. However, these methods are mainly applicable when the observations are cross-sectionally independent, which is not true in our case, especially when the network is dense\label{dense3}, as we assume. Developing a data-driven choice of $K_N$ is beyond the scope of this paper and we leave it for future work.
\bigskip
{\noindent \bf Estimation of $a_i$:} A desired estimator of $a_i$ should satisfy the following high level condition.
\begin{assumption}[Estimation of $a_i$]\label{assumption: estimation of a_i}
We assume that we can estimate $a_i$ with $\widehat{a}_i$ such that $\max_i |\widehat{a}_i-a_i|=O_p\left( \zeta_{a}(N)^{-1} \right)$,
where $\zeta_{a}(N) \rightarrow \infty$ as $N \rightarrow \infty$, satisfying Assumption \ref{assumption: Lipschitz condition} in the Appendix.
\end{assumption}
{\color{black} Here $\zeta_a(N)$ is the order of magnitude that measures the Lipschitz smoothness of the sieve basis. The assumption puts restrictions on the uniform bound of the convergence rate of $\widehat{a}_i$, and we need a more accurate estimator of $a_i$ when the average curvature of the sieve basis is larger.}
For the purpose of our paper, any estimation method that yields an estimator $\widehat{a}_i$ satisfying the restriction in Assumption \ref{assumption: estimation of a_i} can be adopted.
For example, assuming the parametric specification as in (\ref{model.network.formation.parametric}),
\begin{equation}
d_{ij}=\mathbb{I}(t(\mathbf{x}_{2i},\mathbf{x}_{2j})'\lambda+a_i+a_j \geq u_{ij}) \label{model.network.formation.parametric.2}
\end{equation}
with regularity conditions of Assumption \ref{assumption:Graham} in the Appendix, {\color{black} including the error $u_{ij}$ following a logistic distribution,}
\cite{Graham2017} showed that the joint maximum likelihood estimator that solves
\begin{align*}
& (\widehat{a}_1,...,\widehat{a}_N)
\\
&:= \operatorname*{argmax}_{\lambda,(a_1,...,a_N) }
\left(\sum_{i=1}^N\sum_{j<i}d_{ij} \exp\left(t(\mathbf{x}_{2i},\mathbf{x}_{2j})'\lambda+a_i+a_j
\right)-\ln\left[1+\exp(t(\mathbf{x}_{2i},\mathbf{x}_{2j})'\lambda+a_i+a_j)\right]\right)
\end{align*}
satisfies
\begin{equation}
\sup_{1\leq i\leq N} |\widehat{a}_i-a_i|
\leq O\left(\sqrt{\frac{\ln N}{N}}\right) \label{eq.a_ihat.unform.convergence}
\end{equation}
with probability $1-O(N^{-2})$.
In this case we have $\zeta_a(N) = \sqrt{\frac{N}{\ln N}}$.
Notice that the requirement that the network formation in (\ref{model.network.formation.parametric.2}) be dense is necessary for $\widehat{a}_i$ to satisfy the desired uniform convergence rate in (\ref{eq.a_ihat.unform.convergence}). \label{remark dense-limitation}
Examples of other estimation methods include \cite{Fernandez-val}, \cite{jochmans2016modified}, \cite{dzemski2017empirical}, and \cite{jochmans2018semiparametric}.\label{dzemski, jochman}
\subsection{With $(\mathbf{x}_{2i},a_i)$ as Control Function}\label{sec: estimation with x,a as control}
As we assume in Section \ref{section: alternative identification}, we consider the case where $\mathbf{x}_{1i}$ and $\mathbf{x}_{2i}$ do not overlap. When $a_i$ is observed and the conditional expectations $\mathbf{h}_{*}(\mathbf{x}_{2i},a_i) = (h_{*}^y(\mathbf{x}_{2i},a_i),\mathbf{h}_{*}^w(\mathbf{x}_{2i},a_i),\mathbf{h}_{*}^z(\mathbf{x}_{2i},a_i)):= (\mathbb{E}(y_i|\mathbf{x}_{2i},a_i),\mathbb{E}(\mathbf{w}_i|\mathbf{x}_{2i},a_i),\mathbb{E}(\mathbf{z}_i|\mathbf{x}_{2i},a_i))$ are known, we can estimate $\beta^0$ by the 2SLS similar to $\widehat{\beta}^{\inf}_{2SLS}$ in (\ref{def.2sls.inf.ai}),
\begin{align}\label{def.2sls.inf.x2iai}
&\bar{\beta}_{2SLS}^{\text{inf}} \nonumber \\
&= \left[ \sum_{i=1}^N (\mathbf{w}_i - \mathbf{h}_{*}^{w}(\mathbf{x}_{2i},a_i))(\mathbf{z}_i - \mathbf{h}_*^{z}(\mathbf{x}_{2i},a_i))' \left(\sum_{i=1}^N (\mathbf{z}_i - \mathbf{h}_*^{z}(\mathbf{x}_{2i},a_i)) (\mathbf{z}_i - \mathbf{h}_*^{z}(\mathbf{x}_{2i},a_i))'\right)^{-1} \right. \nonumber \\
& \left. \qquad \times \sum_{i=1}^N (\mathbf{z}_i - \mathbf{h}_*^{z}(\mathbf{x}_{2i},a_i)) (\mathbf{w}_i - \mathbf{h}_*^{w}(\mathbf{x}_{2i},a_i))' \right]^{-1} \nonumber \\
&\times
\left[ \sum_{i=1}^N (\mathbf{w}_i - \mathbf{h}_*^{w}(\mathbf{x}_{2i},a_i))(\mathbf{z}_i - \mathbf{h}_*^{z}(\mathbf{x}_{2i},a_i))' \left(\sum_{i=1}^N (\mathbf{z}_i - \mathbf{h}_*^{z}(\mathbf{x}_{2i},a_i)) (\mathbf{z}_i - \mathbf{h}_*^{z}(\mathbf{x}_{2i},a_i))'\right)^{-1} \right. \nonumber \\
& \left. \qquad \times \sum_{i=1}^N (\mathbf{z}_i - \mathbf{h}_*^{z}(\mathbf{x}_{2i},a_i)) (y_i - h_*^{y}(\mathbf{x}_{2i},a_i))' \right]^{-1}.
\end{align}
When $a_i$ is unknown and $\mathbf{x}_{2i}$ is also used in the control function, under the monotonicity condition of the link formation as in (\ref{eq.monotone.link.formation}), we can implement the infeasible estimator using the average node degree without estimating $a_i$.
To be more specific, first we denote
\begin{align*}
&\mathbb{P}( d_{ij} = 1 | \mathbf{x}_{2i},a_i) =: \text{deg}(\mathbf{x}_{2i},a_i) =: \text{deg}_i.
\end{align*}
Under the monotonicity condition in (\ref{eq.monotone.link.formation}), $(\mathbf{x}_{2i},a_i)$ and $(\mathbf{x}_{2i}, \text{deg}_i)$ are one-to-one. This implies that for any $b_i \in \{ y_i,\mathbf{w}_{i},\mathbf{z}_{i} \}$,
\[
h_*^b(\mathbf{x}_{2i},a_i) = \mathbb{E}(b_i| \mathbf{x}_{2i},a_i) = \mathbb{E}(b_i| \mathbf{x}_{2i},{\rm deg}_i) =: h_{**}^b(\mathbf{x}_{2i},{\rm deg}_i).
\]
Notice that the natural estimator of ${\rm deg}_i$ is the node degree of $i$, the number of connections with node (individual) $i$ in the network scaled by the network size:
\[
\widehat{\text{deg}}_i := \frac{1}{N-1} \sum_{j=1, \neq i}^N d_{ij}.
\]
Recall that the link $d_{ij}$ is formed by
\[
d_{ij}=\mathbb{I}( g( t(\mathbf{x}_{2i},\mathbf{x}_{2j}), a_i, a_j ) - u_{ij} \geq 0).
\]
Also recall that the unobserved link-specific error terms $u_{ij}$ are assumed to be independent of all the other variables and randomly drawn. Let $\Phi(\cdot)$ be the cdf of $u_{ij}$. Also let $\pi(\mathbf{x}_2,a)$ be the joint density function of $(\mathbf{x}_{2i},a_i)$.
Then, for each $(\mathbf{x}_{2i},a_i)$, by the WLLN conditioning on $(\mathbf{x}_{2i},a_i)$, we have
\begin{align}
\widehat{\text{deg}}_i &:= \frac{1}{N-1} \sum_{j=1, \neq i}^N \mathbb{I}(g( t(\mathbf{x}_{2i},\mathbf{x}_{2j}), a_i, a_j ) - u_{ij}\geq 0) \nonumber \\
&\rightarrow_p \int \Phi\left( g( t(\mathbf{x}_{2i},\mathbf{x}_{2}), a_i, a ) \right) \pi(\mathbf{x}_2,a) d\mathbf{x}_2 da \nonumber \\
&= \mathbb{P}( d_{ij} = 1 | \mathbf{x}_{2i},a_i) \nonumber \\
&=: \text{deg}_i > 0 \label{eq.deg.limit}
\end{align}
as the network size $N$ grows to infinity. Here the limit of the average network $\rm{deg}_i > 0$ follows since we assume the network is dense.
This shows that $\widehat{\text{deg}}_i$ can be used as an estimator of $\text{deg}_i$.
In fact, we can show that under the regularity conditions in Assumption \ref{assumption:limit.dist} in the Appendix, $\sup_{i} \mathbb{E} [( \sqrt{N} (\widehat{\text{deg}}_i - \text{deg}_i ))^{2B} ] < \infty$ for any finite integer $B \geq 2$, from which we can deduce that
\begin{equation}
\max_{1 \leq i \leq N} |\widehat{\text{deg}}_i - \text{deg}_i | = O_p\left(\zeta_{deg}(N)^{-1}\right), \label{eq.estimation.deg}
\end{equation}
where
\begin{equation*}
\zeta_{deg}(N):= o(1) N^{\frac{B-1}{2B}}.
\end{equation*}
This corresponds to the regularity condition in Assumption \ref{assumption: estimation of a_i}.
Suppose that $\mathbf{r}^{K}(\mathbf{x}_{2i},\text{deg}_i) =(r_1(\mathbf{x}_{2i},\text{deg}_i),\ldots,r_{K}(\mathbf{x}_{2i},\text{deg}_i))'$ is a sieve basis of the unknown function $\mathbf{h}_*(\mathbf{x}_{2i},a_i)$.
For each $b_i \in \{ y_i,\mathbf{w}_{i},\mathbf{z}_{i} \}$, a sieve estimator of $h_{**}^b(\mathbf{x}_{2i},\text{deg}_i) = \mathbb{E}(b_i| \mathbf{x}_{2i},a_i)$ is the OLS projection of $b_i$ on $\mathbf{r}^{K}(\mathbf{x}_{2i},\widehat{\text{deg}}_i)$. For example,
\begin{align*}
\widehat{h}_*^y(\mathbf{x}_{2i},a_i)
&= \widehat{h}_{**}^y(\mathbf{x}_{2i},\text{deg}_i) \\
&= \mathbf{r}^K(\mathbf{x}_{2i},\widehat{\text{deg}}_i)'
\left( \sum_{i=1}^N \mathbf{r}^K(\mathbf{x}_{2i},\widehat{\text{deg}}_i)\mathbf{r}^K(\mathbf{x}_{2i},\widehat{\text{deg}}_i)' \right)^{-1}
\sum_{i=1}^N \mathbf{r}^K(\mathbf{x}_{2i},\widehat{\text{deg}}_i)y_i.
\end{align*}
Then, we have
\begin{align}\label{def.betabar}
&\bar{\beta}_{2SLS} \nonumber \\
&= \left[ \sum_{i=1}^N (\mathbf{w}_i - \widehat{\mathbf{h}}_{*}^{w}(\mathbf{x}_{2i},a_i))(\mathbf{z}_i - \widehat{\mathbf{h}}_*^{z}(\mathbf{x}_{2i},a_i))' \left(\sum_{i=1}^N (\mathbf{z}_i - \widehat{\mathbf{h}}_*^{z}(\mathbf{x}_{2i},a_i)) (\mathbf{z}_i - \widehat{\mathbf{h}}_*^{z}(\mathbf{x}_{2i},a_i))'\right)^{-1} \right. \nonumber \\
& \left. \qquad \times \sum_{i=1}^N (\mathbf{z}_i - \widehat{\mathbf{h}}_*^{z}(\mathbf{x}_{2i},a_i)) (\mathbf{w}_i - \widehat{\mathbf{h}}_*^{w}(\mathbf{x}_{2i},a_i))' \right]^{-1} \nonumber \\
&\times
\left[ \sum_{i=1}^N (\mathbf{w}_i - \widehat{\mathbf{h}}_*^{w}(\mathbf{x}_{2i},a_i))(\mathbf{z}_i - \widehat{\mathbf{h}}_*^{z}(\mathbf{x}_{2i},a_i))' \left(\sum_{i=1}^N (\mathbf{z}_i - \widehat{\mathbf{h}}_*^{z}(\mathbf{x}_{2i},a_i)) (\mathbf{z}_i - \widehat{\mathbf{h}}_*^{z}(\mathbf{x}_{2i},a_i))'\right)^{-1} \right. \nonumber \\
& \left. \qquad \times \sum_{i=1}^N (\mathbf{z}_i - \widehat{\mathbf{h}}_*^{z}(\mathbf{x}_{2i},a_i)) (y_i - \widehat{h}_*^{y}(\mathbf{x}_{2i},a_i))' \right]^{-1}.
\end{align}
For more details see Section \ref{appendix:beta_bar} in the Appendix.
The two different estimators $\widehat{\beta}_{2SLS}$ and $\bar{\beta}_{2SLS}$ are implemented using different control functions, and these two approaches have their own pros and cons. For $\widehat{\beta}_{2SLS}$, a good estimator of $a_i$ is required, which imposes restrictions on the network formation model (\ref{model.network.formation}) in the form of (\ref{model.network.formation.parametric}). Compared to this, the estimator $\bar{\beta}_{2SLS}$ that uses $(\mathbf{x}_{2i}, {\rm deg}_i)$ as control functions does not require a restriction like (\ref{model.network.formation.parametric}). It requires only the monotonicity of the net surplus function as in (\ref{eq.monotone.link.formation}) of Section \ref{sec: model of network formation}. However, $\bar{\beta}_{2SLS}$ has disadvantages: because it uses $x_{2i}$ as a part of the control function, as discussed in Section \ref{section: alternative identification}, this approach cannot identify and estimate the coefficients of the regressor $\mathbf{x}_{2i}$ if $\mathbf{x}_{2i}$ is a relevant regressor of the outcome.\label{remark-comparison-of-approaches}
Later in Section \ref{section: monte carlo}, where we present the Monte Carlo simulations, we compare the finite sample properties of $\widehat{\beta}_{2SLS}$ and $\bar{\beta}_{2SLS}$ in both dense and sparse network setups.
\section{Limit Distribution and Standard Error}\label{subsection: estimation, limiting distributin of estimator}
In this section we present the asymptotic distributions of the two 2SLS estimators $\widehat{\beta}_{2SLS}$ and $\bar{\beta}_{2SLS}$, and show how to estimate standard errors. We also discuss key technical issues in deriving the limits. All details of the technical derivations and proofs can be found in the Appendix.
\subsection{Limiting Distribution and Standard Error of $\widehat{\beta}_{2SLS}$}
Recall the definitions $h^{y}(a_i):= \mathbb{E}[y_i|a_i], \quad h^{\upsilon}(a_i):= \mathbb{E}[\upsilon_i|a_i], \quad \mathbf{h}^{\mathbf{w}}(a_i) := \mathbb{E} (\mathbf{w}_i|a_i), \quad \mathbf{h}^{\mathbf{z}}(a_i) := \mathbb{E} (\mathbf{z}_i|a_i).$
Define
$\eta^{y}_i: = y_i - h^{y}(a_i), \quad \eta^{\upsilon}_i: = \upsilon_i - h^{\upsilon}(a_i), \quad \eta^{\mathbf{w}}_i = \mathbf{w}_i - \mathbf{h}^{\mathbf{w}}(a_i), \quad \eta_i^{\mathbf{z}} = \mathbf{z}_i - \mathbf{h}^{\mathbf{z}}(a_i).$
Let $\bm{\eta}_N^{\upsilon} = (\eta^{\upsilon}_1,...,\eta^{\upsilon}_N)'$ and $\mathbf{H}^{\upsilon}_N(\mathbf{a}_N) = (h^{\upsilon}(a_1),...,h^{\upsilon}(a_N))'$.
Let $\widehat{h}^{\upsilon}(a_i)$, $\widehat{\mathbf{h}}^{\mathbf{w}}(a_i)$, and $\widehat{\mathbf{h}}^{\mathbf{z}}(a_i)$ denote the sieve estimators of $h^{\upsilon}(a_i)$, $h^{\mathbf{w}}(a_i)$ and $h^{\mathbf{z}}(a_i)$, respectively.
In the Appendix, we derive the asymptotic distribution of $\widehat{\beta}_{2SLS}$ in three steps.
First, we show that the sampling error caused by the use of $\hat{a}_i$ instead of $a_i$ is asymptotically negligible (see Lemma \ref{lemma: error from A_hat} of the Supplementary Appendix \ref{appenxid: error from A_hat-A}.).
Next, we control the error introduced by the non-parametric estimation of $h^{l}(a_i)$, where $l \in \{\upsilon,\mathbf{w},\mathbf{z}\}$. In Lemma \ref{lemma: series approximation error} of the Supplementary Appendix \ref{appendix: series approximation error} we show that under the regularity conditions, the estimation error in $\widehat{h}^l(a_i)$ vanishes at a suitable rate. Combining these two, we deduce
\[
\sqrt{N} (\widehat{\beta}_{2SLS} - \widehat{\beta}^{\inf}_{2SLS}) = o_p(1).
\]
The last step is to derive the limiting distribution of the infeasible estimator $\sqrt{N} ( \widehat{\beta}^{\inf}_{2SLS} - \beta^0)$.
In the Supplementary Appendix \ref{appendix: distribution of est} we show the following:
\begin{align}
\frac{1}{N}\sum_{i=1}^N(\mathbf{w}_i-\mathbf{h}^{\mathbf{w}}(a_i))
(\mathbf{z}_i-\mathbf{h}^\mathbf{z}(a_i))' &\xrightarrow{p} \mathbf{S}^{\mathbf{w}\mathbf{z}} \label{eq.WLLN.wz} \\
\frac{1}{N}\sum_{i=1}^N(\mathbf{z}_i-\mathbf{h}^{\mathbf{z}}(a_i))
(\mathbf{z}_i-\mathbf{h}^{\mathbf{z}}(a_i))' &\xrightarrow{p} \mathbf{S}^{\mathbf{z}\mathbf{z}} \label{eq.WLLN.zz} \\
\frac{1}{\sqrt{N}}\sum_{i=1}^N(\mathbf{z}_i-\mathbf{h}^{\mathbf{z}}(a_i))
\eta^{\upsilon}_i
& \Rightarrow \mathcal{N}(0,\mathbf{S}^{\mathbf{z}\mathbf{z}\sigma}), \label{eq.asy.normal.a_i}
\end{align}
where the closed forms of the limits $\mathbf{S}^{\mathbf{w}\mathbf{z}}$ and $\mathbf{S}^{\mathbf{z}\mathbf{z}}$ are found in Lemma \ref{lemma: limit of S^ZZ and S^WZ} and $\mathbf{S}^{\mathbf{z}\mathbf{z}\sigma}$ in Lemma \ref{lemma:numerator.limit variance} of Supplementary Appendix.
Notice that the derivation of the limiting distribution in (\ref{eq.asy.normal.a_i}) allows $\eta^{\upsilon}_i = \upsilon_i - \mathbb{E}(\upsilon_i|a_i)$ to be conditionally heteroskedastic, and so \label{heteroskedasticity} $\sigma^2(\mathbf{x}_i,a_i) := \mathbb{E}[(\upsilon_i - \mathbb{E}[ \upsilon_i|a_i])^2| \mathbf{x}_i,a_i]$ is allowed to depend on $(\mathbf{x}_i,a_i)$.
Combining all the limit results leads to the following theorem.
\begin{theorem}[Limiting Distribution]\label{theorem: central limit theorem}
Suppose that Assumptions \ref{as:basic}, \ref{assumption: rank}, \ref{assumption: estimation of a_i}, \ref{assumption:sieve basis}, \ref{assumption: Lipschitz condition}, and \ref{assumption:limit.dist}(i)-(v) in the Appendix hold. Then, we have
\begin{align*}
\sqrt{N}(\widehat{\beta}_{2SLS}-\beta^0)
& \Rightarrow
\mathcal{N}
\left(0, \Omega \right),
\end{align*}
where
\begin{align}
\Omega &= \left(\mathbf{S}^{\mathbf{w}\mathbf{z}}\left(\mathbf{S}^{\mathbf{z}\mathbf{z}}\right)^{-1}(\mathbf{S}^{\mathbf{w}\mathbf{z}})^{\prime}\right)^{-1}
\left(
\mathbf{S}^{\mathbf{w}\mathbf{z}}\left(\mathbf{S}^{\mathbf{z}\mathbf{z}}\right)^{-1}\mathbf{S}^{\mathbf{z}\mathbf{z}\sigma}
\left(\mathbf{S}^{\mathbf{z}\mathbf{z}}\right)^{-1} (\mathbf{S}^{\mathbf{w}\mathbf{z}})^{\prime}
\right) \left(\mathbf{S}^{\mathbf{w}\mathbf{z}}\left(\mathbf{S}^{\mathbf{z}\mathbf{z}}\right)^{-1}(\mathbf{S}^{\mathbf{w}\mathbf{z}})^{\prime}\right)^{-1}.
\end{align}
\end{theorem}
The theorem requires several regularity conditions which are presented in Appendix \ref{appendix: assumptions}.
In addition to conditions of random sampling of
$(y_i,\mathbf{x}_i,a_i)$ in Assumption \ref{as:basic} and the full rank condition in Assumption \ref{assumption: rank}, we assume conditions that ensure $a_i$ can be consistently estimated, and that the error between $\mathbf{h}(a_i)$ and $\widehat{\mathbf{h}}(\widehat{a}_i)$ converges to zero at a suitable rate (Assumptions \ref{assumption: estimation of a_i}, \ref{assumption:sieve basis} and \ref{assumption: Lipschitz condition}).
We also impose restrictions on the outcome model (\ref{model:outcome}) and the network formation model (\ref{model.network.formation}) (Assumption \ref{assumption:limit.dist}). We assume $|\beta_1^0|$ is bounded below $1$ so that the spillover effect has a unique solution, and $\| \beta_2^0 \|$ is bounded above $0$ so that the IVs are strong. We also assume the observables $(y_i,\mathbf{x}_i)$ and $\mathbf{t}_{ij}$ are bounded, and $a_i$ has a compact support in $[-1,1]$. This boundedness condition is required as a technical regularity condition that simplifies the proofs of the limits in (\ref{eq.WLLN.wz}), (\ref{eq.WLLN.zz}), and (\ref{eq.asy.normal.a_i}), which involves some uniformity in the limit. \label{boundedness1}
The asymptotic variance can be consistently estimated by
\begin{align}
\widehat{\Omega}
= \left(\widehat{\mathbf{S}}^{\mathbf{w}\mathbf{z}} \left(\widehat{\mathbf{S}}^{\mathbf{z}\mathbf{z}}\right)^{-1} (\widehat{\mathbf{S}}^{\mathbf{w}\mathbf{z}})^{\prime}\right)^{-1}
\left(
\widehat{\mathbf{S}}^{\mathbf{w}\mathbf{z}} \left(\widehat{\mathbf{S}}^{\mathbf{z}\mathbf{z}}\right)^{-1} \widehat{\mathbf{S}}^{\mathbf{z}\mathbf{z}\sigma}
\left(\widehat{\mathbf{S}}^{\mathbf{z}\mathbf{z}}\right)^{-1} (\widehat{\mathbf{S}}^{\mathbf{w}\mathbf{z}})^{\prime}
\right)
\left(\widehat{\mathbf{S}}^{\mathbf{w}\mathbf{z}}\left(\widehat{\mathbf{S}}^{\mathbf{z}\mathbf{z}}\right)^{-1}(\widehat{\mathbf{S}}^{\mathbf{w}\mathbf{z}})^{\prime} \right)^{-1},
\end{align}
where
\begin{align*}
\widehat{\mathbf{S}}^{\mathbf{w}\mathbf{z}}
&=\frac{1}{N}\sum_{i=1}^N\left(\mathbf{w}_i - \widehat{\mathbf{h}}^\mathbf{w}(\widehat{a}_i) \right)
\left(\mathbf{z}_i - \widehat{\mathbf{h}}^\mathbf{z}(\widehat{a}_i) \right)' \\
\widehat{\mathbf{S}}^{\mathbf{z}\mathbf{z}}
&=\frac{1}{N}\sum_{i=1}^N\left(\mathbf{z}_i - \widehat{\mathbf{h}}^\mathbf{z}(\widehat{a}_i) \right)
\left(\mathbf{z}_i - \widehat{\mathbf{h}}^\mathbf{z}(\widehat{a}_i) \right)' \\
\widehat{\mathbf{S}}^{ZZ\sigma^2}
&=\frac{1}{N}\sum_{i=1}^N\left(\mathbf{z}_i - \widehat{\mathbf{h}}^\mathbf{z}(\widehat{a}_i) \right)
\left(\mathbf{z}_i - \widehat{\mathbf{h}}^\mathbf{z}(\widehat{a}_i) \right)' (\widehat{\eta}^{\upsilon}_i)^2,
\end{align*}
and $\widehat{\eta}^{\upsilon}_i=y_i-\widehat{h}^y(\widehat{a}_i)-(\mathbf{w}_i-\widehat{\mathbf{h}}^{\mathbf{w}}(\widehat{a}_i))'\widehat{\beta}_{2SLS}.$
\subsection{Limiting Distribution and Standard Error of $\bar{\beta}_{2SLS}$}
The process is analogous to the one presented in the previous section.
Again, let $b_i^l$ be the $l^{th}$ element in $(y_i,\mathbf{w}_i',\mathbf{z}_i')'$.
Recall the definition that
\begin{align*}
h^l_{*}(\mathbf{x}_{2i},a_i) =\mathbb{E}[b^l_i|\mathbf{x}_{2i},a_i] = \mathbb{E}[b^l_i|\mathbf{x}_{2i},\text{deg}_i] =: h^l_{**}(\mathbf{x}_{2i},\text{deg}_i).
\end{align*} Further, let $\eta^l_{*i}=b^l_i-h^l_{*}(\mathbf{x}_{2i},a_i) = b^l - h^l_{**}(\mathbf{x}_{2i},\text{deg}_i)$, and
let $\widehat{h}^l_{**}(\mathbf{x}_{2i},\text{deg}_i)$ denote a sieve estimator of $h^l_{**}(\mathbf{x}_{2i},\text{deg}_i)$.
As in the previous section, we derive the asymptotic distribution of $\bar{\beta}_{2SLS}$ in three steps.
First, we show that the error that stems from the use of the estimate $\widehat{\text{deg}_i}$ for $\text{deg}_i$, $\widehat{h}^l_{**}(\mathbf{x}_{2i},\widehat{\text{deg}}_i) - \widehat{h}^l_{**}(\mathbf{x}_{2i},\text{deg}_i)$, is asymptotically negligible.
In the second step, we control the error introduced by the non-parametric estimation of $h^l_{**}(\mathbf{x}_{2i},\text{deg}_i)$, $\widehat{h}^l_{**}(\mathbf{x}_{2i},\text{deg}_i)-h^l_{**}(\mathbf{x}_{2i},\text{deg}_i)$. This implies
\[
\sqrt{N} (\bar{\beta}_{2SLS} - \bar{\beta}^{\inf}_{2SLS}) = o_p(1).
\]
The last step is to derive the limiting distribution of the infeasible estimator $\sqrt{N} ( \bar{\beta}^{\inf}_{2SLS} - \beta^0)$ by showing
\begin{align*}
\frac{1}{N}\sum_{i=1}^N(\mathbf{w}_i-\mathbf{h}_*^{\mathbf{w}}(\mathbf{x}_{2i},a_i))
(\mathbf{z}_i-\mathbf{h}_*^\mathbf{z}(\mathbf{x}_{2i},a_i))' &\xrightarrow{p} \mathbf{\bar{S}}^{\mathbf{w}\mathbf{z}} \\
\frac{1}{N}\sum_{i=1}^N(\mathbf{z}_i-\mathbf{h}_*^{\mathbf{z}}(\mathbf{x}_{2i},a_i))
(\mathbf{z}_i-\mathbf{h}_*^{\mathbf{z}}(\mathbf{x}_{2i},a_i))' &\xrightarrow{p} \mathbf{\bar{S}}^{\mathbf{z}\mathbf{z}} \\
\frac{1}{\sqrt{N}}\sum_{i=1}^N(\mathbf{z}_i-\mathbf{h}_*^{\mathbf{z}}(\mathbf{x}_{2i},a_i))
\eta^{\upsilon}_{*i}
& \Rightarrow \mathcal{N}(0,\mathbf{\bar{S}}^{\mathbf{z}\mathbf{z}\sigma}),
\end{align*}
Combining all the limit results we have the following theorem.
\begin{theorem}[Limiting Distribution]\label{theorem: central limit theorem beta_bar}
Suppose that Assumptions \ref{as:basic}, \ref{as:basic.alternative}, \ref{assumption:rank.alternative},
\ref{assumption:sieve basis.alternative}, \ref{assumtion:sieve with x2, Lipschitz}, and \ref{assumption:limit.dist} hold. Then, we have
\begin{align*}
\sqrt{N}(\bar{\beta}_{2SLS}-\beta^0)
& \Rightarrow
\mathcal{N}
\left(0, \bar{\Omega} \right),
\end{align*}
where
\begin{align}
\bar{\Omega} &= \left(\mathbf{\bar{S}}^{\mathbf{w}\mathbf{z}}\left(\mathbf{\bar{S}}^{\mathbf{z}\mathbf{z}}\right)^{-1}(\mathbf{\bar{S}}^{\mathbf{w}\mathbf{z}})^{\prime}\right)^{-1}
\left(
\mathbf{\bar{S}}^{\mathbf{w}\mathbf{z}}\left(\mathbf{\bar{S}}^{\mathbf{z}\mathbf{z}}\right)^{-1}\mathbf{\bar{S}}^{\mathbf{z}\mathbf{z}\sigma}
\left(\mathbf{\bar{S}}^{\mathbf{z}\mathbf{z}}\right)^{-1} (\mathbf{\bar{S}}^{\mathbf{w}\mathbf{z}})^{\prime}
\right)
\left(\mathbf{\bar{S}}^{\mathbf{w}\mathbf{z}}\left(\mathbf{\bar{S}}^{\mathbf{z}\mathbf{z}}\right)^{-1} (\mathbf{\bar{S}}^{\mathbf{w}\mathbf{z}})^{\prime}\right)^{-1}. \nonumber
\end{align}
\end{theorem}
The asymptotic result in Theorem \ref{theorem: central limit theorem beta_bar} requires the following regularity conditions which are formally presented in the Appendix. First, Assumption \ref{as:basic.alternative} assumes that the regressors in the outcome equation, $\mathbf{x}_{1i}$ and the observables in the network formation $\mathbf{x}_{2i}$ do not overlap. Assumption \ref{assumption:rank.alternative} is a full rank condition for $\bar{\beta}_{2SLS}$.
Assumptions \ref{assumption:sieve basis.alternative} and \ref{assumtion:sieve with x2, Lipschitz} regard the sieve used in constructing the estimator $\bar{\beta}_{2SLS}$.
Comparing with the assumptions assumed in Theorem \ref{theorem: central limit theorem}, Theorem \ref{theorem: central limit theorem beta_bar} does not require the high level condition of Assumption \ref{assumption: estimation of a_i} because we do not use an estimator of $a_i$. Instead it requires an additional restriction that the net surplus function in the link formation be strictly monotonic in $a_i$ conditional on $(\mathbf{x}_{2i},\mathbf{x}_{2j},a_j)$, which implies the required monotonicity condition in (\ref{eq.monotone.link.formation}).
Like in the case of $\widehat{\beta}_{2SLS}$, we allow $\eta^{\upsilon}_{*i} = \upsilon_i - \mathbb{E}(\upsilon_i|\mathbf{x}_{2i},a_i)$ to be conditionally heteroskedastic, and $\sigma^2_{*}(\mathbf{x}_i,a_i) := \mathbb{E}[(\upsilon_i - \mathbb{E}[ \upsilon_i|\mathbf{x}_{2i},a_i])^2| \mathbf{x}_i,a_i]$ is allowed to depend on $(\mathbf{x}_i,a_i)$.
The asymptotic variance can be consistently estimated by
\begin{align}
\widehat{\bar{\Omega}} &= \left(\widehat{\mathbf{\bar{S}}}^{\mathbf{w}\mathbf{z}}\left(\widehat{\mathbf{\bar{S}}}^{\mathbf{z}\mathbf{z}}\right)^{-1} (\widehat{\mathbf{\bar{S}}}^{\mathbf{w}\mathbf{z}})^{\prime})\right)^{-1}
\left(
\widehat{\mathbf{\bar{S}}}^{\mathbf{w}\mathbf{z}}\left(\widehat{\mathbf{\bar{S}}}^{\mathbf{z}\mathbf{z}}\right)^{-1} \widehat{\mathbf{\bar{S}}}^{\mathbf{z}\mathbf{z}\sigma}
\left(\widehat{\mathbf{\bar{S}}}^{\mathbf{z}\mathbf{z}}\right)^{-1} (\widehat{\mathbf{\bar{S}}}^{\mathbf{w}\mathbf{z}})^{\prime})
\right)
\left(\widehat{\mathbf{\bar{S}}}^{\mathbf{w}\mathbf{z}}\left(\widehat{\mathbf{\bar{S}}}^{\mathbf{z}\mathbf{z}}\right)^{-1}(\widehat{\mathbf{\bar{S}}}^{\mathbf{w}\mathbf{z}})^{\prime})\right)^{-1},
\end{align}
where
\begin{align*}
\widehat{\mathbf{\bar{S}}}^{\mathbf{w}\mathbf{z}}
&=\frac{1}{N}\sum_{i=1}^N\left(\mathbf{w}_i - \widehat{\mathbf{h}}_{**}^\mathbf{w}(\mathbf{x}_{2i},\widehat{\text{deg}}_i) \right)
\left(\mathbf{z}_i - \widehat{\mathbf{h}}_{**}^\mathbf{z}(\mathbf{x}_{2i},\widehat{\text{deg}}_i) \right)' \\
\widehat{\mathbf{\bar{S}}}^{\mathbf{z}\mathbf{z}}
&=\frac{1}{N}\sum_{i=1}^N\left(\mathbf{z}_i - \widehat{\mathbf{h}}_{**}^\mathbf{z}(\mathbf{x}_{2i},\widehat{\text{deg}}_i) \right)
\left(\mathbf{z}_i - \widehat{\mathbf{h}}_{**}^\mathbf{z}(\mathbf{x}_{2i},\widehat{\text{deg}}_i) \right)' \\
\widehat{\mathbf{\bar{S}}}^{\mathbf{z}\mathbf{z}\sigma^2}
&=\frac{1}{N}\sum_{i=1}^N\left(\mathbf{z}_i - \widehat{\mathbf{h}}_{**}^\mathbf{z}(\mathbf{x}_{2i},\widehat{\text{deg}}_i) \right)
\left(\mathbf{z}_i - \widehat{\mathbf{h}}_{**}^\mathbf{z}(\mathbf{x}_{2i},\widehat{\text{deg}}_i) \right)' (\widehat{\eta}^{\upsilon}_{**i})^2,
\end{align*}
and $\widehat{\eta}_{**i}^{\upsilon}=y_i-\widehat{h}_{**}^y(\mathbf{x}_{2i},\widehat{\text{deg}}_i)-(\mathbf{w}_i-\widehat{\mathbf{h}}_{**}^{\mathbf{w}}(\mathbf{x}_{2i},\widehat{\text{deg}}_i))'\bar{\beta}_{2SLS}.$
\section{Monte Carlo}\label{section: monte carlo}
We consider both dense and sparse network Monte Carlo designs. In the dense network case links are formed according to\footnote{This follows the approach of \cite{Graham2017}.}
\[
d_{ij} = \mathbb{I}\left\{ x_{2i}x_{2j}\lambda_d + a_i + a_j -u_{ij} \geq 0 \right\},
\]
where $x_{2i}\in\{-1,1\}$, $\lambda_d=1$ and $u_{ij}$ follows a logistic distribution. This link rule implies that agents have a strong taste for homophilic matching since $x_{2i}x_{2j}\lambda_d=1$ when $x_{2i}=x_{2j}$ and $x_{2i}x_{2j}\lambda_d=-1$ when $x_{2i}\neq x_{2j}$.
In the sparse network case links are formed according to
\[
d_{ij} = \mathbb{I}\left\{(|x_{2i}-x_{2j}|+3)\lambda_s + a_i + a_j -u_{ij} \geq 0 \right\},
\]
with $\lambda_s=-1$. This rule also implies homophily on observable characteristics.
Individual-level degree heterogeneity is generated according to
\[
a_i=\varphi(\alpha_L\mathbb{I}\left\{ x_{2i}=-1 \right\} +\alpha_H\mathbb{I} \left\{ x_{2i}=1 \right\} + \xi_i),
\]
with $\alpha_L \leq \alpha_H$ and $\xi_i$ a centered Beta random variable $
\xi_i|x_{2i}\sim \left\{Beta(\mu_0,\mu_1)-\frac{\mu_0}{\mu_0+\mu_1}\right\}$
so that $a_i\in \left[\alpha_L-\frac{\mu_0}{\mu_0+\mu_1},\alpha_H+\frac{\mu_1}{\mu_0+\mu_1}\right]$. We choose values of the network formation parameters so that $a_i \in[-1,1]$.
In the main text we present results based on the following parameter values. In the dense network case we set
$\mu_0=1/4$, $\mu_1=3/4$, $\alpha_L=\alpha_H=-3/4$, which yields an average node degree $=23$ when $N=100$. The sparse network formation design is generated by setting $\mu_0=1$, $\mu_1=1$, $\alpha_L=\alpha_H=-1/4$, which gives an average degree $=1.78$ when $N=100$.\footnote{Results for 14 other network formation designs can be found in Section \ref{appendix: supplementary monte carlo} of the online appendix. Most results are similar to the ones presented in the main text.}
Individual outcomes are generated according to
\[
y_i=\beta_1\sum_{j=1 \atop j\neq i}^N g_{ij}y_j+\beta_2 x_{1i}+\beta_3
\sum_{j=1 \atop j\neq i}^N g_{ij}x_{1j}+h(a_i)+\varepsilon_i.
\]
In the simulations, we set $\beta_1=0.8$, $\beta_2=\beta_3=5$, $x_{1i}=3q_1+\cos(q_2)/0.8+\epsilon_i$, where $q_1,q_2\sim\mathcal{N}(x_{2i},1)$, and $\varepsilon_i,\epsilon_i\sim\mathcal{N}(0,1)$.
For $h(a_i)$ we use the following functional forms: $h(a_i)=\exp (3 a_i)$, $h(a_i)=\cos(3 a_i)$, $h(a_i)=\sin(3 a_i)$. A plot of $h(a_i)$ for these functional forms is presented in Figure \ref{figure: h}. We can see that the exponential function yields a strongly increasing impact on the individual outcome, and with the cosine functions the returns are increasing up to a certain point and then decreasing; however the sine function gives a more irregular pattern.
We estimate the outcome equation coefficients $(\beta_1,\beta_2, \beta_3)$ using the standard 2SLS estimator for peer effects and the Hermite polynomial sieve as well as a polynomial sieve.
For the dense network case, we estimate $a_i$ using $\widehat{a}_i$ and implement the following control functions:
using a control function linear in $\widehat{a}_i$,
$\widehat{h}(\widehat{a}_i)$, $\widehat{h}(a_i)$,
$\widehat{h}(\widehat{\rm deg}_i,x_{2i})$\footnote{
{Note that since $x_{2i}$ is discrete with a finite support, $\{ x_1,...,x_M \}$,} we have
$
r(x_{2i},{\rm deg_i}) = \sum_{m=1}^M r(x_m,{\rm deg_i}) \mathbb{I}\{ x_{2i} = x_m \}.
$
We can then approximate
$
r(x_{2i}, {\rm deg_i}) \simeq \sum_{k=1}^{K_N} \left\{ \sum_{m=1}^M \alpha_{m,k} q_k^d(\rm deg_i) \mathbb{I}\{ x_{2i} = x_m \} \right\}.
$}, and
$h(a_i)$. For the sparse network case the estimator of $a_i$ is not reliable\footnote{To estimate $a_i$, we use the JMLE proposed in \cite{Graham2017}. As \cite{Graham2017} states, in sparse designs the JMLE rarely even exists, rendering it unusable in practice
when the network is too sparse. See \cite{Graham2017} for more details.} and we implement the following control functions: linear in $a_i$, $\widehat{h}(a_i)$, $\widehat{h}(\widehat{\rm deg}_i,x_{2i})$ and $h(a_i)$. In both the dense and sparse setup we also implement a benchmark model with no control for the endogeneity of the network.
In the paper, due to space limitations, we present Monte Carlo results obtained using the Hermite polynomial sieve with $K_N=4$.
Specifically, Tables \ref{table: MC dense main text} and \ref{MC sparse main text} include results for the dense and sparse network specifications, respectively. Results for the other orders of $K_N$ are not notably different; in the Online Supplement we provide results for fourteen other network formation designs, for $K_N=4,8$ and for the Hermite polynomial and polynomial sieve functions.
\begin{figure}[!h]\caption{\bf $h(a_i)$ for selected functional forms of $h(a_i)$}\label{figure: h}
\includegraphics[scale=.4]{ha.png}
\end{figure}
We also perform conventional leave-one-out cross validation to find data-dependent $K_N$ (chosen as the $K_N$ that minimizes the Root Mean Square Error (RMSE) of the prediction based on the leave-one-out estimator, see for example \cite{andrews1991asymptotic}, \cite{hansen2014nonparametric}). We report the statistics on the cross-validation in Table \ref{table: CV}.
The differences in RMSE are very small between the different values of $K_N$.\\
\label{remark: MC discussion}
Analyzing the Monte Carlo results for the dense network specification in Table \ref{table: MC dense main text}, we can see that, as expected from our asymptotic theories, the control functions $\widehat{h}(\widehat{a}_i)$ and $\widehat{h}(\widehat{deg}_i,x_{2i})$ perform better than the estimator with a linear control function, as well as the estimator that does not control for the endogeneity of the network in terms of mean bias. This difference is more pronounced in the case when $h(a_i)$ is the sine or cosine function.
Both the control for degree approach and the control function that uses $\widehat{h}(\widehat{a}_i)$ yield a low bias and have the correct size on all coefficients in all cases. In the simulations we also implemented the control function $\widehat{h}(a_i)$, that is, using the true $a_i$ instead of $\widehat{a}_i$. These results are very similar to the ones obtained using $\widehat{h}(\widehat{a}_i)$, which is in line with the estimator $\widehat{a}_i$ having a very low bias, as detailed in the table footnotes. This suggest that the approach of using $\widehat{h}(\widehat{a}_i)$ as a control function works very well when a highly precise estimator of $a_i$ is available (for example when the network size $N$ is large.).
Looking at Table \ref{MC sparse main text} and the results for the sparse design, we can see that the control for degree approach performs very well across all functional forms of $h(a_i)$. In the sparse setup, the bias of all estimates, including those that do not control for the endogeneity of the network, is small. However, the size of the no control and linear control estimates is not correct. If a precise estimator of $a_i$ is available, the control function $\widehat{h}(a_i)$ also performs well with low bias and correct size in all cases.
Table \ref{table: CV} shows that the performance of the estimators does not differ notably for different values of $K_N$. As for the choice of $K_N$ we present in the tables, we have run simulations for a range of values of $K_N$ and the results did not differ significantly. As deriving a theory for a data driven choice of $K_N$ is beyond the scope of this paper, for applied researchers we suggest estimating the model over a range of $K_N$ and seeing whether the results vary significantly. As shown in our Monte Carlo simulations, the control function approach yields results robust to the choice of $K_N$ for different non-linear functions.
\begin{table}[!h]\caption{\footnotesize {\bf Design 4 dense network: Parameter values across 1000 Monte Carlo replications with $K_N=4$ and Hermite polynomial sieve}} \label{table: MC dense main text}
\begin{threeparttable}
\centering \footnotesize
\scalebox{.8}{\begin{tabular}{cccccccccccccc}\toprule
\multicolumn{14}{c}{$h(a_i) = \exp(a_i)$}\\
\cellcolor{yellow}$N$&\multicolumn{6}{|c|}{\cellcolor{yellow}$100$}&\multicolumn{6}{|c|}{\cellcolor{yellow}$250$}&\\\hline
CF&$(0)$&$(1)$&$(2)$&$(3)$&$(4)$&$(5)$& $(0)$ &$(1)$&$(2)$&$(3)$&$(4)$&$(5)$&\\\hline
\multirow{4}{*}{$\beta_1=0.8$}& 0.002 & 0.004 &-0.000 &-0.000 &0.000 &-0.000& 0.004 &0.007 &-0.001& -0.001 & -0.001& -0.000 & \textit{mean bias} \\
&(0.010 )&(0.013 )&(0.015 )&(0.015 )&(0.024 )&(0.010 )&(0.009 )&(0.013 )&(0.015 )&(0.015 )&(0.025 )&(0.009 )&\textit{std}\\
& 0.133 & 0.115 &0.056 &0.061 &0.058 &0.058& 0.306 &0.225 &0.057& 0.057 &0.064& 0.050 &\textit{size} \\ \midrule
\multirow{4}{*}{$\beta_2=5$}& -0.003 & -0.004 &-0.000 &-0.000 &0.000 &-0.000& -0.002 &-0.004 &0.000& 0.000 &-0.000& 0.000 &\textit{mean bias} \\
&(0.031 )&(0.032 )&(0.034 )&(0.033 )&(0.035 )&(0.031 )&(0.020 )&(0.021 )&(0.020 )&(0.020 )&(0.021 )&(0.020 )&\textit{std}\\
& 0.058 & 0.069 &0.074 &0.068 &0.074 &0.057& 0.069 &0.079& 0.055 &0.059 & 0.058 &0.061 &\textit{size} \\\midrule
\multirow{4}{*}{$\beta_3=5$}& -0.032 & -0.048 &0.006& 0.008 &0.006 &0.006 &-0.066 &-0.107 &0.009& 0.013 &0.012& 0.009 &\textit{mean bias} \\
&(0.178 )&(0.217 )&(0.251 )&(0.250 )&(0.269 )&(0.174 )&(0.163 )&(0.219 )&(0.249 )&(0.248 )&(0.270 )&(0.152 )&\textit{std}\\
& 0.078 & 0.078 &0.055 &0.060 &0.061 &0.061& 0.156& 0.172 &0.051 &0.054 & 0.062 &0.050 &\textit{size} \\\midrule
\multicolumn{14}{c}{$h(a_i) = \sin(a_i)$}\\
\cellcolor{yellow}$N$&\multicolumn{6}{|c|}{\cellcolor{yellow}$100$}&\multicolumn{6}{|c|}{\cellcolor{yellow}$250$}&\\\hline
CF&$(0)$&$(1)$&$(2)$&$(3)$&$(4)$&$(5)$& $(0)$ &$(1)$&$(2)$&$(3)$&$(4)$&$(5)$&\\\hline
\multirow{4}{*}{$\beta_1=0.8$}& -0.008 & -0.005 &-0.000 &-0.000 &-0.000 &-0.000& -0.015 &-0.010 &-0.001& -0.001 & -0.001& -0.001 & \textit{mean bias} \\
&(0.014 )&(0.014 )&(0.016 )&(0.015 )&(0.025 )&(0.011 )&(0.017 )&(0.015 )&(0.015 )&(0.015 )&(0.026 )&(0.010 )&\textit{std}\\
& 0.464 & 0.160 &0.058 &0.061 &0.059 &0.045& 0.753 &0.293 &0.054& 0.057 &0.071& 0.053 &\textit{size} \\ \midrule
\multirow{4}{*}{$\beta_2=5$}& 0.007 & 0.005 &-0.001 &-0.000 &0.000 &-0.000& 0.007 &0.005 &-0.000& 0.000 &-0.000& 0.000 &\textit{mean bias} \\
&(0.033 )&(0.034 )&(0.035 )&(0.033 )&(0.036 )&(0.031 )&(0.022 )&(0.022 )&(0.021 )&(0.020 )&(0.021 )&(0.020 )&\textit{std}\\
& 0.075 & 0.072 &0.067 &0.068 &0.072 &0.060& 0.076 &0.071& 0.056 &0.059 & 0.060 &0.055 &\textit{size} \\\midrule
\multirow{4}{*}{$\beta_3=5$}& 0.113 & 0.078 &0.009& 0.008 &0.009 &0.005 &0.236 &0.165 &0.010& 0.013 &0.012& 0.012 &\textit{mean bias} \\
&(0.222 )&(0.231 )&(0.258 )&(0.250 )&(0.277 )&(0.191 )&(0.268 )&(0.249 )&(0.255 )&(0.248 )&(0.276 )&(0.177 )&\textit{std}\\
& 0.237 & 0.100 &0.057 &0.060 &0.053 &0.053& 0.646& 0.248 &0.056 &0.054 & 0.055 &0.048 &\textit{size} \\\midrule
\multicolumn{14}{c}{$h(a_i) = \cos(a_i)$}\\
\cellcolor{yellow}$N$&\multicolumn{6}{|c|}{\cellcolor{yellow}$100$}&\multicolumn{6}{|c|}{\cellcolor{yellow}$250$}&\\\hline
CF&$(0)$&$(1)$&$(2)$&$(3)$&$(4)$&$(5)$& $(0)$ &$(1)$&$(2)$&$(3)$&$(4)$&$(5)$&\\\hline
\multirow{4}{*}{$\beta_1=0.8$}& -0.009 & 0.004 &-0.000 &-0.000 &0.000 &-0.000& -0.017 &0.010 &-0.000& -0.001 & -0.001& -0.001 & \textit{mean bias} \\
&(0.016 )&(0.014 )&(0.017 )&(0.015 )&(0.025 )&(0.010 )&(0.018 )&(0.016 )&(0.015 )&(0.015 )&(0.026 )&(0.009 )&\textit{std}\\
& 0.459 & 0.104 &0.055 &0.061 &0.057 &0.053& 0.745 &0.318 &0.059& 0.057 &0.059& 0.046 &\textit{size} \\ \midrule
\multirow{4}{*}{$\beta_2=5$}& 0.009 & -0.004 &-0.000 &-0.000 &0.001 &-0.000& 0.008 &-0.005 &0.000& 0.000 &0.000& 0.000 &\textit{mean bias} \\
&(0.040 )&(0.034 )&(0.036 )&(0.033 )&(0.037 )&(0.031 )&(0.026 )&(0.022 )&(0.021 )&(0.020 )&(0.021 )&(0.020 )&\textit{std}\\
& 0.075 & 0.061 &0.062 &0.068 &0.070 &0.060& 0.084 &0.077& 0.053 &0.059 & 0.055 &0.062 &\textit{size} \\\midrule
\multirow{4}{*}{$\beta_3=5$}& 0.123 & -0.051 &0.004& 0.008 &0.004 &0.004 &0.264 &-0.161 &0.008& 0.013 &0.010& 0.011 &\textit{mean bias} \\
&(0.257 )&(0.232 )&(0.266 )&(0.250 )&(0.286 )&(0.176 )&(0.292 )&(0.258 )&(0.256 )&(0.248 )&(0.276 )&(0.157 )&\textit{std}\\
& 0.224 & 0.074 &0.053 &0.059 &0.055 &0.056& 0.640& 0.256 &0.055 &0.054 & 0.057 &0.047 &\textit{size} \\\midrule
\end{tabular}}
\begin{tablenotes}\tiny
\item CF - control function. $(0)$ - none, $(1)$ - $\lambda_a\hat{a}_i$, $(2)$ - $\hat{h}(\hat{a}_i)$, $(3)$ - $\hat{h}(a_i)$, $(4)$ - $\hat{h}(\widehat{deg}_i,x_{2i})$, $(5)$ - $h(a_i)$.
\item The network design parameters are $\mu_0=0.25$, $\mu_1=0.75$, $\alpha_L=-0.75$, $\alpha_H=-0.75$
\item Average number of links for $N=100$ is $23.0$, for $N=250$ it is $57.8$.
\item Average skewness for $N=100$ is $0.66$, for $N=250$ it is $0.89$.
\item Size is the empirical size of t-test against the truth.
\item N$=100$, $corr(a_i,\bm{x}_{2i})=0.004$,N$=250$, $corr(a_i,\bm{x}_{2i})=0.001$
\item The bias of $\hat{a}_i$ is calculated as $a_i-\hat{a}_i$.
\item For $N=100$, $\hat{a}_i$ mean bias$=0.018$, median bias$=0.008$, std$=0.271$.
\item For $N=250$, $\hat{a}_i$ mean bias$=0.007$, median bias$=0.004$, std$=0.167$.
\end{tablenotes}
\end{threeparttable}
\end{table}
\begin{table}[!h]\caption{\footnotesize {\bf Design 4 sparse network: Parameter values across 1000 Monte Carlo replications with $K_N=4$ and Hermite polynomial sieve}} \label{MC sparse main text}
\begin{threeparttable}
\centering \footnotesize
\scalebox{.8}{\begin{tabular}{cccccccccccc}\toprule
\multicolumn{12}{c}{$h(a_i) = \exp(a_i)$}\\
\cellcolor{yellow}$N$&\multicolumn{5}{|c|}{\cellcolor{yellow}$100$}&\multicolumn{5}{|c|}{\cellcolor{yellow}$250$}&\\\hline
CF&$(0)$&$(1)$&$(2)$&$(3)$&$(4)$& $(0)$ &$(1)$&$(2)$&$(3)$&$(4)$&\\\hline
\multirow{4}{*}{$\beta_1=0.8$}& 0.001 & 0.001 &0.000 &0.000 &0.000 &0.002& 0.003 &0.000 &-0.000& 0.000 &\textit{mean bias} \\
&(0.003 )&(0.003 )&(0.002 )&(0.003 )&(0.002 )&(0.004 )&(0.004 )&(0.002 )&(0.003 )&(0.002 )&\textit{std}\\
& 0.089 & 0.090 &0.052 &0.056 &0.049 &0.269& 0.257 &0.072 &0.055& 0.064 &\textit{size} \\ \midrule
\multirow{4}{*}{$\beta_2=5$}& -0.001 & -0.002 &-0.003 &-0.002 &-0.003 &-0.007& -0.008 &0.000 &0.001& 0.001 &\textit{mean bias} \\
&(0.039 )&(0.039 )&(0.033 )&(0.041 )&(0.032 )&(0.027 )&(0.027 )&(0.021 )&(0.025 )&(0.021 )&\textit{std}\\
& 0.043 & 0.046 &0.065 &0.061 &0.060 &0.078& 0.084 &0.055& 0.066 &0.049 &\textit{size} \\\midrule
\multirow{4}{*}{$\beta_3=5$}& -0.004 & -0.004 &-0.002& 0.002 &-0.002 &-0.027 &-0.028 &-0.001 &-0.000& -0.001 &\textit{mean bias} \\
&(0.076 )&(0.077 )&(0.066 )&(0.075 )&(0.065 )&(0.063 )&(0.064 )&(0.052 )&(0.058 )&(0.051 )&\textit{std}\\
& 0.034 & 0.038 &0.063 &0.063 &0.047 &0.085& 0.090& 0.056 &0.068 &0.060 &\textit{size} \\\midrule
\multicolumn{12}{c}{$h(a_i) = \sin(a_i)$}\\
\cellcolor{yellow}$N$&\multicolumn{5}{|c|}{\cellcolor{yellow}$100$}&\multicolumn{5}{|c|}{\cellcolor{yellow}$250$}&\\\hline
CF&$(0)$&$(1)$&$(2)$&$(3)$&$(4)$& $(0)$ &$(1)$&$(2)$&$(3)$&$(4)$&\\\hline
\multirow{4}{*}{$\beta_1=0.8$}& -0.000 & 0.000 &0.000 &0.000 &0.000 &-0.002& -0.000 &0.000 &-0.000& 0.000 &\textit{mean bias} \\
&(0.003 )&(0.002 )&(0.002 )&(0.003 )&(0.002 )&(0.003 )&(0.002 )&(0.002 )&(0.003 )&(0.002 )&\textit{std}\\
& 0.059 & 0.048 &0.052 &0.057 &0.051 &0.170& 0.068 &0.072 &0.059& 0.071 &\textit{size} \\ \midrule
\multirow{4}{*}{$\beta_2=5$}& -0.007 & -0.002 &-0.003 &-0.002 &-0.003 &0.005& 0.001 &0.000 &0.001& 0.000 &\textit{mean bias} \\
&(0.039 )&(0.032 )&(0.033 )&(0.041 )&(0.032 )&(0.026 )&(0.022 )&(0.021 )&(0.025 )&(0.021 )&\textit{std}\\
& 0.052 & 0.061 &0.066 &0.062 &0.059 &0.083& 0.061 &0.055& 0.073 &0.048 &\textit{size} \\\midrule
\multirow{4}{*}{$\beta_3=5$}& -0.001 & -0.001 &-0.002& 0.002 &-0.002 &0.016 &-0.001 &-0.001 &-0.000& -0.001 &\textit{mean bias} \\
&(0.078 )&(0.067 )&(0.066 )&(0.076 )&(0.065 )&(0.064 )&(0.052 )&(0.052 )&(0.058 )&(0.051 )&\textit{std}\\
& 0.059 & 0.053 &0.063 &0.067 &0.049 &0.079& 0.057& 0.056 &0.065 &0.057 &\textit{size} \\\midrule
\multicolumn{12}{c}{$h(a_i) = \cos(a_i)$}\\
\cellcolor{yellow}$N$&\multicolumn{5}{|c|}{\cellcolor{yellow}$100$}&\multicolumn{5}{|c|}{\cellcolor{yellow}$250$}&\\\hline
CF&$(0)$&$(1)$&$(2)$&$(3)$&$(4)$& $(0)$ &$(1)$&$(2)$&$(3)$&$(4)$&\\\hline
\multirow{4}{*}{$\beta_1=0.8$}& 0.001 & 0.001 &0.000 &0.000 &0.000 &0.002& 0.002 &0.000 &-0.000& 0.000 &\textit{mean bias} \\
&(0.003 )&(0.003 )&(0.002 )&(0.003 )&(0.002 )&(0.003 )&(0.003 )&(0.002 )&(0.003 )&(0.002 )&\textit{std}\\
& 0.073 & 0.081 &0.052 &0.053 &0.049 &0.197& 0.216 &0.072 &0.067& 0.068 &\textit{size} \\ \midrule
\multirow{4}{*}{$\beta_2=5$}& -0.002 & -0.002 &-0.003 &-0.002 &-0.003 &-0.005& -0.006 &0.000 &0.001& 0.000 &\textit{mean bias} \\
&(0.038 )&(0.038 )&(0.033 )&(0.041 )&(0.032 )&(0.025 )&(0.025 )&(0.021 )&(0.025 )&(0.021 )&\textit{std}\\
& 0.047 & 0.051 &0.066 &0.061 &0.062 &0.062& 0.074 &0.055& 0.065 &0.047 &\textit{size} \\\midrule
\multirow{4}{*}{$\beta_3=5$}& -0.003 & -0.003 &-0.002& 0.002 &-0.002 &-0.020 &-0.022 &-0.001 &0.000& -0.001 &\textit{mean bias} \\
&(0.073 )&(0.073 )&(0.066 )&(0.074 )&(0.065 )&(0.061 )&(0.062 )&(0.052 )&(0.059 )&(0.051 )&\textit{std}\\
& 0.038 & 0.036 &0.063 &0.065 &0.049 &0.069& 0.079& 0.056 &0.070 &0.062 &\textit{size} \\\midrule
\end{tabular}}
\begin{tablenotes}\tiny
\item CF - control function. $(0)$ - none, $(1)$ - $\lambda_a a_i$, $(2)$ - $\hat{h}(a_i)$, $(3)$ - $\hat{h}(\widehat{deg}_i,x_{2i})$, $(4)$ - $h(a_i)$.
\item The network design parameters are $\mu_0=1.00$, $\mu_1=1.00$, $\alpha_L=-0.25$, $\alpha_H=-0.25$
\item Average number of links for $N=100$ is $1.8$, for $N=250$ it is $4.5$.
\item Average skewness for $N=100$ is $0.81$, for $N=250$ it is $0.62$.
\item Size is the empirical size of t-test against the truth.
\item N$=100$, $corr(a_i,{\bf{x}}_{2i})=-0.001$,N$=250$, $corr(a_i,{\bf{x}}_{2i})=-0.002$
\end{tablenotes}
\end{threeparttable}
\end{table}
\begin{table}\caption{\footnotesize {\bf Cross-Validation results}: Parameter values across 1000 Monte Carlo replications for dense network design 4 and Hermite polynomial sieve} \label{table: CV}
\begin{threeparttable}
\centering \footnotesize
\scalebox{.9}{\begin{tabular}{cccccccccccccc}\toprule
&\cellcolor{yellow}$N$&\multicolumn{6}{|c|}{\cellcolor{yellow}$100$}&\multicolumn{6}{|c|}{\cellcolor{yellow}$250$}\\\hline
&&\multicolumn{6}{c}{$K_N$}&\multicolumn{6}{c}{$K_N$}\\\hline
&$\beta_0-\hat{\beta_0}$&$3$&$4$&$5$&$6$&$7$& $8$ &$3$&$4$&$5$&$6$&$7$& $8$\\\midrule
\multicolumn{14}{c}{\bf Control function: $\widehat{h}(a_i)$}\\ \midrule
\multirow{4}{*}{\rotatebox[origin=c]{90}{$\exp(a_i)$}}&mean&1.287 &1.247 &1.279 &1.280 &1.288 &1.264 &1.172 &1.181 &1.188 &1.200 &1.185 &1.181 \\
&median&0.576 &0.551 &0.561 &0.562 &0.569 &0.568 &0.530 &0.534 &0.534 &0.538 &0.531 &0.531 \\
&std&1.864 &1.813 &1.887 &1.905 &1.889 &1.816 &1.673 &1.691 &1.702 &1.733 &1.698 &1.681 \\
&iqr&1.553 &1.499 &1.532 &1.543 &1.546 &1.537 &1.427 &1.436 &1.442 &1.450 &1.442 &1.441 \\\midrule
\multirow{4}{*}{\rotatebox[origin=c]{90}{$\cos(a_i)$}}&mean&1.877 &1.898 &1.883 &1.866 &1.925 &1.884 &1.793 &1.810 &1.809 &1.797 &1.800 &1.795 \\
&median&0.921 &0.931 &0.916 &0.922 &0.940 &0.916 &0.896 &0.904 &0.901 &0.897 &0.896 &0.889 \\
&std&2.528 &2.538 &2.528 &2.490 &2.608 &2.537 &2.351 &2.380 &2.380 &2.362 &2.373 &2.373 \\
&iqr&2.333 &2.402 &2.357 &2.344 &2.407 &2.357 &2.274 &2.282 &2.292 &2.274 &2.274 &2.261 \\\midrule
\multirow{4}{*}{\rotatebox[origin=c]{90}{$\sin(a_i)$}}&mean&1.433 &1.450 &1.454 &1.452 &1.490 &1.483 &1.360 &1.375 &1.362 &1.369 &1.367 &1.375 \\
&median&0.647 &0.653 &0.652 &0.665 &0.675 &0.666 &0.624 &0.632 &0.620 &0.631 &0.619 &0.625 \\
&std&2.050 &2.071 &2.072 &2.051 &2.144 &2.140 &1.911 &1.936 &1.915 &1.920 &1.946 &1.940 \\
&iqr&1.730 &1.762 &1.783 &1.773 &1.803 &1.799 &1.666 &1.680 &1.672 &1.680 &1.673 &1.673 \\\bottomrule
\toprule
\multicolumn{14}{c}{\bf Control function: $\widehat{h}(\widehat{deg_i},a_i)$}\\ \midrule
\multirow{4}{*}{\rotatebox[origin=c]{90}{$\exp(a_i)$}}&mean&1.930 &1.908 &1.879 &1.898 &1.977 &1.891 &1.625 &1.584 &1.666 &1.705 &1.636 &1.601 \\
&median&0.784 &0.775 &0.749 &0.756 &0.797 &0.762 &0.700 &0.682 &0.714 &0.711 &0.701 &0.689 \\
&std&3.124 &3.105 &3.203 &3.198 &3.371 &3.015 &2.437 &2.409 &2.538 &2.716 &2.482 &2.444 \\
&iqr&2.181 &2.166 &2.085 &2.120 &2.193 &2.152 &1.926 &1.860 &1.954 &1.965 &1.922 &1.889 \\\midrule
\multirow{4}{*}{\rotatebox[origin=c]{90}{$\cos(a_i)$}}&mean&2.555 &2.522 &2.527 &2.535 &2.576 &2.570 &2.225 &2.220 &2.268 &2.244 &2.219 &2.236 \\
&median&1.137 &1.135 &1.125 &1.148 &1.154 &1.142 &1.060 &1.054 &1.062 &1.056 &1.043 &1.043 \\
&std&3.854 &3.931 &3.956 &3.893 &3.900 &4.009 &3.088 &3.106 &3.235 &3.159 &3.143 &3.189 \\
&iqr&3.039 &2.957 &2.956 &2.990 &3.066 &3.014 &2.745 &2.724 &2.763 &2.749 &2.713 &2.721 \\\midrule
\multirow{4}{*}{\rotatebox[origin=c]{90}{$\sin(a_i)$}}&mean&2.058 &2.033 &2.053 &1.996 &2.093 &2.085 &1.755 &1.799 &1.768 &1.742 &1.805 &1.845 \\
&median&0.861 &0.838 &0.860 &0.846 &0.877 &0.878 &0.780 &0.797 &0.773 &0.774 &0.782 &0.795 \\
&std&3.244 &3.392 &3.216 &3.119 &3.315 &3.317 &2.560 &2.677 &2.622 &2.574 &2.769 &2.935 \\
&iqr&2.380 &2.317 &2.383 &2.327 &2.416 &2.416 &2.108 &2.144 &2.105 &2.080 &2.137 &2.156 \\\bottomrule
\end{tabular}}
\begin{tablenotes}\tiny
\item The statistics are based on conventional leave one out cross-validation.
\end{tablenotes}
\end{threeparttable}
\end{table}
\section{Conclusions}\label{section: conclusions}
In this paper we show that, whenever the network is likely endogenous, it is important to control for this endogeneity when estimating peer effects.
Failing to control for the endogeneity of the connections matrix in general leads to biased estimates of peer effects. We show that under specific assumptions, we can use the control function approach to deal with the endogeneity problem. We assume that unobserved individual characteristics directly affect link formation and individual outcomes. We leave the functional form through which unobserved individual characteristics enter the outcome equation unspecified and estimate it using a non-parametric approach. The estimators we propose are easy to use in applied work, and Monte Carlo results show that they perform well compared to a linear control function estimator. Erroneously assuming that unobserved characteristics enter the outcome equation in a linear fashion can lead to a serious bias in the estimated parameters.
{
\bibliographystyle{chicago}
\begin{thebibliography}{}
\bibitem[\protect\citeauthoryear{Abadie and Imbens}{Abadie and
Imbens}{2006}]{AbadieImbens2006}
Abadie, A. and G.~W. Imbens (2006).
\newblock Large sample properties of matching estimators for average treatment
effects.
\newblock {\em Econometrica\/}~{\em 74\/}(1), 235--267.
\bibitem[\protect\citeauthoryear{Ahn and Powell}{Ahn and
Powell}{1993}]{Ahn1993}
Ahn, H. and J.~L. Powell (1993).
\newblock Semiparametric estimation of censored selection models with a
nonparametric selection mechanism.
\newblock {\em Journal of Econometrics\/}~{\em 58\/}(1-2), 3--29.
\bibitem[\protect\citeauthoryear{Arduini, Patacchini, Rainone, et~al.}{Arduini
et~al.}{2015}]{Arduini2015a}
Arduini, T., E.~Patacchini, E.~Rainone, et~al. (2015).
\newblock Parametric and semiparametric iv estimation of network models with
selectivity.
\newblock Technical report, Einaudi Institute for Economics and Finance (EIEF).
\bibitem[\protect\citeauthoryear{Auerbach}{Auerbach}{2016}]{Auerbach2016}
Auerbach, E. (2016).
\newblock Identification and estimation of models with endogenous network
formation.
\newblock {\em Working paper\/}.
\bibitem[\protect\citeauthoryear{Banerjee, Chandrasekhar, Duflo, and
Jackson}{Banerjee et~al.}{2013}]{Banerjee2013}
Banerjee, A., A.~G. Chandrasekhar, E.~Duflo, and M.~O. Jackson (2013).
\newblock The diffusion of microfinance.
\newblock {\em Science\/}~{\em 341\/}(6144), 1236498.
\bibitem[\protect\citeauthoryear{Blume, Brock, Durlauf, and Ioannides}{Blume
et~al.}{2011}]{Blume2010}
Blume, L.~E., W.~A. Brock, S.~N. Durlauf, and Y.~M. Ioannides (2011).
\newblock Identification of social interactions.
\newblock In J.~Benhabib, A.~Bisin, and M.~Jackson (Eds.), {\em Handbook of
social economics}, Volume~1, pp.\ 853--964. Amsterdam: Elsevier.
\bibitem[\protect\citeauthoryear{Blume, Brock, Durlauf, and Jayaraman}{Blume
et~al.}{2015}]{blume2015linear}
Blume, L.~E., W.~A. Brock, S.~N. Durlauf, and R.~Jayaraman (2015).
\newblock Linear social interactions models.
\newblock {\em Journal of Political Economy\/}~{\em 123\/}(2), 444--496.
\bibitem[\protect\citeauthoryear{Bramoull\'{e}, Djebbari, and
Fortin}{Bramoull\'{e} et~al.}{2009}]{Bramoulle2009}
Bramoull\'{e}, Y., H.~Djebbari, and B.~Fortin (2009).
\newblock {Identification of peer effects through social networks}.
\newblock {\em Journal of Econometrics\/}~{\em 150\/}(1), 41--55.
\bibitem[\protect\citeauthoryear{Brock and Durlauf}{Brock and
Durlauf}{2001}]{brock2001interactions}
Brock, W.~A. and S.~N. Durlauf (2001).
\newblock Interactions-based models.
\newblock In J.~J. Heckman and E.~Leamer (Eds.), {\em Handbook of
econometrics}, Volume~5, pp.\ 3297--3380. Amsterdam: Elsevier.
\bibitem[\protect\citeauthoryear{Chen, Fern{\'a}ndez-Val, and Weidner}{Chen
et~al.}{2014}]{ChenFernandez-ValWeidner2018}
Chen, M., I.~Fern{\'a}ndez-Val, and M.~Weidner (2014).
\newblock Nonlinear factor models for network and panel data.
\newblock {\em arXiv preprint arXiv:1412.5647\/}.
\bibitem[\protect\citeauthoryear{De~Paula}{De~Paula}{2017}]{dePaula2016}
De~Paula, A. (2017).
\newblock Econometrics of network models.
\newblock In {\em Advances in Economics and Econometrics: Theory and
Applications, Eleventh World Congress}, pp.\ 268--323. Cambridge University
Press Cambridge.
\bibitem[\protect\citeauthoryear{De~Weerdt and Fafchamps}{De~Weerdt and
Fafchamps}{2011}]{fafchamps2011}
De~Weerdt, J. and M.~Fafchamps (2011).
\newblock Social identity and the formation of health insurance networks.
\newblock {\em Journal of Development Studies\/}~{\em 47\/}(8), 1152--1177.
\bibitem[\protect\citeauthoryear{Ductor, Fafchamps, Goyal, and van~der
Leij}{Ductor et~al.}{2014}]{Ductor2011}
Ductor, L., M.~Fafchamps, S.~Goyal, and M.~J. van~der Leij (2014).
\newblock Social networks and research output.
\newblock {\em Review of Economics and Statistics\/}~{\em 96\/}(5), 936--948.
\bibitem[\protect\citeauthoryear{Dzemski}{Dzemski}{2018}]{dzemski2017empirical}
Dzemski, A. (2018).
\newblock An empirical model of dyadic link formation in a network with
unobserved heterogeneity.
\newblock {\em forthcoming in Review of Economics and Statistics\/}.
\bibitem[\protect\citeauthoryear{Epple and Romano}{Epple and
Romano}{2011}]{Epple2011}
Epple, D. and R.~E. Romano (2011).
\newblock Peer effects in education: A survey of the theory and evidence.
\newblock In J.~Benhabib, A.~Bisin, and M.~Jackson (Eds.), {\em Handbook of
social economics}, Volume~1, pp.\ 1053--1163. Amsterdam: Elsevier.
\bibitem[\protect\citeauthoryear{Escanciano, Jacho-Ch{\'a}vez, and
Lewbel}{Escanciano et~al.}{2014}]{Escancianoetal2014}
Escanciano, J.~C., D.~T. Jacho-Ch{\'a}vez, and A.~Lewbel (2014).
\newblock Uniform convergence of weighted sums of non and semiparametric
residuals for estimation and testing.
\newblock {\em Journal of Econometrics\/}~{\em 178\/}(3), 426--443.
\bibitem[\protect\citeauthoryear{Fafchamps and Gubert}{Fafchamps and
Gubert}{2007}]{fafchamps2007risk}
Fafchamps, M. and F.~Gubert (2007).
\newblock Risk sharing and network formation.
\newblock {\em American Economic Review\/}~{\em 97\/}(2), 75--79.
\bibitem[\protect\citeauthoryear{Fern\'{a}ndez-Val and
Weidner}{Fern\'{a}ndez-Val and Weidner}{2013}]{Fernandez-val}
Fern\'{a}ndez-Val, I. and M.~Weidner (2013).
\newblock {Individual and time effects in nonlinear panel models with large N,
T}.
\newblock {\em arXiv preprint arXiv:1311.7065\/}.
\bibitem[\protect\citeauthoryear{Goldsmith-Pinkham and
Imbens}{Goldsmith-Pinkham and Imbens}{2013}]{GoldsmithP2013}
Goldsmith-Pinkham, P. and G.~W. Imbens (2013).
\newblock {Social Networks and the Identification of Peer Effects}.
\newblock {\em Journal of Business \& Economic Statistics\/}~{\em 31\/}(3),
253--264.
\bibitem[\protect\citeauthoryear{Graham}{Graham}{2011}]{Graham2011}
Graham, B.~S. (2011).
\newblock Econometric methods for the analysis of assignment problems in the
presence of complementarity and social spillovers.
\newblock In J.~Benhabib, A.~Bisin, and M.~Jackson (Eds.), {\em Handbook of
social economics}, Volume~1, pp.\ 965--1052. Amsterdam: Elsevier.
\bibitem[\protect\citeauthoryear{Graham}{Graham}{2017}]{Graham2017}
Graham, B.~S. (2017).
\newblock An econometric model of network formation with degree heterogeneity.
\newblock {\em Econometrica\/}~{\em 85\/}(4), 1033--1063.
\bibitem[\protect\citeauthoryear{Hahn and Ridder}{Hahn and
Ridder}{2013}]{HahnRidder2013}
Hahn, J. and G.~Ridder (2013).
\newblock Asymptotic variance of semiparametric estimators with generated
regressors.
\newblock {\em Econometrica\/}~{\em 81\/}(1), 315--340.
\bibitem[\protect\citeauthoryear{Hall and Heyde}{Hall and
Heyde}{2014}]{HallHeyde2014}
Hall, P. and C.~C. Heyde (2014).
\newblock {\em Martingale limit theory and its application}.
\newblock New York: Academic press.
\bibitem[\protect\citeauthoryear{Hansen}{Hansen}{2014}]{hansen2014nonparametric}
Hansen, B.~E. (2014).
\newblock Nonparametric sieve regression: Least squares, averaging least
squares, and cross-validation.
\newblock In J.~Racine, L.~Su, and A.~Ullah (Eds.), {\em Handbook of Applied
Nonparametric and Semiparametric Econometrics and Statistics}, pp.\
215--248. Oxford: Oxford University Press.
\bibitem[\protect\citeauthoryear{Heckman, Ichimura, and Todd}{Heckman
et~al.}{1998}]{Heckmanetal1998}
Heckman, J.~J., H.~Ichimura, and P.~Todd (1998).
\newblock Matching as an econometric evaluation estimator.
\newblock {\em Review of Economic Studies\/}~{\em 65\/}(2), 261--294.
\bibitem[\protect\citeauthoryear{Hsieh and Lee}{Hsieh and
Lee}{2016}]{hsieh2016social}
Hsieh, C.-S. and L.~F. Lee (2016).
\newblock A social interactions model with endogenous friendship formation and
selectivity.
\newblock {\em Journal of Applied Econometrics\/}~{\em 31\/}(2), 301--319.
\bibitem[\protect\citeauthoryear{Jackson}{Jackson}{2005}]{jackson2005survey}
Jackson, M.~O. (2005).
\newblock A survey of network formation models: stability and efficiency.
\newblock In G.~Demange and M.~Wooders (Eds.), {\em Group Formation in
Economics: Networks, Clubs, and Coalitions}, pp.\ 11--49. New York:
Cambridge University Press.
\bibitem[\protect\citeauthoryear{Jochmans}{Jochmans}{2016}]{jochmans2016modified}
Jochmans, K. (2016).
\newblock Modified-likelihood estimation of the b-model.
\newblock Technical report, Sciences Po Departement of Economics.
\bibitem[\protect\citeauthoryear{Jochmans}{Jochmans}{2018}]{jochmans2018semiparametric}
Jochmans, K. (2018).
\newblock Semiparametric analysis of network formation.
\newblock {\em Journal of Business \& Economic Statistics\/}~{\em 36\/}(4),
705--713.
\bibitem[\protect\citeauthoryear{Johnsson and Moon}{Johnsson and
Moon}{2019}]{JohnssonMoon2019}
Johnsson, I. and H.~R. Moon (2019).
\newblock Estimation of peer effects in endogenous social networks: Control
function approach.
\newblock {\em Working Paper, available from http://www-bcf.usc.edu/~moonr/\/}.
\bibitem[\protect\citeauthoryear{Kelejian and Prucha}{Kelejian and
Prucha}{1998}]{Kelejian1998}
Kelejian, H.~H. and I.~R. Prucha (1998).
\newblock A generalized spatial two-stage least squares procedure for
estimating a spatial autoregressive model with autoregressive disturbances.
\newblock {\em The Journal of Real Estate Finance and Economics\/}~{\em
17\/}(1), 99--121.
\bibitem[\protect\citeauthoryear{Lee}{Lee}{2003}]{Lee2003}
Lee, L. (2003).
\newblock {Best Spatial Two‐Stage Least Squares Estimators for a Spatial
Autoregressive Model with Autoregressive Disturbances}.
\newblock {\em Econometric Reviews\/}~{\em 22\/}(4), 307--335.
\bibitem[\protect\citeauthoryear{Lee}{Lee}{2007a}]{Lee2007a}
Lee, L.-F. (2007a).
\newblock {GMM and 2SLS estimation of mixed regressive, spatial autoregressive
models}.
\newblock {\em Journal of Econometrics\/}~{\em 137\/}(2), 489--514.
\bibitem[\protect\citeauthoryear{Lee}{Lee}{2007b}]{Lee2007}
Lee, L.-F. (2007b).
\newblock {Identification and estimation of econometric models with group
interactions, contextual factors and fixed effects}.
\newblock {\em Journal of Econometrics\/}~{\em 140\/}(2), 333--374.
\bibitem[\protect\citeauthoryear{Lee, Liu, and Lin}{Lee et~al.}{2010}]{Lee2010}
Lee, L.-f., X.~Liu, and X.~Lin (2010).
\newblock Specification and estimation of social interaction models with
network structures.
\newblock {\em The Econometrics Journal\/}~{\em 13\/}(2), 145--176.
\bibitem[\protect\citeauthoryear{Li}{Li}{1987}]{li1987asymptotic}
Li, K.-C. (1987).
\newblock Asymptotic optimality for cp, cl, cross-validation and generalized
cross-validation: discrete index set.
\newblock {\em Annals of Statistics\/}, 958--975.
\bibitem[\protect\citeauthoryear{Li et~al.}{Li
et~al.}{1987}]{andrews1991asymptotic}
Li, K.-C. et~al. (1987).
\newblock Asymptotic optimality for $ c\_p, c\_l $, cross-validation and
generalized cross-validation: Discrete index set.
\newblock {\em The Annals of Statistics\/}~{\em 15\/}(3), 958--975.
\bibitem[\protect\citeauthoryear{Li and Racine}{Li and Racine}{2007}]{Li2008a}
Li, Q. and S.~J. Racine (2007).
\newblock {\em {Nonparametric Econometrics: Theory and Practice}}.
\newblock Princeton University Press, Princeton.
\bibitem[\protect\citeauthoryear{Mammen, Rothe, Schienle, et~al.}{Mammen
et~al.}{2012}]{Mammenetal2012}
Mammen, E., C.~Rothe, M.~Schienle, et~al. (2012).
\newblock Nonparametric regression with nonparametrically generated covariates.
\newblock {\em The Annals of Statistics\/}~{\em 40\/}(2), 1132--1170.
\bibitem[\protect\citeauthoryear{Manski}{Manski}{1993}]{Manski1993a}
Manski, C.~F. (1993).
\newblock Identification of endogenous social effects: The reflection problem.
\newblock {\em Review of Economic Studies\/}~{\em 60\/}(3), 531--542.
\bibitem[\protect\citeauthoryear{Manski}{Manski}{2000}]{Manski2000}
Manski, C.~F. (2000).
\newblock Economic analysis of social interactions.
\newblock Technical report, NBER.
\bibitem[\protect\citeauthoryear{Newey}{Newey}{1997}]{Newey1997}
Newey, W.~K. (1997).
\newblock {Convergence rates and asymptotic normality for series estimators}.
\newblock {\em Journal of Econometrics\/}~{\em 79\/}(1), 147--168.
\bibitem[\protect\citeauthoryear{Newey}{Newey}{2009}]{Newey2009}
Newey, W.~K. (2009).
\newblock Two-step series estimation of sample selection models.
\newblock {\em The Econometrics Journal\/}~{\em 12\/}(s1), S217--S229.
\bibitem[\protect\citeauthoryear{Powell}{Powell}{1987}]{Powell1987}
Powell, J. (1987).
\newblock {\em Semiparametric estimation of bivariate latent variable models}.
\newblock University of Wisconsin--Madison, Social Systems Research Institute,
Madison.
\bibitem[\protect\citeauthoryear{Qu and Lee}{Qu and Lee}{2015}]{Qu2015}
Qu, X. and L.-F. Lee (2015).
\newblock {Estimating a spatial autoregressive model with an endogenous spatial
weight matrix}.
\newblock {\em Journal of Econometrics\/}~{\em 184\/}(2), 209--232.
\bibitem[\protect\citeauthoryear{Robinson}{Robinson}{1988}]{Robinson1988}
Robinson, P. (1988).
\newblock {Root-N-consistent semiparametric regression}.
\newblock {\em Econometrica: Journal of the Econometric Society\/}~{\em
56\/}(4), 931--954.
\bibitem[\protect\citeauthoryear{Shalizi}{Shalizi}{2012}]{Shalizi2012}
Shalizi, C.~R. (2012).
\newblock Comment on ``why and when `flawed' social network analyses still
yield valid tests of no contagion''.
\newblock {\em Statistics, politics, and policy\/}~{\em 3\/}(1), 5.
\bibitem[\protect\citeauthoryear{Sheng}{Sheng}{2012}]{Sheng2012}
Sheng, S. (2012).
\newblock Identification and estimation of network formation games.
\newblock {\em Unpublished manuscript\/}.
\bibitem[\protect\citeauthoryear{Wahba et~al.}{Wahba
et~al.}{1985}]{wahba1985comparison}
Wahba, G. et~al. (1985).
\newblock A comparison of gcv and gml for choosing the smoothing parameter in
the generalized spline smoothing problem.
\newblock {\em The Annals of Statistics\/}~{\em 13\/}(4), 1378--1402.
\bibitem[\protect\citeauthoryear{Weinberg}{Weinberg}{2007}]{Weinberg2007}
Weinberg, B.~A. (2007).
\newblock Social interactions with endogenous associations.
\newblock Working Paper 13038, National Bureau of Economic Research.
\end{thebibliography}
}