Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
84,545 characters · 12 sections · 87 citation commands
A Pairwise Strategic Network Formation Model with Group Heterogeneity: With an Application to International Travel
Empirical modeling of network formation is an important research topic that has been studied for several decades. While most of these models has been developed in the mathematical statistics literature, as the importance of network structure in many economic activities has been increasingly recognized, there is currently a growing number of econometric studies that focus on network formation in conjunction with the significant advancement in the related econometric techniques.\footnote{ For recent developments regarding econometric approaches for analyzing network formation, we refer readers to, for example, chandrasekhar2016econometrics and de2020econometric. }
Econometric studies on network formation can be classified into two types: those that attempt to explicitly incorporate the interaction of individuals in the realizing network structure endogenously affecting the network formation behavior (e.g., leung2015two; mele2017structural; sheng2020structural) and those that do not account for such simultaneous interactions but emphasize modeling a flexible form of individual heterogeneity (e.g., graham2017econometric; jochmans2018semiparametric; dzemski2019empirical). For the former type, network formation is modeled as a game in which agents strategically form links to maximize their payoffs. Although this game-theoretic approach is (economic) theoretically well-underpinned, we often encounter serious analytical difficulties due to the presence of multiple equilibria. To circumvent these difficulties, we typically need to introduce some ad hoc behavioral assumptions into the network formation process, or we simply discontinue point-identifying the models and resort to partial identification. Compared with the former, the latter is more “descriptive” than “structural”, but has great flexibility in the model specification. These models are relatively easy to implement and, thus, are appealing to empirical researchers. However, they are not suitable for analyzing the interactions of agents in network formation, which should be an essential factor in economic and social network data.
These two types of econometric models have their own advantages. Hence, it is an ingenuous idea to construct a new model that has the advantages of both approaches by combining them. However, to my best knowledge, there are only a few papers that address this way of model extension (e.g., graham2016homophily; graham2020testing; pelican2020optimal). graham2016homophily gets around the multiple equilibria issue by considering the “dynamic” (rather than the instantaneous) interdependencies in network formation. The latter two papers consider a quite general framework that incorporates both a general form of strategic interaction and unobserved degree heterogeneity into a single model. However, they mainly focus on testing the presence of interactions, but not on the estimation of the models.
In this paper, we propose a new “pairwise” network formation model that is empirically tractable while retaining the nice properties of the above-mentioned approaches. More specifically, we assume that each link connection is determined by the strategic interaction solely between the corresponding dyad of agents, without affecting or being affected by other dyads, rather than regarding the realized network as a consequence of a large $n$-player game. Although ignoring such network externalities limits the range of applications of our model, it would still cover a fairly large number of interesting empirical situations, and, most importantly, the multiple equilibria problem can be greatly mitigated. In addition, we allow each agent's payoff to depend on two unobserved preference parameters that represent his/her outgoing and incoming propensity, which we refer to as the sender and the receiver effect, respectively.
For estimating our pairwise network formation model, assuming that the model error terms follow some parametric distribution, we propose using the maximum likelihood (ML) method. Although multiple equilibria can occur even in our binary game situation, we can address this issue in the same manner as in bresnahan1990entry and berry1992estimation. As in dzemski2019empirical and yan2019statistical, we treat the agent-specific preference heterogeneity parameters as the fixed-effect parameters to be estimated. However, note that if we include these heterogeneity parameters directly into the model, as the dimension of the parameters increases proportionally to the sample size, our ML estimator suffers from the incidental parameter problem. As will be confirmed in the numerical simulations reported later, the incidental parameter bias can be severe. Thus, to avoid the incidental parameter problem, as a second novel part of this study, we focus on a situation where the individuals can be classified into several unknown groups and their specific effects are homogeneous within each group. If the number of the latent groups is fixed, we can expect that the ML estimator becomes asymptotically unbiased at the parametric rate and hence the standard inference procedure can be applied.
In the literature on statistical network data analysis, uncovering such latent group structures in networks has been intensively studied (e.g., newman2004finding; bickel2009nonparametric; fortunato2010community; karrer2011stochastic; rohe2011spectral; abbe2017community). Indeed, putting the strategic interaction effect aside, our network formation model can be regarded as a type of the stochastic block model, one of the major modeling approaches in the above literature, with additional covariates, as in zhang2019node. In panel data analysis also, identification of unobserved grouped heterogeneity is one of the most active research areas (e.g., bonhomme2015grouped; ke2016structure; su2016identifying; wang2018homogeneity). Although the application of these grouping methods to econometric network models has been relatively limited thus far, it is a promising approach, as discussed in bonhomme2020econometric. If we can recover the true group structure with probability approaching one (w.p.a.1) by using any method, the individual fixed-effect parameters can be estimated at a faster rate than the case of no grouping. In this study, among several alternative methods, we adopt the binary segmentation (BS) method (see, e.g., bai1997estimating; ke2016structure; lian2019homogeneity; wang2020identifying). Compared to the other grouping methods, the BS method has several favorable properties including fast computation speed and robustness as it does not require us to set initial values.
The whole estimation procedure is divided into three steps. The first step is to obtain the ML estimator without considering the latent group structures. Although this estimator suffers from the incidental parameter bias, it is still possible to produce consistent estimates for the heterogeneity parameters. In the second step, we apply the BS method with respect to the estimated sender effect parameters and the receiver effect parameters separately to identify each agent's group memberships. In the third step, we re-estimate the model using the ML method given the estimated group structure. Under certain regularity conditions, we show that the proposed estimator is asymptotically unbiased and normal at the parametric rate. Furthermore, the estimator is asymptotically equivalent to the “oracle” estimator that is obtained based on the (unknown) true group memberships.
To illustrate our model empirically, we investigate the formation of international visa-free travel networks, where the dependent variable of interest is defined as follows: $g_{i,j} = 1$ if country $i$ allows the citizens in country $j$ to visit $i$ without visas and $g_{i,j} = 0$ if not. We apply our model framework to the network of 57 countries selected mainly from Asia, the Middle East, the former USSR, and Oceania. As expected, we find the presence of a certain level of degree heterogeneity in terms of both the sender and the receiver effects. Interestingly, there seems to be a negative correlation between the sender effects and the receiver effects (in other words, there is a tendency that a country's sender effect increases as its receiver effect decreases). Our estimation result also suggests that there is a significant strategic complementarity in the network formation behavior. Another interesting finding is that the countries are homophilous -- tend to connect with similar others -- in terms of the political system. These findings would highlight the usefulness of the proposed model and method.
\paragraph{Organization of the paper:} The remainder of the paper is organized as follows. In Section (ref), we formally introduce the model investigated in this study. In this section, we demonstrate that our pairwise model exhibits multiple equilibria and discuss the conditions under which the model can be point-identified. Section (ref) provides a detailed explanation about our three-step ML estimator. We also investigate the asymptotic properties of the proposed estimator in this section. In Section (ref), we present a set of Monte Carlo experiments to evaluate the finite sample performance of the proposed estimator. Section (ref) presents our empirical analysis, and, finally, Section (ref) concludes. All the technical details are relegated to Appendix.
\paragraph{Notation:} For a natural number $n$, $I_n$ denotes an $n \times n$ identity matrix. $\mathbf{1}\{\cdot\}$ denotes the indicator function, which is one if its argument is true and zero otherwise. For a matrix $A$, we use $||A||$ to denote its Frobenius norm: $||A|| =\sqrt{\mathrm{tr}\{AA^\top\}}$, where $\mathrm{tr}\{\cdot\}$ is a trace of a matrix. When $A$ is a square matrix, we use $\lambda_{\mathrm{min}}(A)$ to denote its smallest eigenvalue. For a vector $\mathbf{a} = (a_1, \ldots, a_k)^\top$, $||a||_\infty$ denotes its maximum norm: $||\mathbf{a}||_\infty = \max_{1 \le i \le k}|a_i|$. For a general set $\mathcal{X}$, we use $\mathcal{X}^{\mathrm{int}}$ to denote its interior. In addition, $|\mathcal{X}|$ denotes the cardinality of $\mathcal{X}$. $c$ (possibly with subscript) denotes a generic positive constant whose exact value may vary per case.
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
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. Following dzemski2019empirical and yan2019statistical, 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:
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 (ref). 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 exist). 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., tamer2003incomplete; lewbel2007coherency; ciliberto2009market; chesher2020structural).
There are several approaches to handle this incompleteness issue in the literature.\footnote{ Refer to de2013econometric for a comprehensive survey on this topic. } Among them, this study adopts the traditional approach developed by bresnahan1990entry and berry1992estimation that focuses only on the unique equilibrium outcomes. That is, we consider estimating the model based only on the information about “one-way links” in the network.
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
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 $\boldsymbol \gamma_0 = (\mathbf{A}_0^\top, \mathbf{B}_0^\top)^\top$. Then, the conditional probabilities of $\{y_{i,j}^{(1,0)} = 1\}$ and $\{y_{i,j}^{(0,1)} = 1\}$ are respectively given as follows:
Here, note that the equalities $y^{(1,0)}_{i,j} = y^{(0,1)}_{j,i}$ and $P^{(1,0)}_{i,j}(\theta, \boldsymbol \gamma) = P^{(0,1)}_{j,i}(\theta, \boldsymbol \gamma)$ hold. Thus, the likelihood function can be concentrated with respect to $(y_{i,j}^{(1,0)}, P_{i,j}^{(1,0)}(\theta, \boldsymbol \gamma))$; hereinafter, we omit the superscripts and denote $y_{i,j} = y_{i,j}^{(1,0)}$ and $P_{i,j}(\theta, \boldsymbol \gamma) = P_{i,j}^{(1,0)}(\theta, \boldsymbol \gamma)$ when there is no confusion. Then, the log-likelihood function can be written as
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.
As mentioned in Introduction section, 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:
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 $\{A_{0,i}\}$, 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 $\{B_{0,i}\}$, the agents can be grouped as $\mathcal{C}^B_0 \equiv \{\mathcal{C}^B_{0,1}, \ldots, \mathcal{C}^B_{0,K^B}\}$. When an individual is a member of the intersection $\mathcal{C}_{0,k}^A \bigcap \mathcal{C}_{0,l}^B$, his/her sender and receiver effects are equal to $a_{0,k}$ and $b_{0,l}$, respectively. This intersection set would correspond to one “community” in the community detection literature. 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., bonhomme2015grouped; okui2020heterogeneous). Then, following these studies, we treat $K^A$ and $K^B$ as known values and assume that $K^A, K^B \ge 2$.\footnote{ This assumption should be debatable. Rather than directly specifying the number of groups, in the literature on the BS algorithm for example, ke2016structure and lian2019homogeneity propose introducing an additional threshold parameter to detect the group structure. However, since the number of groups can vary only discretely with the threshold value, introducing the threshold parameter is essentially equivalent to selecting the number of groups. For the use of the BS algorithm, wang2020identifying formally shows that Bayesian Information Criterion (BIC)-type criterion can consistently select the correct number of groups as the sample size increases. } We discuss how to choose $K^A$ and $K^B$ in practice in Section (ref). 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 (ref). Thus, without loss of generality, we assume that $a_{0,1} = 0$ for location normalization.
Under this setup, the full ML estimator solves
where $\mathbf{a} \equiv (a_1, \ldots, a_{K^A})$ with $a_1 = 0$, $\mathbf{b} \equiv (b_1, \ldots, b_{K^B})$, $\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 (ref) is clearly a combinatorial (NP-hard) optimization problem. In the context of panel data models, several authors have proposed iterative k-means (like) algorithms to computationally obtain a (local) solution efficiently to the problems similar to (ref) (e.g., bonhomme2015grouped; liu2020identification). 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 decompose the maximization problem in (ref) into three steps. The first step is to estimate $\boldsymbol \gamma_0 = (\mathbf{A}_0^\top, \mathbf{B}_0^\top)^\top$ using the full ML estimator based on the log-likelihood function in (ref) 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 (e.g., bai1997estimating; ke2016structure; lian2019homogeneity; wang2020identifying). The final step is to solve (ref) with the group structure replaced by the estimated $\mathcal{C}_0^A$ and $\mathcal{C}_0^B$.
Before presenting the estimation procedure in detail, we discuss identification conditions for the parameters in model (ref). 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 practical reason for this is that our three-step estimator based on the BS algorithm requires preliminary consistent estimates of $\mathbf{A}_0$ and $\mathbf{B}_0$ in the estimation of the group structure. Note that we can always estimate “pseudo-true” group memberships based on the maximum likelihood principle even when some elements of $\mathbf{A}_0$ and $\mathbf{B}_0$ are not point-identified. However, they may not necessarily coincide with the true group memberships in general.\footnote{ A more formal investigation on this issue is left as a future work. For a related discussion, see bonhomme2015grouped. }
Here, again, we need some location normalization to identify $\mathbf{A}_0$ and $\mathbf{B}_0$. In the following, 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 (ref)(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 khan2018information, 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 (ref)(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 in link formation behavior 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 (ref) are standard in that they are satisfied in most of commonly used parametric models (such as logit and probit). In Assumption (ref)(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., yan2016asymptotics; yan2019statistical). The admissible parameter space for the correlation parameter $\rho$, $\mathcal{R}$, depends on the choice of the functional form of $H$, which is typically $\mathcal{R} = [-1,1]$. Assumption (ref) 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 (ref)--(ref) is that the one-way link-formation probabilities $\{P_{i,j}(\theta, \boldsymbol \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 (ref) and Section (ref) may be regarded as dense networks.\footnote{ graham2016homophily and jochmans2018semiparametric developed conditional likelihood methods that can be used to estimate the homophily parameters (i.e., $\beta_0$ in our context), even when the networks are sparse. However, their approach cannot be directly applied to our case because of the interdependence between $g_{i,j}$ and $g_{j,i}$. }
Assumption (ref) 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 (ref) does not hold. It should be noted that this assumption is not inconsistent with Assumption (ref) under the compactness of the parameter space (i.e., Assumption (ref)). 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., tamer2003incomplete; kline2015identification), 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 aradillas2019inference. That is, the model parameters can be point-identified either (i) if the distribution of unobserved payoff disturbances is fully known or (ii) if $\rho_0$ uniquely maximizes the concentrated log-likelihood function (ref). In the literature on network formation models, it is often assumed that $\epsilon_{i,j}$ and $\epsilon_{j,i}$ are independent (e.g., hoff2002latent; jochmans2018semiparametric; yan2019statistical). 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 (ref)(iii). We provide an explicit form of $\partial^2 \mathcal{L}_n^*(\rho)/(\partial \rho)^2$ in Appendix (ref). Even when neither of condition (i) nor (ii) is satisfied, it remains possible to partially identify the parameters, as in aradillas2019inference. For example, if $\mathcal{L}^*_n(\rho)$ has two peaks at $\rho_1$ and $\rho_2$, the resulting identified set is directly given by $\{(\widetilde \beta_0(\rho), \widetilde \alpha_0(\rho), \rho, \widetilde \boldsymbol \gamma_0(\rho)) : \rho \in\{\rho_1, \rho_2\}\}$.
The first step of the ML estimation aims to obtain consistent estimates of $\boldsymbol \gamma_0 = (\mathbf{A}_0^\top, \mathbf{B}_0^\top)^\top$. Let
where $\widehat \theta_n = (\widehat \beta_n^\top, \widehat \alpha_n, \widehat \rho_n)^\top$, and $\widehat \boldsymbol \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 (ref) with the main focus on $\widehat \boldsymbol \gamma_n$. Instead of introducing particular identification conditions, for generality, we directly assume that the true parameter $(\theta_0, \boldsymbol \gamma_0)$ is a unique maximizer of $\mathbb{E} \mathcal{L}_n(\theta, \boldsymbol \gamma)$.
We first establish several consistency results in the next theorem. In particular, we show that the individual specific effects can be uniformly consistently estimated.
Next, we derive the convergence rate of $\widehat \boldsymbol \gamma_n$. To this end, it is convenient to re-define $\widehat \theta_n$ and $\theta_0$ as
respectively, where $\widetilde \boldsymbol \gamma_n(\theta) \equiv \operatorname*{\arg\!\max}_{\boldsymbol \gamma \in \mathbb{C}_n} \mathcal{L}_n(\theta, \boldsymbol \gamma)$ and $\widetilde \boldsymbol \gamma_0(\theta) \equiv \operatorname*{\arg\!\max}_{\boldsymbol \gamma \in \mathbb{C}_n} \mathbb{E} \mathcal{L}_n(\theta, \boldsymbol \gamma)$ for any given $\theta \in \Theta$, assuming that they are well-defined. Further, we define $\boldsymbol \gamma_{-1} = (A_2, \ldots, A_n, B_1, \ldots, B_n)^\top$,
and $\mathcal{H}_{n, \theta\boldsymbol \gamma} (\theta, \boldsymbol \gamma) \equiv \mathcal{H}_{n, \boldsymbol \gamma\theta} (\theta, \boldsymbol \gamma)^\top$. The exact form of $\mathcal{H}_{n, \boldsymbol \gamma\boldsymbol \gamma} (\theta, \boldsymbol \gamma)$ can be found in (ref) in Appendix (ref). Finally, let
This $\mathcal{I}_{n,\theta\theta}(\theta, \boldsymbol \gamma)$ matrix serves as the Hessian matrix for the concentrated ML estimator $\widehat \theta_n$ (see, e.g., amemiya1985advanced). Now, we introduce the following assumptions.
Assumptions (ref) and (ref) 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 \boldsymbol \gamma_n$, as shown in the next theorem.
The max-norm convergence rate obtained in Theorem (ref) is consistent with the result of Theorem 3 in graham2017econometric and that of Theorem 3.1 in yan2019statistical.
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
In the following, we mainly describe the estimation of the group membership for the sender effects, $\mathcal{C}^A_0$. The exactly the same procedure described below can be used to estimate $\mathcal{C}^B_0$.
An 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)}\}$; namely,
Further, we define
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:
We find the first break point, say $\widehat t_1$, by
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.
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$. In other words, either one of the two converges to a sequence of constants as the sample size increases. Then, when we compute $\widehat S^A_{1,\widehat t_1}(\widehat t_1)$ and $\widehat S^A_{\widehat t_1 + 1, n}(n)$, and 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
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 given by $\{\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.
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 = \operatorname*{\arg\!\min}_{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
for $k = 1, \ldots, K^A$. 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 the true group structure $\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 (ref) is parallel to Assumption A2 in wang2020identifying. 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 (ref) is almost analogous to ke2016structure and wang2020identifying, for completeness, we provide it in Appendix (ref).
In the final step, we solve (ref) 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
where $\widehat{\mathcal{L}}_n(\delta) \equiv \mathcal{L}_n\left(\theta, \left\{ \sum_{k = 1}^{K^A} a_k \cdot \mathbf{1}\{ i \in \widehat{\mathcal{C}}^A_{n,k}\} \right\}, \left\{ \sum_{k = 1}^{K^B} b_k \cdot \mathbf{1}\{ i \in \widehat{\mathcal{C}}^B_{n,k}\} \right\} \right)$. Similarly, we define
where $\mathcal{L}_n(\delta) \equiv \mathcal{L}_n\left(\theta, \left\{ \sum_{k = 1}^{K^A} a_k \cdot \mathbf{1}\{ i \in \mathcal{C}^A_{0,k}\} \right\}, \left\{ \sum_{k = 1}^{K^B} b_k \cdot \mathbf{1}\{ i \in \mathcal{C}^B_{0,k}\} \right\} \right)$; 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 the standard parametric ML estimator, the estimator follows a normal distribution asymptotically at the parametric rate, and its asymptotic covariance matrix is given by the inverse Fisher Information matrix. Meanwhile, we have shown in Theorem (ref) 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, we can claim that the final ML estimator $\widehat \delta_n$ has asymptotically the same statistical performance as the oracle estimator $\widehat \delta_n^{\text{oracle}}$, and, thus, it is asymptotically fully efficient. 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)$.
In this section, we examine the finite sample performance of the three-step ML estimator. We consider the following data-generating process for the Monte Carlo experiments:
where $Z_{i,j,1} = |X_i - X_j|$ with $X_i \overset{i.i.d.}{\sim} \mathrm{Uniform}[-1,1]$, $Z_{i,j,2} \overset{i.i.d.}{\sim} N(0, 1)$, $(\epsilon_{i,j}, \epsilon_{j,i})$ is i.i.d. across dyads as the standard bivariate normal with correlation coefficient $\rho_0 = 0.6$, $(\beta_{0,1}, \beta_{0,2}, \alpha_0) = (-1.2, 1.6, 0.6)$, and $K^A = K^B = 3$. For the groupwise heterogeneity parameters, we consider $(a_{0,1},a_{0,2},a_{0,3}) = (0, r, 2r)$ and $ (b_{0,1},b_{0,2},b_{0,3}) = (-0.4 - r, -0.4, -0.4 + r)$ for $r \in \{0.4, 0.7, 1.0\}$. The smaller (larger) $r$ becomes, the more difficult (easier) the identification of the group structure. The group memberships are determined randomly while maintaining the equal size of each group. Exceptionally for observation $1$, $1 \in \mathcal{C}_{0,1}^A$ (so that $A_{0,1} = 0$) is fixed throughout the experiments. For each model setup, we consider two sample sizes: $n \in \{54,75\}$; thus, the size of each group is 18 in the former case and is 25 in the latter. The number of Monte Carlo repetitions is set to 500 for each single experiment. For the estimation of the group memberships, for comparison, we use both the standard BS method without repartitions and the repartitioned BS method.
We first report the simulation results of estimating the common parameters $(\alpha_0, \beta_{0,1}, \beta_{0,2}, \rho_0)$. Table (ref) presents the bias and RMSE (root mean squared error) for the following four estimators: the initial ML estimator given in (ref) (1st-step ML), the three-step ML estimator based on the BS method with no iterations (BS0) and that with two iterations (BS2), and the oracle estimator based on the true group membership (Oracle). For estimating $(\beta_{0,1}, \beta_{0,2})$, as expected, the 1st-step ML estimator is largely biased for all scenarios due to the incidental parameter problem. Although the three-step estimators (i.e., BS0 and BS2) also have some biases when $r = 0.4$ and $n = 54$, the biases disappear as either $r$ or $n$ increases. Thus, these biases are probably due to frequent misclassification of group memberships under small $r$ and $n$. In terms of RMSE, although we can observe a certain gap between the oracle estimator and the three-step estimators, the gaps can be reduced by increasing $r$ and $n$, which is consistent with our theory. Interestingly, even when using the 1st-step ML estimator, the strategic interaction effect and the error correlation parameter can be estimated with almost no bias.
The simulation results of estimating the group memberships are summarized in Table (ref). Here, we compare the performance of BS0 and BS2 in terms of the ratio of correct group classification. First of all, the results indicate that the repartitioned BS method (i.e., BS2) clearly outperforms the standard BS method without repartitions (i.e., BS0). As expected, as $r$ gets smaller, correctly predicting the group membership becomes significantly more difficult. If the gaps between the values of the group effects are sufficiently large and the sample size is not small, BS2 can attain almost 90% of correct classification.\footnote{ One might view that the results reported in Table (ref) are not particularly good for the BS algorithm. A main reason for that would be that our model is a bivariate binary response model, whereas most of the previous studies using the BS algorithm has focused on models with a continuous outcome. } We cannot observe any clear difference between the estimation of $\mathcal{C}^A_0$ and that of $\mathcal{C}^B_0$.
As an empirical application of our model and method, we analyze the network of international visa-free travels. The dependent variable of interest is $G_n = (g_{i,j})_{1 \le i,j \le n}$, where $g_{i,j} = 1$ if country $i$ allows the citizens in country $j$ to visit $i$ without visas, and $g_{i,j} = 0$ otherwise. Since the bilateral relationship about visa-free policy is expected to be complementary, this would fit into our model framework.
In this empirical study, we consider 57 countries selected mainly from Asia, the Middle East, the former USSR, and Oceania.\footnote{ The list of countries used in this empirical study is as follows: Armenia, Australia, Azerbaijan, Bahrain, Bangladesh, Belarus, Bhutan, Brunei, Cambodia, China, Cyprus, Estonia, Fiji, Georgia, Hong Kong, India, Indonesia, Iran, Iraq, Israel, Japan, Jordan, Kazakhstan, Kiribati, Kuwait, Kyrgyzstan, Laos, Latvia, Lebanon, Lithuania, Malaysia, Moldova, Mongolia, Myanmar, Nauru, Nepal, New Zealand, Oman, Pakistan, Papua New Guinea, Philippines, Qatar, Russia, Saudi Arabia, Singapore, South Korea, Sri Lanka, Tajikistan, Thailand, Tonga, Turkey, UAE, Ukraine, Uzbekistan, Vanuatu, Viet Nam, and Yemen. These countries are selected based on geographical proximity and ease of data collection. } The information about the visa policy of each country is taken from Henly and Partners: Passport Index 2020 (\url{https://www.henleypassportindex.com/passport}).\footnote{ Based on their definition, we categorize electronic travel authorization (eTA) and on-arrival visa as visa-free access. } The total number of dyads in this network is $57(57 - 1)/2 = 1596$. From Table (ref), which summarizes the distribution of the link connections, we can observe that the number of country pairs with one-way links is smaller than that with mutual links or no links. This would suggest the presence of complementarity in the network formation process. According to the above-mentioned passport index, Japan, Singapore, and South Korea are the top three countries among the 57 countries in terms of the number of all countries with visa-free access. For our restricted sample network, South Korea has the largest in-degree $\sum_i g_{i,\text{South Korea}} = 49$. For the out-degree, Nepal has the largest value $\sum_j g_{\text{Nepal},j} = 55$; that is, Nepal allows 55 countries (out of 57) to visit Nepal only with on-arrival visas. More detailed information can be found in Table (ref) in Appendix (ref).
The network for all the 57 countries is quite complicated and difficult to grasp the entire picture. As one illustration of our data, Figure (ref) presents the sub-network obtained by restricting the vertices to the Eastern and Southeastern Asian countries. The left panel in the figure shows the whole shape of this sub-network. (Note that the direction of the arrows in the figure is “not” the direction of visa-free access, but it represents that the target country is allowed to visit the country at the arrow's origin without visas.) The right panel shows the network created by leaving only one-way links from the left one; in other words, this is $(g_{i,j}\cdot (1 - g_{j,i}))_{i,j \in \text{Brunei}, \ldots, \text{Viet Nam}}$. From this figure, we can expect the existence of a certain level of degree heterogeneity. More specifically, Cambodia, for example, has five outgoing one-way links in this sub-network, suggesting that this country would have a larger sender effect $A$. In contrast, countries such as Japan and South Korea would exhibit a larger receiver effect $B$.
For estimating the network formation model, we consider five covariates; for their definitions, see Table (ref). The summary statistics of the covariates are provided in Table (ref) in Appendix (ref). With these variables, we consider the following payoff function:
where we assume that $(\epsilon_{i,j}, \epsilon_{j,i})$ have the standard bivariate normal distribution with correlation coefficient $\rho_0$. To estimate our network formation model with grouped degree heterogeneity, we first need to determine the number of groups for the sender effects $\{A_{0,i}\}$ and that for the receiver effects $\{B_{0,i}\}$, $K^A$ and $K^B$, respectively. Then, following ke2016structure and wang2020identifying, the optimal $(K^A, K^B)$ is selected as the minimizer of the BIC criterion: $-2\widehat{\mathcal{L}}_n(\widehat \delta_n) + (6 + K^A + K^B)\ln (1596)$. Then, as a result of searching over the models with $(K^A, K^B) \in \{2, \ldots, 7\}^2$, we find that the model with $(K^A, K^B) = (7, 6)$ achieves the smallest BIC (see Table (ref) in Appendix (ref) for more detailed information), and this is the model reported here. For comparison, we estimate not only our proposed model, which we call the grouped heterogeneity model, but also a dyadic bivariate probit model without strategic interaction and degree heterogeneities as a benchmark. For the estimation of the group memberships, we employ the repartitioning method with two iterations (i.e., BS2 in the previous section).
The estimation results are summarized in Table (ref). First of all, as expected, our proposed model suggests that there is a significant strategic complementarity in the network formation behavior. We can also find a certain level of degree heterogeneity in terms of both the sender and the receiver effects. Comparing the grouped heterogeneity model and the benchmark model, the log-likelihood value for the former is apparently significantly larger than that for the latter. This large difference in the degree of model fitting also demonstrates the significance of strategic effect and unobserved heterogeneity (note however that these models are not nested). For specific parameter estimates, we can observe several non-negligible differences between the two models. For example, the effect of the export amount is predicted to be positive in the grouped heterogeneity model, whereas the benchmark model predicts a significantly negative impact. The error correlation parameter is not significantly different from zero in our model but is weakly positively significant in the benchmark model. This result would be understandable since the benchmark model can account for the interdependence of the links only through the error correlation. Additionally, there are several interesting findings. For both models, if countries $i$ and $j$ are located in the same region, they become more likely to allow visa-free access, as expected. Not only in terms of geographical proximity, but we can also observe significant homophily in terms of the political system.
For the estimation results of country-specific effects, we report the estimated group memberships in Table (ref). As expected from the above discussion, countries such as Cambodia are indeed classified into the highest group (i.e., Group 7) in terms of the sender effect. The other two countries that have Group-7 sender effect are Nepal and Sri Lanka. For the receiver effect, these two countries are classified as Group 2, and Cambodia is in Group 3. Overall, interestingly, there seems to be a weak negative correlation between the sender effects and the receiver effects. As expected from the above discussion, Japan and South Korea indeed belong to the group with the highest receiver effect (i.e., Group 6). The magnitudes of the receiver effects seem to roughly correlate with the size of the countries' economies (with some exceptions, such as China, India, and Russia).
This paper proposed a network formation model with pairwise strategic interaction and grouped degree heterogeneity. Assuming some parametric form for the error distribution, we proved that the model parameters can be identified under the availability of agent-specific covariates that have large supports and also have variations across all potential partners. For estimating the model, based on the same idea as in bresnahan1990entry and berry1992estimation, we proposed the three-step ML procedure: in the first-step, the model is estimated without considering the group structure; subsequently, we estimate the group memberships using the BS algorithm given the estimates for the heterogeneity parameters obtained in the first step; and, finally, based on the estimated group memberships, we re-estimate the model. Under certain regularity conditions, we showed that the proposed estimator is asymptotically unbiased and distributed as normal at the parametric rate. The results of the Monte Carlo simulations show that our estimator performs reasonably well in finite samples. An empirical application to international visa-free travel networks indicates the usefulness of the proposed model.
Several limitations and extensions are as follows. First, our approach can be used only in pairwise network formation games with no network externalities to/from the rest of the links, and this limits the empirical applicability. Therefore, it would be worthwhile to extend our results to network formation models with general network externalities involving more than two agents. However, we conjecture that we would resort to partial identification to achieve this. Second, our approach requires that the degree heterogeneity parameters have discrete support, although, in reality, it is possible that they are continuous. To address this issue, it is of interest to modify our model in a similar manner to bonhomme2017discretizing and investigate the three-step ML estimator in which $K^A$ and $K^B$ grows slowly to infinity. Third, as our model is a dyadic binary game model, where a pairwise network formation model is its special case, we can consider its ordered-response game version as a natural extension. For example, we might be interested in analyzing bilateral military relations: non-alliance, quasi-alliance, or alliance. We expect that such extension can be relatively easily achieved by adopting the ML estimator discussed in aradillas2019inference. Finally, related to the empirical application in this study, we might be interested in investigating the causal effect of visa policies between two countries on the flows of tourists between them; this is a dyadic treatment evaluation problem when the treatment variable is determined strategically. To deal with such situations, combining the results of this study and the marginal treatment effect framework developed in hoshino2020treatment would be beneficial. We leave these topics for future research.