2.1 Pairwise strategic network formation model

Suppose that we have a sample of n agents that form social networks whose connections are represented by an $n \times n$ adjacency matrix $G_n = (g_{i,j})_{1 \le i,j \le n}$ . These agents can be individuals, firms, municipalities, or nations depending on the context. The network is directed; that is, regardless of the value of $g_{j,i}$ , we observe $g_{i,j} = 1$ if agent i links to j and $g_{i,j} = 0$ otherwise. There are no self-loops; that is, the diagonal elements of $G_n$ are all zero. Throughout the paper, we assume that the status of $(g_{i,j}, g_{j,i})$ is determined solely by the pair of agents (i,j), without considering the status of other network links. Specifically, for each pair (i,j), suppose that i’s marginal payoff of forming a link with j given $g_{j,i} = q$ is written as

\begin{align*} u_{i,j}(q) = Z_{i,j}^\top \beta_0 + \alpha_0 q + A_{0,i} + B_{0,j} - \epsilon_{i,j}, \;\; \text{for} \;\; i \neq j.\end{align*}

Here, $Z_{i,j} \in \mathbb{R}^{d_z}$ is a vector of observed covariates, $A_{0,i} \in \mathbb{R}$ is agent i’s individual specific effect as a “sender,” $B_{0,j}\in\mathbb{R}$ is j’s individual specific effect as a “receiver,” $\epsilon_{i,j} \in \mathbb{R}$ is an unobservable payoff component, and $\beta_0 \in \mathbb{R}^{d_z}$ and $\alpha_0 \in \mathbb{R}$ are unknown coefficient vector and the interaction effect parameter, respectively. The covariates $Z_{i,j}$ and $Z_{j,i}$ may contain common elements; however, in the later discussion, we require that they must have some agent-specific elements that can vary across the partners. The individual specific effects $A_{0,i}$ and $B_{0,j}$ , which we call the sender and the receiver effect, respectively, can be interpreted as the level of i’s willingness to create connections with others and the popularity of j, respectively, that generate degree heterogeneity across the agents. We treat $\{(A_{0,i}, B_{0,i})\}$ as fixed-effect parameters to be estimated.

We assume that the agents have complete information; that is, the realizations of $(Z_{i,j}, Z_{j,i})$ and $(\epsilon_{i,j}, \epsilon_{j,i})$ are common knowledge to both i and j. Then, if we assume that the observed network $G_n$ is formed by a collection of Nash equilibrium actions, we obtain the following econometric model:

(1) \begin{align}\begin{split} g_{i,j} & = \mathbf{1}\left\{ Z_{i,j}^\top \beta_0 + \alpha_0 g_{j,i} + A_{0,i} + B_{0,j} \ge \epsilon_{i,j} \right\} \\ g_{j,i} & = \mathbf{1}\left\{ Z_{j,i}^\top \beta_0 + \alpha_0 g_{i,j} + A_{0,j} + B_{0,i} \ge \epsilon_{i,j} \right\}, \;\; \text{for} \;\; i \neq j\end{split}\end{align}

The following are the two examples to which the above framework can be potentially applied.

In these examples, we can naturally imagine that the strategic interaction effect $\alpha_0$ is positive (i.e., strategic complements). We assume that strategic complementarity would be reasonable for most empirical situations of network formation games. Then, throughout the paper, we impose this assumption: $\alpha_0 > 0$ . Under strategic complementarity, each pair’s Nash equilibrium action can be summarized in Figure 1. As shown in the figure, the space of $(\epsilon_{i,j}, \epsilon_{j,i})$ cannot be partitioned into non-overlapping regions associated with the four alternative realizations of $(g_{i,j}, g_{j,i})$ . That is, both $(g_{i,j}, g_{j,i}) = (1,1)$ and $(g_{i,j}, g_{j,i}) = (0,0)$ can occur in the shaded area in the figure, and the link status is not uniquely determined in this area (i.e., multiple equilibria). This non-uniqueness of model-consistent decisions is called incompleteness and has been extensively studied in the literature on simultaneous equation models for discrete outcomes (e.g., Tamer, Reference Tamer2003; Lewbel, Reference Lewbel2007; Ciliberto & Tamer, Reference Ciliberto and Tamer2009; Chesher & Rosen, Reference Chesher and Rosen2020).

Figure 1.

Pure strategy of Nash equilibrium.

There are several approaches to handle this incompleteness issue in the literature.Footnote ² Among them, this study adopts the traditional approach developed by Bresnahan & Reiss (Reference Bresnahan and Reiss1990) and Berry (Reference Berry1992) that focuses only on the unique equilibrium outcomes. That is, we consider estimating the model using only the information about “one-way links” in the network. Specifically, we develop an ML estimator based on the probabilities of the three regions partitioned by the bold solid lines in Figure 1.

In the following, we assume that $\{\epsilon_{i,j}\}$ are identically distributed with a known cumulative distribution function (CDF) F. Further, we assume that the pairs $\{(\epsilon_{i,j}, \epsilon_{j,i})\}$ are independent and identically distributed (i.i.d.) across pairs, and their joint distribution is represented by $H(\cdot, \cdot;\; \rho_0)$ such that $\Pr(\epsilon_{i,j} \le a_1, \epsilon_{j,i} \le a_2) = H(a_1, a_2 ;\; \rho_0)$ , where $\rho_0 \in \mathbb{R}$ is a parameter controlling the correlation between $\epsilon_{i,j}$ and $\epsilon_{j,i}$ . Define $y^{(1,0)}_{i,j} \equiv \mathbf{1}\{(g_{i,j}, g_{j,i}) = (1,0)\}$ and $y^{(0,1)}_{i,j} \equiv \mathbf{1}\{(g_{i,j}, g_{j,i}) = (0,1)\}$ . Let $\chi_{n,a}$ be the a-th column of $I_n$ , $\mathbf{A}_0 = (A_{0,1}, \ldots, A_{0,n})^\top$ , and $\mathbf{B}_0 = (B_{0,1}, \ldots, B_{0,n})^\top$ , so that we can write $Z_{i,j}^\top \beta_0 + A_{0,i} + B_{0,j}= W_{i,j}^\top \Pi_0$ , where $W_{i,j} = (Z_{i,j}^\top, \chi_{n,i}^\top, \chi_{n,j}^\top)^\top$ and $\Pi_0 = (\beta_0^\top, \mathbf{A}_0^\top, \mathbf{B}_0^\top)^\top$ . In addition, we denote $\theta_0 = (\beta_0^\top, \alpha_0, \rho_0)^\top$ and $\gamma_0 = (\mathbf{A}_0^\top, \mathbf{B}_0^\top)^\top$ . Then, the probabilities of $\{y_{i,j}^{(1,0)} = 1\}$ and $\{y_{i,j}^{(0,1)} = 1\}$ are respectively given as follows:

\begin{align*} P^{(1,0)}_{i,j}(\theta_0, \gamma_0) & \equiv F(W_{i,j}^\top \Pi_0) - H(W_{i,j}^\top \Pi_0, W_{j,i}^\top \Pi_0 + \alpha_0;\; \rho_0), \\ P^{(0,1)}_{i,j}(\theta_0, \gamma_0) & \equiv F(W_{j,i}^\top \Pi_0) - H(W_{i,j}^\top \Pi_0 + \alpha_0, W_{j,i}^\top \Pi_0;\; \rho_0)\end{align*}

Here, note that the equalities $y^{(1,0)}_{i,j} = y^{(0,1)}_{j,i}$ and $P^{(1,0)}_{i,j}(\theta, \gamma) = P^{(0,1)}_{j,i}(\theta, \gamma)$ hold. Thus, the likelihood function can be concentrated with respect to $(y_{i,j}^{(1,0)}, P_{i,j}^{(1,0)}(\theta, \gamma))$ ; hereinafter, we omit the superscripts and denote $y_{i,j} = y_{i,j}^{(1,0)}$ and $P_{i,j}(\theta, \gamma) = P_{i,j}^{(1,0)}(\theta, \gamma)$ when there is no confusion. Then, the log-likelihood function can be written as

(2) \begin{align} \begin{split} \mathcal{L}_n(\theta, \gamma) & = \frac{2}{N}\sum_{i = 1}^n \sum_{j > i} \left[ y_{i,j} \ln P_{i,j}(\theta, \gamma) + y_{j,i} \ln P_{j,i}(\theta, \gamma) + (1 - y_{i,j} - y_{j,i}) \ln (1 - P_{i,j}(\theta, \gamma) - P_{j,i}(\theta, \gamma) )\right] \\ & = \frac{1}{N}\sum_{i = 1}^n \sum_{j \neq i} \left[ y_{i,j} \ln P_{i,j}(\theta, \gamma) + y_{j,i} \ln P_{j,i}(\theta, \gamma) + (1 - y_{i,j} - y_{j,i}) \ln (1 - P_{i,j}(\theta, \gamma) - P_{j,i}(\theta, \gamma) )\right] \\ & = \frac{1}{N} \sum_{i = 1}^n \sum_{j \neq i} \left[ 2 y_{i,j} \ln P_{i,j}(\theta, \gamma) + (1 - 2y_{i,j}) \ln (1 - P_{i,j}(\theta, \gamma) - P_{j,i}(\theta, \gamma) ) \right] \end{split}\end{align}

where $N \equiv n(n-1)$ . As above, we can consider three equivalent representations for the log-likelihood function and switch between them according to analytical convenience.

**framework (cf. Abadie et al., Reference Abadie, Athey, Imbens and Wooldridge2020). Particularly, Athey et al. (Reference Athey, Eckles and Imbens2018) considered the hypothesis testing of the causal effects under network interference using a finite-population network model. However, to the best of our knowledge, finite-population inference in the context of network formation models is still a theoretically unexplored topic.

2.2 Discrete heterogeneity

This study considers situations where the agents are grouped into several sub-samples, and the individual fixed effects are heterogeneous across these groups but are homogeneous within the groups in the following manner:

(3) \begin{align} \begin{split} A_{0,i} = \sum_{k = 1}^{K^A} a_{0,k} \cdot \mathbf{1}\{ i \in \mathcal{C}^A_{0,k}\}, \quad B_{0,i} = \sum_{k = 1}^{K^B} b_{0,k} \cdot \mathbf{1}\{ i \in \mathcal{C}^B_{0,k}\} \end{split}\end{align}

That is, the agents can be classified into $K^A$ groups $\mathcal{C}^A_0 \equiv \{\mathcal{C}^A_{0,1}, \ldots, \mathcal{C}^A_{0,K^A}\}$ in terms of the sender effects, where $K^A$ is the total number of groups, which form a partition of $\{1, \ldots, n\}$ into $K^A$ subsets. Similarly, in terms of the receiver effects, the agents can be grouped as $\mathcal{C}^B_0 \equiv \{\mathcal{C}^B_{0,1}, \ldots, \mathcal{C}^B_{0,K^B}\}$ . The group where each individual belongs to remains unknown to us. Meanwhile, it is often assumed in the literature that the number of groups is known to researchers (e.g., Bonhomme & Manresa, Reference Bonhomme and Manresa2015; Okui & Wang, Reference Okui and Wang2021). Then, following these studies, we treat $K^A$ and $K^B$ as known values and assume that $K^A, K^B \ge 2$ .Footnote ³ We discuss how to choose $K^A$ and $K^B$ in practice in Section 5. Note that transforming $(\mathbf{A}_0, \mathbf{B}_0)$ to $(\mathbf{A}_0 + c, \mathbf{B}_0 - c)$ for any constant c does not change the model (1). Thus, without loss of generality, we assume that $a_{0,1} = 0$ for normalization.

Under this setup, the full ML estimator solves

(4) \begin{align} \max_{\left(\theta, a_2, \ldots, a_{K^A}, b_1, \ldots, b_{K^B}, \mathcal{C}^A, \mathcal{C}^B \right) }\mathcal{L}_n(\theta, \mathbf{A}, \mathbf{B})\end{align}

where $\mathcal{C}^A \equiv \{\mathcal{C}^A_1, \ldots, \mathcal{C}^A_{K^A}\}$ , $\mathcal{C}^B \equiv \{\mathcal{C}^B_1, \ldots, \mathcal{C}^B_{K^B}\}$ , $A_i = \sum_{k = 1}^{K^A} a_k \cdot \mathbf{1}\{ i \in \mathcal{C}^A_k\}$ , and $B_i = \sum_{k = 1}^{K^B} b_k \cdot \mathbf{1}\{ i \in \mathcal{C}^B_k\}$ . The maximization problem in (4) is clearly a combinatorial (NP-hard) optimization problem. In the context of panel data models, several authors have proposed iterative algorithms to obtain a (local) solution efficiently to the problems similar to (4) (e.g., Bonhomme & Manresa, Reference Bonhomme and Manresa2015; Liu et al., Reference Liu, Shang, Zhang and Zhou2020). However, the iterative algorithm is still computationally demanding. More importantly, it cannot be directly applied to our network model where each agent’s heterogeneity parameters affect not only the value of his/her own likelihood function but also that of the others.

Hence, in this paper, we propose to break down the maximization problem in (4) into three steps. The first step is to estimate $\gamma_0 = (\mathbf{A}_0^\top, \mathbf{B}_0^\top)^\top$ using the full ML estimator based on the log-likelihood function in (2) without explicitly considering the group structure. Given the consistent estimates of these parameters, the second step is to estimate the group memberships $\mathcal{C}_0^A$ and $\mathcal{C}_0^B$ using the BS algorithm. The final step is to solve (4) with the group structure replaced by the estimated $\mathcal{C}_0^A$ and $\mathcal{C}_0^B$ .

2.3 Identification

Before presenting the estimation procedure in detail, we discuss identification conditions for the parameters in model (1). It is important to note that, even when individual heterogeneity parameters have only finite variations groupwisely, each individual’s parameters must be point-identified separately to estimate the group structure consistently. A reason for this is that our three-step estimator requires preliminary consistent estimates of $\mathbf{A}_0$ and $\mathbf{B}_0$ in the estimation of the group structure.Footnote ⁴ Here, again, we need some location normalization to identify $\mathbf{A}_0$ and $\mathbf{B}_0$ . Similarly as above, we set $A_{0,1} = 0$ without loss of generality. To facilitate the discussion, we also introduce several simplifying assumptions, some of which are mentioned previously.

Let $\Theta \equiv \mathcal{B} \times \mathcal{A} \times \mathcal{R}$ , $\mathbb{A}_n \equiv \{0\} \times \mathbb{A}^{n-1}$ , $\mathbb{B}_n \equiv \mathbb{B}^n$ , and $\mathbb{C}_n \equiv \mathbb{A}_n \times \mathbb{B}_n$ , where $\mathcal{B} \subset \mathbb{R}^{d_z}$ , $\mathcal{A} \subset \mathbb{R}_{++}$ , $\mathcal{R} \subset \mathbb{R}$ , $\mathbb{A} \subset \mathbb{R}$ , and $\mathbb{B} \subset \mathbb{R}$ are parameter spaces for $\beta$ , $\alpha$ , $\rho$ , $A_i$ ’s, and $B_i$ ’s, respectively.

In Assumption 2.1(i), we assume that the marginal CDF of the error term is known. This assumption is typically adopted in the estimation of complete information games. As shown by Khan & Nekipelov (Reference Khan and Nekipelov2018), when the marginal CDFs of $(\epsilon_{i,j}, \epsilon_{j,i})$ are unknown, it is generally impossible to estimate the interaction effect $\alpha_0$ at the parametric rate. Assumption 2.1(ii) requires that the error terms are independent across dyads. Note that the parameters $(A_{0,i}, B_{0,i})$ have the role of accommodating all unobserved payoff components specific to i. Therefore, assuming the independence within the remainders $\{\epsilon_{i,1} ,\ldots, \epsilon_{i,n} \}$ should not be too restrictive. The other requirements in Assumption 2.1 are standard in that they are satisfied in most of commonly used parametric models. In Assumption 2.2(ii), we assume that the fixed-effect parameters $\{(A_{0,i}, B_{0,i})\}$ are bounded. Although imposing boundedness on the degree heterogeneity parameters is commonly accepted in the literature on network formation models, some studies consider a more general framework where $||\mathbf{A}_0||_\infty$ and $||\mathbf{B}_0||_\infty$ can grow slowly (e.g., Yan et al., Reference Yan, Leng and Zhu2016; Yan et al., Reference Yan, Jiang, Fienberg and Leng2019). The parameter space $\mathcal{R}$ for the correlation parameter $\rho$ depends on the choice of the functional form of H, which is typically $\mathcal{R} = [-1,1]$ . Assumption 2.3 should not be restrictive in practice. Hereinafter, we fix the values of $\{Z_{i,j}\}$ ; that is, we interpret the following analysis as being conditional on the realization of $\{Z_{i,j}\}$ . Thus, any randomness in the model is considered to be due to the randomness of $\{\epsilon_{i,j}\}$ .

An important implication from Assumptions 2.1–2.3 is that the one-way link probabilities $\{P_{i,j}(\theta, \gamma)\}$ are uniformly bounded away from 0 and 1 for all possible parameter values. In other words, we are assuming that our networks are dense such that the number of one-way links per agent will be about proportional to the number of sampled agents. The plausibility of this assumption depends on the context of application. For example, international trade networks and the visa-free travel networks given in Example 2.2 and Section 5 may be regarded as dense networks.Footnote ⁵

Assumption 2.4 For any $(\beta_1, \alpha_1, \gamma_1), (\beta_2, \alpha_2, \gamma_2) \in \mathcal{B} \times \mathcal{A} \times \mathbb{C}_n$ such that $(\beta_1, \alpha_1, \gamma_1) \neq (\beta_2, \alpha_2, \gamma_2)$ , either or both (a) and (b) hold:

\begin{align*} a. \;\; & \liminf_{n \to \infty} \frac{1}{N} \sum_{i = 1}^n \sum_{j \neq i} \mathbf{1}\left\{ W_{i,j}^\top (\Pi_1 - \Pi_2) > 0, \;\; W_{j,i}^\top (\Pi_1 - \Pi_2) + \alpha_1 - \alpha_2 < 0 \right\} > 0\\ b. \;\; & \liminf_{n \to \infty} \frac{1}{N} \sum_{i = 1}^n \sum_{j \neq i} \mathbf{1}\left\{ W_{i,j}^\top (\Pi_1 - \Pi_2) < 0, \;\; W_{j,i}^\top (\Pi_1 - \Pi_2) + \alpha_1 - \alpha_2 > 0 \right\} > 0 \end{align*}

where $\Pi_1 = (\beta_1^\top, \gamma_1^\top)^\top$ , and $\Pi_2 = (\beta_2^\top, \gamma_2^\top)^\top$ .

Assumption 2.4 is our main identification condition, which basically requires the following two conditions. The first condition is a standard full-rank condition for $\{W_{i,j}\}$ and $\{W_{j,i}\}$ . The second condition is that at least either $Z_{i,j}$ or $Z_{j,i}$ should contain agent-specific covariates that have large enough supports and also have variations across all potential partners. If no such variables exist, since the signs of $Z_{i,j}^\top (\beta_1 - \beta_2)$ and $Z_{j,i}^\top (\beta_1 - \beta_2)$ cannot differ for some parameter values for all (i,j)’s, Assumption 2.4 does not hold. It should be noted that this assumption is not inconsistent with Assumption 2.3 under the compactness of the parameter space (i.e., Assumption 2.2). While the existence of player-specific continuous variables with unbounded supports is typically required in the identification of non/semiparametric game models—the so-called identification-at-infinity argument (see, e.g., Tamer, Reference Tamer2003; Kline, Reference Kline2015), we can develop our identification result under less stringent conditions owing to the full-parametric model specification.

These identification results are similar to those in Theorem 2 of Aradillas-Lopez & Rosen (Reference Aradillas-Lopez and Rosen2019). In the literature on network formation models, it is often assumed that $\epsilon_{i,j}$ and $\epsilon_{j,i}$ are independent (e.g., Hoff et al., Reference Hoff, Raftery and Handcock2002; Jochmans, Reference Jochmans2018; Yan et al., Reference Yan, Jiang, Fienberg and Leng2019). If they are independent, since $\rho_0 = 0$ is known, condition (i) is satisfied. Condition (ii) clearly depends on the choice of H function and is difficult to verify in general; however, this is directly empirically testable. A more primitive sufficient condition for this is that $\mathcal{L}_n^*(\rho)$ is strictly concave, which is satisfied when $\partial^2 \mathcal{L}_n^*(\rho)/(\partial \rho)^2$ is strictly negative under Assumption 2.1(iii). We provide an explicit form of $\partial^2 \mathcal{L}_n^*(\rho)/(\partial \rho)^2$ in Appendix C.1. For another identification condition for $\rho_0$ , see, for example, Theorem B.2 of Hoshino & Yanagi (Reference Hoshino and Yanagi2021), that is based on the monotonicity of the likelihood function with respect to $\rho$ .

**(e.g., Rothenberg, Reference Rothenberg1971; Bjorn & Vuong, Reference Bjorn and Vuong1984), although it is not easy to verify in practice. If one admits the existence of a player-specific continuous variable that has a positive density on the whole $\mathbb{R}$ , then the model can be easily identified by the identification-at-infinity approach in the same manner as in Tamer (Reference Tamer2003). Since assuming the existence of such unbounded variables is restrictive in practice, several authors have proposed other approaches based on some shape restrictions on the distribution of unobservables (e.g., Kline, Reference Kline2016). Investigating whether their approaches can be applied to our model is an interesting topic, but it is left for a future work.

3.1 First-step ML estimator

The first step of the ML estimation aims to obtain consistent estimates of $\gamma_0 = (\mathbf{A}_0^\top, \mathbf{B}_0^\top)^\top$ . Let

(5) \begin{align}\begin{split} (\widehat \theta_n, \widehat \gamma_n) = \mathrm{argmax}_{(\theta, \gamma) \: \in \: \Theta \times \mathbb{C}_n} \mathcal{L}_n(\theta, \gamma)\end{split}\end{align}

where $\widehat \theta_n = (\widehat \beta_n^\top, \widehat \alpha_n, \widehat \rho_n)^\top$ , and $\widehat \gamma_n = (\widehat{\mathbf{A}}_n^\top, \widehat{\mathbf{B}}_n^\top)^\top$ . Below, we present the asymptotic properties of the initial full ML estimator in (5) with the main focus on $\widehat \gamma_n$ . Instead of introducing particular identification conditions, for generality, we directly assume that the true parameter $(\theta_0, \gamma_0)$ is a unique maximizer of $\mathbb{E} \mathcal{L}_n(\theta, \gamma)$ .

We first establish several consistency results in the next theorem.

Next, we derive the convergence rate of $\widehat \gamma_n$ . To this end, it is convenient to re-define $\widehat \theta_n$ and $\theta_0$ as

\begin{align*} \widehat \theta_n = \mathrm{argmax}_{\theta \in \Theta} \mathcal{L}_n(\theta, \widetilde{\gamma}_n(\theta)) \;\; \text{and} \;\; \theta_0 = \mathrm{argmax}_{\theta \in \Theta} \mathbb{E} \mathcal{L}_n(\theta, \widetilde{\gamma}_0(\theta))\end{align*}

respectively, where $\widetilde \gamma_n(\theta) \equiv \mathrm{argmax}_{\gamma \in \mathbb{C}_n} \mathcal{L}_n(\theta, \gamma)$ and $\widetilde \gamma_0(\theta) \equiv \mathrm{argmax}_{\gamma \in \mathbb{C}_n} \mathbb{E} \mathcal{L}_n(\theta, \gamma)$ for any given $\theta \in \Theta$ . Further, we define $\gamma_{-1} = (A_2, \ldots, A_n, B_1, \ldots, B_n)^\top$ ,

\begin{align*} \underset{(2n-1) \times (2n-1)}{\mathcal{H}_{n, \gamma\gamma} (\theta, \gamma)} \equiv \frac{\partial^2 \mathcal{L}_n(\theta, \gamma)}{ \partial \gamma_{-1} \partial \gamma_{-1}^\top }, \qquad \underset{(d_z + 2) \times (d_z + 2)}{\mathcal{H}_{n, \theta\theta} (\theta, \gamma)} \equiv \frac{\partial^2 \mathcal{L}_n(\theta, \gamma)}{ \partial \theta \partial \theta^\top} , \qquad \underset{(2n-1) \times (d_z + 2)}{\mathcal{H}_{n, \gamma\theta} (\theta, \gamma)} \equiv \frac{\partial^2 \mathcal{L}_n(\theta, \gamma)}{ \partial \gamma_{-1} \partial \theta^\top}\end{align*}

and $\mathcal{H}_{n, \theta\gamma} (\theta, \gamma) \equiv \mathcal{H}_{n, \gamma\theta} (\theta, \gamma)^\top$ . The exact form of $\mathcal{H}_{n, \gamma\gamma} (\theta, \gamma)$ is given in (A.2) in Appendix A. Finally, let

\begin{align*} \underset{(d_z + 2) \times (d_z + 2)}{\mathcal{I}_{n,\theta\theta}(\theta, \gamma)} \equiv \mathcal{H}_{n, \theta \theta}(\theta, \gamma) - \mathcal{H}_{n, \theta \gamma} (\theta, \gamma) \left[ \mathcal{H}_{n,\gamma\gamma} (\theta, \gamma) \right]^{-1} \mathcal{H}_{n, \gamma \theta } (\theta, \gamma)\end{align*}

This $\mathcal{I}_{n,\theta\theta}(\theta, \gamma)$ serves as the Hessian matrix for the concentrated ML estimator $\widehat \theta_n$ (cf. Amemiya, Reference Amemiya1985).

Assumptions 3.2 and 3.3 should be fairly reasonable in practice. Then, under these additional assumptions, we can derive the $\ell_1$ -norm and max-norm convergence rates for $\widehat \gamma_n$ , as shown in the next theorem.

The max-norm convergence rate obtained in Theorem 3.2 is consistent with the result of Theorem 3 in Graham (Reference Graham2017) and that of Theorem 3.1 in Yan et al. (Reference Yan, Jiang, Fienberg and Leng2019).

3.2 Binary segmentation algorithm

Given the consistent estimates of $\mathbf{A}_0$ and $\mathbf{B}_0$ , we use the BS algorithm to estimate the group structure. We first sort $\widehat{\mathbf{A}}_n$ and $\widehat{\mathbf{B}}_n$ in ascending order and write the order statistics as $\widehat A_{n,(1)} \le \widehat A_{n,(2)} \le \cdots \le \widehat A_{n,(n)}$ , and $\widehat B_{n,(1)} \le \widehat B_{n,(2)} \le \cdots \le \widehat B_{n,(n)}$ . In the following, we mainly describe the estimation of the group membership for the sender effects, $\mathcal{C}^A_0$ . The exactly same procedure described below can be used to estimate $\mathcal{C}^B_0$ .

A concept behind the BS algorithm is quite simple. If $\{A_{0,i}\}$ are heterogeneous across $K^A$ latent groups but are homogeneous within the groups, there should exist $K^A - 1$ “break points” in the sorted $\{A_{0,i}\}$ . Since $\widehat{\mathbf{A}}_n$ is uniformly consistent for $\mathbf{A}_0$ , these break points appear also in the following sequence: $\widehat A_{n,(1)}, \ldots , \widehat A_{n,(n)}$ w.p.a.1. For $1 \le i < j \le n$ , we define $\widehat \Delta^A(i,j)$ as the sum of squared variations over $\{\widehat A_{n,(i)}, \ldots, \widehat A_{n,(j)}\}$ :

\begin{align*} \widehat \Delta^A(i,j) & \equiv \sum_{l = i}^j (\widehat A_{n,(l)} - \bar A_{n,i,j})^2, \;\; \text{where} \;\; \bar A_{n,i,j} \equiv \frac{1}{j - i + 1}\sum_{l = i}^j \widehat A_{n,(l)}\end{align*}

Further, we define

\begin{align*} \widehat S^A_{i,j}(\kappa) & \equiv \left\{ \begin{array}{ll} \frac{1}{j - i + 1} \left( \widehat \Delta^A(i,\kappa) + \widehat \Delta^A(\kappa + 1, j) \right) & \text{if} \;\; j < \kappa \\ \frac{1}{j - i + 1} \widehat \Delta^A(i,j) & \text{if} \;\; j = \kappa \end{array}\right.\end{align*}

That is, $\widehat S^A_{i,j}(\kappa)$ provides the total variance of $\{\widehat A_{n,(i)}, \ldots, \widehat A_{n,(j)}\}$ when a break point is placed at $\kappa$ . Assuming that $K^A \ge 2$ , the BS algorithm proceeds as follows:

Step 1 $(K^A = 2)$ :

We find the first break point, say $\widehat t_1$ , by

\begin{align*} \widehat t_1 = \mathrm{argmin}_{1 \le \kappa < n} \widehat S^A_{1,n}(\kappa)\end{align*}

Then, we can partition $\{\widehat A_{n,(i)}\}$ into the following two subsets: $\{\widehat A_{n,(i)}\} = \{\widehat A_{n,(i)}\}_{i=1}^{\widehat t_1} \bigcup \{\widehat A_{n,(i)}\}_{i = \widehat t_1 + 1}^n$ . If $K^A = 2$ , assuming that $a_{0,1} < a_{0,2}$ without loss of generality, we obtain $\widehat{\mathcal{C}}_{n,1}^A \equiv \{i : \widehat A_{n,(1)} \le \widehat A_{n,i} \le \widehat A_{n,(\widehat t_1)}\}$ and $\widehat{\mathcal{C}}_{n,2}^A \equiv \{i : \widehat A_{n,(\widehat t_1 + 1)} \le \widehat A_{n,i} \le \widehat A_{n,(n)}\}$ as the estimators of $\mathcal{C}_{0,1}^A$ and $\mathcal{C}_{0,2}^A$ , respectively, and the algorithm stops. If $K^A > 2$ , we proceed to the next step.

Step 2 $(K^A = 3)$ :

Now, if $K^A = 3$ , there exists one more break point either in $\{\widehat A_{n,(i)}\}_{i=1}^{\widehat t_1}$ or in $\{\widehat A_{n,(i)}\}_{i = \widehat t_1 + 1}^n$ . Either one of the two converges to a sequence of identical constants as the sample size increases. Therefore, if $\widehat S^A_{1,\widehat t_1}(\widehat t_1) > \widehat S^A_{\widehat t_1 + 1, n}(n)$ for example, the break point is likely to lie in the former subset. Thus, the second break point, say $\widehat t_2$ , can be found by

\begin{align*} \widehat t_2 = \begin{cases} \mathrm{argmin}_{1 \le \kappa < \widehat t_1} \widehat S^A_{1,\widehat t_1}(\kappa) & \text{if} \;\; \widehat S^A_{1,\widehat t_1}(\widehat t_1) > \widehat S^A_{\widehat t_1 + 1, n}(n) \\[5pt] \mathrm{argmin}_{\widehat t_1 + 1 \le \kappa < n} \widehat S^A_{\widehat t_1 + 1, n}(\kappa) & \text{if} \;\; \widehat S^A_{1,\widehat t_1}(\widehat t_1) \le \widehat S^A_{\widehat t_1 + 1, n}(n) \end{cases}\end{align*}

Then, we add $\widehat t_2$ to the set of break points, and (with a little abuse of notation) we sort and re-label the points to ensure that $\widehat t_1 < \widehat t_2$ . The resulting partition is $\{\widehat A_{n,(i)}\} = \{\widehat A_{n,(i)}\}_{i=1}^{\widehat t_1} \bigcup \{\widehat A_{n,(i)}\}_{i = \widehat t_1 + 1}^{\widehat t_2}\bigcup \{\widehat A_{n,(i)}\}_{i = \widehat t_2 + 1}^n$ . Setting $a_{0,1} < a_{0,2} < a_{0,3}$ without loss of generality, we can define $\widehat{\mathcal{C}}_{n,1}^A \equiv \{i : \widehat A_{n,(1)} \le \widehat A_{n,i} \le \widehat A_{n,(\widehat t_1)}\}$ , $\widehat{\mathcal{C}}_{n,2}^A \equiv \{i : \widehat A_{n,(\widehat t_1 + 1)} \le \widehat A_{n,i} \le \widehat A_{n,(\widehat t_2)}\}$ , and $\widehat{\mathcal{C}}_{n,3}^A \equiv \{i : \widehat A_{n,(\widehat t_2 + 1)} \le \widehat A_{n,i} \le \widehat A_{n,(n)}\}$ as the estimators of $\mathcal{C}_{0,1}^A$ , $\mathcal{C}_{0,2}^A$ , and $\mathcal{C}_{0,3}^A$ , respectively.

Step 3 $(K^A \ge 4)$ :

When $K^A = 4$ , if $\max\{\widehat S^A_{1,\widehat t_1}(\widehat t_1), \widehat S^A_{\widehat t_1 + 1, \widehat t_2}(\widehat t_2), \widehat S^A_{\widehat t_2 + 1, n}(n)\} = \widehat S^A_{1,\widehat t_1}(\widehat t_1)$ for example, the third break point should exist in $\{\widehat A_{n,(i)}\}_{i=1}^{\widehat t_1}$ . Then, following the same procedure as above, we can obtain $\widehat t_3 = \mathrm{argmin}_{1 \le \kappa < \widehat t_1} \widehat S^A_{1,\widehat t_1}(\kappa) $ . We repeat these steps until $K^A$ groups are detected with $K^A - 1$ break points. Finally, letting $\widehat t_0 = 0$ and $\widehat t_{K^A} = n$ so that $\widehat t_0 < \widehat t_1 < \cdots < \widehat t_{K^A - 1} < \widehat t_{K^A}$ , each k-th group can be estimated by $\widehat{\mathcal{C}}^A_{n,k} \equiv \{i : \widehat A_{n,(\widehat t_{k - 1} + 1)} \le \widehat A_{n,i} \le \widehat A_{n,(\widehat t_k)}\}$ . Then, $\widehat{\mathcal{C}}^A_n \equiv \{\widehat{\mathcal{C}}^A_{n,1}, \ldots, \widehat{\mathcal{C}}^A_{n, K^A}\}$ is our estimator for $\mathcal{C}_0^A$ .

In the same manner as above, we define $\widehat{\mathcal{C}}^B_n \equiv \{\widehat{\mathcal{C}}^B_{n,1}, \ldots, \widehat{\mathcal{C}}^B_{n, K^B}\}$ for the estimation of $\mathcal{C}_0^B$ . Note that the identification of group membership can be achieved only up to “label swapping.” Thus, without loss of generality, we can label the groups according to the values of the group effects such that $a_{0,1} < a_{0,2} <\cdots < a_{0,K^A}$ ; that is, we define the k-th group as the group with the k-th smallest sender effect. Similarly, we set $b_{0,1} < b_{0,2} <\cdots < b_{0,K^B}$ . To investigate the asymptotic properties of the BS algorithm, we introduce the following assumption.

Assumption 3.4 is parallel to Assumption A2 in Wang & Su (Reference Wang and Su2021). Similar assumptions are commonly used in the literature of panel data models with latent group structure. The following theorem provides the consistency result of the BS algorithm:

Although the proof of Theorem 3.3 is almost analogous to Ke et al. (Reference Ke, Li and Zhang2016) and Wang & Su (Reference Wang and Su2021), for completeness, we provide it in Appendix B.

**by Su et al. (Reference Su, Shi and Phillips2016). As discussed in Wang & Su (Reference Wang and Su2021), the main advantage of the BS algorithm over these alternatives is its computational and interpretational simplicity. The k-means algorithm is NP-hard, and C-Lasso requires solving a large nonlinear optimization problem, whereas our BS method does not involve any computational optimization.Footnote ⁶ Additionally, the C-Lasso may end up leaving some unclassified agents. Therefore, we need to classify such observations afterward based on some criterion. Another possibility is to use finite-mixture modeling, which is more common and traditional in statistical literature. However, the mixture-based approach is more suitable for cross-sectional data. It only estimates the “share” of each group, requiring an additional classification step to estimate the group membership of each agent, as in C-Lasso. Additionally, when $(K^A, K^B)$ is large, the mixture approach may be less practical because the marginal probability of each dyad is represented as a mixture of $(K^A K^B)^2$ conditional probabilities. As mentioned in the introduction, if the group structure is not of interest, the standard fixed-effect bias correction method should be also considered. For more details about these approaches, including some simulation results, see Appendices C.3 and C.4.

3.3 Post-grouping ML estimator

In the final step, we solve (4) approximately by using $(\widehat{\mathcal{C}}^A_n, \widehat{\mathcal{C}}^B_n)$ in the place of $(\mathcal{C}_0^A, \mathcal{C}_0^B)$ . Recalling that $a_1$ is pinned at $a_1 = 0$ , let $\delta \equiv (\theta^\top, a_2, \ldots, a_{K^A}, b_1, \ldots, b_{K^B})^\top$ , and $\mathbb{D}\equiv \Theta \times \mathbb{A}^{K^A - 1} \times \mathbb{B}^{K^B}$ for the parameter space of $\delta$ . We denote $\delta_0$ as the true value of $\delta$ . Then, our final ML estimator for $\delta_0$ is defined as

\begin{align*} \widehat \delta_n = \mathrm{argmax}_{\delta \in \mathbb{D}} \widehat{\mathcal{L}}_n(\delta)\end{align*}

where $\widehat{\mathcal{L}}_n(\delta) \equiv \mathcal{L}_n\bigl(\theta, \bigl\{ \sum_{k = 1}^{K^A} a_k \cdot \mathbf{1}\{ i \in \widehat{\mathcal{C}}^A_{n,k}\} \bigr\}, \bigl\{ \sum_{k = 1}^{K^B} b_k \cdot \mathbf{1}\{ i \in \widehat{\mathcal{C}}^B_{n,k}\} \bigr\} \bigr)$ . Similarly, we define

\begin{align*} \widehat \delta_n^{\text{oracle}} = \mathrm{argmax}_{\delta \in \mathbb{D}} \mathcal{L}_n(\delta)\end{align*}

where $\mathcal{L}_n(\delta) \equiv \mathcal{L}_n\bigl(\theta, \bigl\{ \sum_{k = 1}^{K^A} a_k \cdot \mathbf{1}\{ i \in \mathcal{C}^A_{0,k}\} \bigr\}, \bigl\{ \sum_{k = 1}^{K^B} b_k \cdot \mathbf{1}\{ i \in \mathcal{C}^B_{0,k}\} \bigr\} \bigr)$ ; that is, $\widehat \delta_n^{\text{oracle}}$ is the “oracle” estimator that is computed based on the true $\mathcal{C}_0^A$ and $\mathcal{C}_0^B$ . Since $\widehat \delta_n^{\text{oracle}}$ is a standard parametric ML estimator, the estimator follows a normal distribution asymptotically with its asymptotic covariance matrix given by the inverse Fisher Information matrix. Meanwhile, we have shown in Theorem 3.3 that the estimated group memberships $(\widehat{\mathcal{C}}^A_n, \widehat{\mathcal{C}}^B_n)$ are equal to $(\mathcal{C}_0^A, \mathcal{C}_0^B)$ w.p.a.1. Therefore, the final ML estimator $\widehat \delta_n$ has asymptotically the same statistical performance as the oracle estimator $\widehat \delta_n^{\text{oracle}}$ . We formally state this result in the next theorem.

Recall that the asymptotic equivalence between $\widehat \delta_n$ and $\widehat \delta_n^{\text{oracle}}$ relies on the dense network structure where each agent’s specific effects can be point-identified. When the networks are not dense, $\widehat \delta_n$ is generally inconsistent, while $\widehat \delta_n^{\text{oracle}}$ may be still consistent (potentially with a slower convergence rate). Finally, note that the above discussions hold true if the repartitioned estimator $(\widehat{\mathcal{C}}^{A, \mathrm{repart}}_n, \widehat{\mathcal{C}}^{B, \mathrm{repart}}_n)$ is used instead of $(\widehat{\mathcal{C}}^A_n, \widehat{\mathcal{C}}^B_n)$ .