An Econometric Model of Network Formation

(1)

Bachelor Thesis - Specialisation Econometrics

An Econometric Model of Network Formation

Anna-Liisa Distefano (11371722) supervised by Sanna Stephan

(2)

Abstract

This paper is a comparison study on link functions within the binary estimation program of the homophily parameter of an undirected dyadic link formation network model as proposed by Graham (2017). Estimation of the homophily parameter is conducted by means of a heuristically derived estimator, namely the tetrad logit estimator (or TLE), which main property is that it can estimate the homophily coefficient within an undirected network, independent of unrestricted agent level degree heterogeneity. This paper contributes to research by introducing the possibility of adaptation of the TLE to networks which underlying error term is distributed non-logistically, and to open the doors for advanced derivation of an analogous tetrad probit estimator.

(3)

1 Introduction

The Information Age and the associated Digital Revolution are impacting human reality in terms of the way we communicate, commute, form relationship and understand complex natural events. As we interact, collaborate and exchange ideas faster and more frequently than ever, it becomes increasingly unrealistic to picture individuals as independent entities. Instead, we can view agents as part of their "social network". As a result, the science of Network Formation is gaining attention from a broad variety of fields. Initially, it was predominantly developed within the fields of psychology and social sciences. In recent years however, increasing interest has been shown by the scientific communities of mathematicians, statisticians and physicians, who take a deeper analytical and numerical approach to the analysis of network formation (Newman, Watts, &Strogatz, 2002). Specifically, as the field of economic research has advanced to favor mathematical modeling, the study of socio-economic networks has become a field of interest for application of the network formation models put forward by mathematicians. Indeed, there is an increasing demand for models that explain the formation of various socio-economic networks we observe such as online and offline social networks.

Currently, network formation models are generally classified into two branches, namely the dynamic models, and the agent-based models. The latter involve a statistical analysis in which the network dynamics are often simulated based on theoretical assumptions on the topography of the network after observation of a certain network feature (such as a degree distribution). The aim of this type of models is to put forward a data generating process that can mimic this feature. The former take a micro-founded approach, often based on the microeconomic notions of game-theory. Attempts to merge the two approaches have been done in the hope to build a more general framework (Mele, 2017; Borshchev&Filippov, 2004), but the research on the comparison and combination of the two models is poor (Figueredo&Aickelin, 2011). For this reason, the econometric model proposed by Graham (2017) is of particular interest. Some of the fundamental assumptions of his research stem from microeconomic theory, while his methodology suggests an approach close to an agent-based one. Grahams research can be interpreted as an attempt to elevate a purely dynamic or purely agent-based model to a more complex form, one that merges important characteristic of a variety of network formation models. This is important because of the applicability of network theory across a broad variety fields.

(5)

(2017). The model takes into account parameters for homophily and for individual degree heterogeneity, due to the intuitive and literary evidence (as discussed later) that agents tend to link according to these properties. Further, the paper introduces and studies two estimators for the homophily parameter, namely the Tetrad Logit Estimator and the Joint Maximum Likelihood Estimator, under the assumption of the random error distribution being an i.i.d. logistic one. This assumption about the anatomy of the network is crucial because Graham uses it to derive the estimation program of the TLE. A question arises then if this TLE is robust to the random variable having a different distribution link, for example a probit one. The purpose of this paper is to test if the TLE is usable in networks which idiosyncratic component does not appear to be distributed according to a logistic distribution. To achieve this, the results of Graham (2017) are directly replicated, and compared with the simulation of the TLE under a different assumption regarding the link function. Lastly, a test is performed for the difference between the simulated estimates. The paper is structured as follows: in the first section the main theories about network formation are discussed along with the model proposed by Graham (2017), and background theory about the types of link functions. Secondly, the methodology used for the simulation of the Montecarlo experiments and the estimation procedures are covered. Lastly, the results are displayed and analysed, and a conclusion is drawn with respect to the research question, explicitly stated below.

I aim at contributing to the fast-growing field of Network Formation with a study on the comparison of link functions, specifically within the context of a logit estimator within a class of network formation models, by answering the question: Is the tetrad logit estimator as derived by Graham (2017) applicable to networks in which the shock term is not assumed to be a logistically distributed?

2 Theoretical Framework

2.1 Random Graph Models and Strategic Network Formation Models

A random graph is a mathematical structure which objects link according to a probability distribution. One of the first contributions to the literature of network simulation by random graph modelling is given by Erdos and Rényi (1960). Their research proposes a probability model for the structure of a random graph and a study of the evolution

(6)

of it as the number of edges is increased. This model however fails to replicate clustering1_{, hubs}2_{and the formation}

of triads, which are nevertheless observed in reality. Newman, Watts and Strogatz (2002) brought the model forward by developing an exactly solvable network model with a clustering coefficient and varying degree distributions. This approach is still somewhat simplistic as it poses a series of obstacles to the generalisation and prediction of models, one of them being the inability to model (continuous) value-edge networks (Desmarais&Cranmer, 2012).

For this reason, the Generalised Exponential Random Graph Model (GERM) was introduced within the science of network modelling. The GERGM is characterised by a probability model of an observed network that conditions on relevant observed structures that affect the topography of a network. The main advantages of the GERGM are its capability to model networks with continuous valued edges (bounded or unbounded) and their ability to include both endogenous and exogenous factors as explanators for prediction. Furthermore, the network statistics can be specified such that they capture important characteristics such as homophily3, degree heterogeneity4, transitivity5 and reciprocity6. (Desmarais&Cranmer, 2012). Examples of econometric network formation models that make use of GERM are given by Gualdani (2018), and Lu, Jerath and Singh (2013).

An alternative way to model networks is by recognising a network as a set of individuals whose intent is to maximise their utility, hence the choices of which can be described by a utility function. Models of this nature are called strategic network formation models and are based on the assumption that individuals form links with others based on their own utility maximisation. In this context, a network is then viewed as the equilibrium outcome of strategic interactions within sets of these utility maximising individuals (Sheng, 2012). For example, the formation of a network can be described as a simultaneous-move game with incomplete information (Ridder&Sheng, 2017; Leung, 2015), or as a sequential-move game (Mele, 2017). More specifically, the game-theoretic model proposed by Ridder and Sheng (2017) allows for non-separable utility7 _{and it proposes a two-step estimation procedure that is}

claimed to be applicable to both directed and undirected networks, as well as discrete and continuous variables. The latter estimation procedure was imported from Leung (2015), which inspired Ridder and Sheng (2017) also for the property for which linking decisions are allowed to depend on both agent attributes and network topography.

1_{grouping of a set of nodes so that the nodes in the same cluster resemble each other more than to those in other clusters} 2_{nodes with a greatly higher than average number of links}

3_{homophily is the theory that similar nodes may be more likely to attach to each other than dissimilar ones} 4_{the degree of a node is the number of connections it has to other nodes}

5_{if there is a link from i to j , and also from j to h, then there is also a link from i to h. eg. friends of friends are friends} 6_{measure of the likelihood of vertices in a directed network to be mutually linked}

(7)

2.2 The Graham Model

The network formation model as proposed by Graham (2017) is an example of an econometric model which bases some of its assumptions on microeconomic theory. The model is not exactly a strategic network formation model nor a RGM or a GERGM, which makes it a network model of peculiar nature and complexity. It is a model of undirected dyadic network formation, the agents of which are allowed to link by homophily and degree heterogeneity. The aim of the model is to introduce and compare two estimators for the homophily parameter, namely the Tetrad Logit Estimator (TLE) and the Joint Maximum Likelihood Estimator (JMLE). The choice to focus on homophily stems from the fact that the social interactions between individuals have been shown to heavily depend on the similarity between themselves and the targeted potential connection (McPherson, Smith-Lovin, &Cook, 2001). The estimators are developed and studied on the grounds of Monte-Carlo experiments that differ in network design, and are analysed within the asymptotic theory of large networks.

The network is simulated according to the following rule. The interaction of agents i and j yields to a total surplus, which is defined by a sum of observed dyad attributes, unobserved agent-level attributes, and a random component. Graham states on the basis of utility theory that the agents i and j will only form a link when the total surplus from linking is greater than 0. Based on this, the network is simulated according to the adjacency matrix8being defined as follows:

Dij =1(XiXjβ0+ Ai+ Aj− Uij≥ 0) (1)

With Dij being an element the adjacency matrix, and where for this Monte Carlo experiment, the homophily

parameter β0 is set to 1 and the agent level attributes X are set to take on the values of Xi = {−1, 1} . The

expression XiXj is the 1x1 dyad level element that captures the agent level attributes of the network. Its form

implies a strong propensity to homophilic linking because the product can be either −1 (when Xi6= Xj) or 1 when

the agent attributes are equal. Further, the components Ai and Aj are fixed effects that determine the unobserved

agent-level attribute of degree heterogeneity. Finally, the last component Uij is the random component, assumed

to be i.i.d. across dyads.

Having defined the matrix D as above, Graham (2017) imposes five baseline assumptions, the first one being that the conditional likelihood of the network D is:

(8)

P r(Dij = d|X, A0) = ₁ 1 + exp{XiXjβ0+ Ai0+ Aj0} 1−d _exp{X iXjβ0+ Ai0+ Aj0 1 + exp{XiXjβ0+ Ai0+ Aj0 d (2)

By the definition of a logit likelihood, this assumption is an implication that the random component Uij follows

a logistic distribution. Furthermore, U is assumed to be symmetric such that Uij = Uji. Given this simulation

structure for the network D and the conditional likelihood of it, the estimator of the homophily parameter β0 is

obtained by maximising the conditional log-likelihood derived from (2) and (1).

The conclusions of his findings are that assuming a logistic idiosyncratic component Uij, all networks have the

property of transitivity and exhibit a single giant component9. The estimated ˆβT Lare found to remain well-defined

within all of the Monte-Carlo simulations, while the alternative estimator ˆβJ M L is found to fail to exist in certain

designs. This suggests that the ˆβT L remains insensitive to the topography of the network it is subjected to, while

usage of the ˆβJ M L estimation heavily depends on the density of the network. The decision to focus on the ˆβT L

comes in view of these results.

2.3 Types of Binary Choice Model Link Functions

A generalised linear model for binary response data is defined as: P r(y = 1|x) = g−1(x0β), with the variable y being the response variable that takes on the values of either 1 or 0, depending on wether a specified event occurs or not, β being the regression coefficients and g being the link function (Van Horn, 2015).

Econometric literature has proposed a wide range of options of link functions for binary choice models (Koenker&Yoon, 2009; Bazan, Romeo, &Rodrigues, 2014; Gilchrist&Green, 1983; Guerrero&Johnson, 1982), and has criticised the widespread idea that links are interchangeable, bringing forward the idea that there exists a link function that best describes the data (Koenker&Yoon, 2009). In face of this proposition, a test to examine the adequacy of a hypothesised link was developed by Pregibon (1980). He constructed a procedure based on the assumption of the existence of a true but unknown link function g∗, being fitted by another g0. His hypotheses are formulated as

follows:

H0: g0(µ) = g(µ; α0, δ0) and Ha: g∗(µ) = g(µ; α∗, δ∗) (3)

with the first equation being the hypothesised link, and the second being the assumed to be true link function.

(9)

Pregibon (1980) concludes that the choice of the link function is one that should indeed be made with statistical consideration. In this regard, four link functions are discussed, namely the well-known (1) logit and (2) probit links, aswell as the (3) cloglog and (4) logloglink.

Firstly, the logistic model is introduced because of its notoriety and widespread applications. The reason is that the logit link yields to an exactly solvable cumulative distribution function and likelihood (Heij et al., 2004). Its probability distribution is defined according to the following transformation of the dependent variable yi:

logit: λ(yi) =

eyi

1 + eyi (4)

Secondly, the other widely-used link function is given by the probit link. The probit transformation is given by the probability density function of the standard normal distribution10:

probit: φ(yi) = 1 √ 2πe −1 2y 2 i ₍₅₎

In contrast to the logit link, the CDF of the latter cannot be calculated exactly and must be computed numerically by approximation of the integral over the pdf. Unless the tails of the distribution are significantly heavy, the integration algorithms for the approximation of the CDF lead to very similar results to the exact results of the logit model (Heij et al., 2004). However, the normal and logistic densities lead to different marginal effects, such that approximately βLGT = 1.6βP BT. Lastly, the loglog link and the complementary-log link are described by the

following pdf: loglog: g1(yi) = e−e (α+βx) (6) cloglog: g2(yi) = 1 − e−e (α+βx) (7)

These links are often used for extreme asymmetric distributions. For example, the cloglog link arises when y describes a nonzero count that can be modelled with a Poisson distribution (Van Horn, 2015). On the other hand, the loglog link is rarely used in practice because the transform is inappropriate for pi < 1₂ (ChambersCox, 1967).

Furthermore, the above links can also be defined as functions of pi, the probability of the occurrence of an event

yi = 1, rather than the event yi= 0 (Koenker&Yoon, 2009):

logit: g(pi) = log(pi/(1 − pi))

probit: g(pi) = Φ−1(pi)

(10)

cloglog: g(pi) = log(−log(1 − pi))

loglog: g(pi) = log(log(pi))

All the links described above are parametric link functions, that is they allow for estimation of the parameters of a given model based on a probability distribution. It is important to note however that estimation by nonparametric identification also exists, for example in the index model of link formation as proposed by Gao (2018). The study of these links is beyond the scope of this paper.

3 Methodology

3.1 Simulation Frameworks

For estimation of the TLE, networks are simulated as follows. The number of nodes N describes the size of the network, and is taken at N = 40 so that the number of dyads n is N₂ = 780 and the number of tetrads N₄ = 91390. For comparison, six different types of networks are generated: three that are symmetrically distributed with degree heterogeneity uncorrelated with agent level attributes, and three that are right skewed with correlation between degree heterogeneity and the agent level attributes. For each type of network 100 Monte Carlo replication studies are performed. The networks are simulated based on the adjacency matrix in (1) and the unobserved individual level degree-heterogeneity A is generated as follows:

Ai= αL1(Xi = −1) + αH1(Xi= 1) + Vi (8)

where Vi|Xi ∼ {Beta(λ0, λ1) − λ0/(λ0+ λ1)} and αL ≤ αH. Ai takes continuous values on the interval [αL − λ0

λ0+λ1, αH−

λ1

λ0+λ1] and has conditional expectations E[Ai|Xi= −1] = αL and E[Ai|Xi= 1] = αH .

The parameters αL, αH, λ0, λ1are used to calibrate the topography of the networks by setting their values in six

combinations of the form [αL, αH, λ0, λ1], with αL∈ {−12, −1, −2, − 2 3, − 7 6, − 13 6}, αH∈ {−12, −1, −2, − 1 6, − 2 3, − 5 6},

λ0∈ {1, 1, 1,1₄,1₄,1₄} and λ0∈ {1, 1, 1,3₄,3₄,3₄} so that the simulation designs can be classified into the two distinct

families mentioned above. The distribution is assumed to be a centered Beta one because of its flexibility and its opportunity to control skewness and kurtosis. Specifically, the λ values control the shape of the distribution so that by definition of the skewness of a Beta distribution, the network distribution is symmetric when λ0= λ1 and

(11)

contribution of the agent level attributes Xi, so that the first two terms in (8) can never contribute simultaneously

since Xi can either be 1 or −1. Furthermore, when αL = αH, the sum of the first two terms is the same for any

value of Xi, which implies that there is no correlation between agent level attributes and degree heterogeneity.

When αL 6= αH, the sum of the first two terms depends on the value of Xi, which implies correlation (Graham,

2017).

3.2 Estimation of the TLE

It is natural to attempt the estimation of β0directly from optimising the conditional likelihood of the network. This

likelihood function is assumed to have the following form for a given cumulative distribution function F, observed agents characteristics Xi, Xj and unobserved agent characteristics Ai, Aj:

P r(Dij = d|X, A) = F (XiXjβ0+ Ai+ Aj)d(1 − F (XiXjβ0+ Ai+ Aj))(1−d) (9)

which form depends on the choice of link function.

However, estimation of β0based on (9) leads to an optimisation problem that is not exactly solvable due to the

identification problem that occurs because a different unknown degree is randomly assigned to each agent. This arbitrary nature of the degree sequences vector A results in a model with a large number of unknown parameters that are hardly identified. Consequently, the empirical estimation procedure of β0 is based on probabilities of

the observable configurations of the tetrad subgraphs that are independent of A, rather than on the conditional likelihood as defined in (9). Graham (2017) derives this methodology by approaching a heuristic implementation for TLE, based on the relative probability of observing a dyadic link within a tetrad subgraph Sij,kl, conditional

on observing one in another permutation of it. A graphical representation is provided in the figure below: when the dyadic link does not appear in neither of the forms illustrated below, the value of Sij,kl= 0.

To understand how Sij,kl is defined, recall the definition of Dij in (1). The tetrad subgraph Sij,kl = 1 when

a dyadic link is formed between nodes i, j and k, l but not i, k and j, l, which implies that Sij,kl = 1 when

Dij = 1 ∩ Dkl = 1 ∩ Dik = 0 ∩ Djl = 0. Analogously, Sij,kl = −1 when a link is formed between nodes i, k and

j, l but not i, j and k, l, which means that Sij,kl = −1 when Dij = 0 ∩ Dkl = 0 ∩ Dik = 1 ∩ Djl = 1. The tetrad

subgraph configuration element Sij,kl is then defined as follows:

(12)

with Dij being the element of the matrix as defined in (1), where Dij = 1 if the contribution of i and j forming

a link is greater or equal to the contribution of the error term for edge ij, so when XiXjβ0+ Ai+ Aj ≥ Uij, and

similarly Dij= 0 when XiXjβ0+ Ai+ Aj< Uij.

This definition allows for the introduction of a supporting variable Wij,kl that describes observable utility of edge

formation independently of degree heterogeneity vector A. For Sij,kl = 1 it is true that XiXjβ0≥ Uij− Ai− Aj,

XkXlβ0 ≥ Ukl− Ak − Al, and XiXkβ0 < Uik− Ai− Ak and XjXlβ0 < Ujl− Aj − Al. For Sij,kl = −1 it is

true that XiXkβ0 ≥ Uik− Ai− Ak, XjXlβ0 ≥ Ujl− Aj− Al, and XiXjβ0 < Uij − Ai − Aj and XkXlβ0 <

Ukl − Ak − Al. Based on this, let the 1x1 element of agent attributes that does not depend on Ai and Aj be

defined as Wij,kl = XiXj + XkXl− (XiXk + XjXl), so that Wij,klβ0 > 0 means that (XiXj + XkXl)β0 >

(XiXk+ XjXl)β0. The condition Wij,klβ0> 0 then identifies the case in which the observable surplus of Sij,kl= 1

is larger than that of Sij,kl= −1, which probability can be simplified to a sum of the error terms by substitution of

the aforementioned conditions as follows. Wij,klβ0= (XiXj+ XkXl− (XiXk+ XjXl))β0≥ 0 implies by comparison

that (Uij− Ai− Aj) + (Ukl− Ak− Al) − (Uik− Ai− Ak) − (Ujl− Aj− Al) ≥ 0 which simplifies to an expression

that is independent of degree heterogeneity A shown below:

P r(β0Wij,kl> 0) = P r(β0(XiXj+ XkXl− (XiXk+ XjXl)) > 0) = P r(β0(Uij+ Ukl− Uik− Ujl) ≥ 0) (11)

If U is logistically distributed, this probability is a sum of logistic probabilities which is also logistically distributed. If U follows a normal distribution, this probability is normally distributed because a sum of gaussian densities is gaussian as well (Bain&Engelhardt, 1987). By symmetry of the logistic distribution and the normal distribution, and based on the latter definitions, the ˆβT Lcan then be derived as the solution to a new log-likelihood maximisation

(13)

problem based on the conditional probability independent of degree heterogeneity as:

P r(Sij,kl= 1|X, A, Sij,kl∈ {−1, 1}) = F (Wij,klβ0) (12)

which leads to the conditional log-likelihood for the distinct configurations Sij,lk:

lij,kl(β0) = log{(F (Sij,klWij,klβ0))|Sij,kl|)} = |Sij,kl| log{F (Sij,klWij,klβ0)} (13)

with F (.) being the CDF function of the distribution in question. From this log-likelihood, Graham (2017) develops the estimation procedure of the TLE around the theory of U-statistics. A U-Statistic is is defined as a statistic that averages across all permutations of the variable arguments of a symmetric function which expectation results in an unbiased estimable parameter (Ferguson&Thomas, 2005). This is of important relevance in this case because of the fact that a different permutation of the tetrad {ij, lk} leads to a different value of the conditional likelihood (13). The criterion function is then defined as the average of the contributions of the likelihoods for all the 4! permutations of the nodes across each of the N₄ possible tetrads. However, Graham (2017) finds that only 3 of those 4! permutations of the indices of (13) are relevant, namely {ij, kl}, {ij, lk} and {ik, lj}, so that the ultimate criterion function is given by:

LN(β) = N 4 −1₁ 3 X i<j<k<l

[lij,kl(β) + lij,lk(β) + lik,lj(β)] (14)

This definition is important because it suggests that at least one term in (14) will be nonzero, and that consequently the tetrad corresponding to that permutation will have a nontrivial contribution to the TLE (Graham, 2017).

The TLE is found solving the f.o.c of LN(β) by the means of a simplified 4 step logit estimation program

based on a fit of the simulations of all Sij,kl6= 0 onto the corresponding simulated Wij,kl. The first step consist of

calculating S and W for all N₄ for the three relevant permutations of the indices mentioned above. Secondly, a dataset is created by concatenating these estimates so to obtain a 3 N₄ × 2 matrix, where the elements of the first column are the values S and those of the second are the values W . Thirdly, all rows i with S = 0 for i = 1, 2, ..., 3 N₄ are dropped so that only the relevant cases Si = 1 and Si= −1 are kept, so that the dimensions of the data matrix

become R × 2, where R = 3 N₄ − q and q is the number of rows in which S = 0. Lastly, the elements Si= −1 are

set to 0, so that the first column can be modelled as a dependent random variable ˜S ∈ {0, 1} ∼ Binomial. The assumption of the binomial distribution is based on the notion that if a random variable has a binary outcome of either success or not, it can be described as such (Bain&Engelhardt, 1987). The TLE is found by estimation of the

(14)

regression:

˜

Sj = βWj (15)

for all j = 1, 2, ..., R. The constant is not included in the fit because inclusion of an intercept in the model equation is nonsensical, and the regression of ˜S onto W can be performed with either logit or probit link. The estimate of β0by tetrad logit estimation is the estimated coefficient on W .

3.3 Testing

The aim of this research paper is to test the robustness of the TLE to changes in the distribution of the error term of the network simulation expression. This can be achieved by setting up a hypothesis testing procedure on link functions inspired by the idea proposed by Pregibon (1980). In face of this, two DGPs are considered, in combination with two different estimation programs. The first part of the analysis consists of the estimation of β0,

with Uij set to follow the standard logistic distribution, and the link being logit, so that β0 is estimated by the

correct estimation program. The second part consists of estimating β0 with Uij set to follow the standard normal

distribution, and of applying the "wrong" logit estimation program. The third part consists of estimating β0 with

Uij set to follow the standard normal distribution, using the adapted estimation program in which the link is set to

be probit, the assumed to be correct estimation program. The framework of the testing procedure can be structured as follows. Define βLas the estimator which network has a logistic error term Uij and is estimated under the correct

logit estimation program; define βN as the estimator of a network which Uij follows the normal distribution and

is estimated by using the logit estimation program; lastly, define βP as the estimator which network has gaussian

error and is estimated with the correct probit estimation program. Since the true value of β0used to generate the

data equals β0= 1, the tests can be formulated as simple two sided tests on the single homophily parameter:

H0: βL= 1 Ha: βL6= 1 (16)

H0: βN = 1 Ha: βN 6= 1 (17)

H0: βP = 1 Ha: βP 6= 1 (18)

The hypotheses in (16) can be interpreted as a test of the validity of the TLE, meaning that if βL is found to not

be significantly different from its true value, then it can be concluded that the TLE derived by Graham (2017) can be replicated correctly and is valid. The hypotheses in (17) on the other hand are to be interpreted as a test for

(15)

the robustness of the logit estimator to changes in the distribution of the shock term. Hypotheses (18) investigate whether − given a different error function − readily adjusting the last step of the estimation program by solely changing the link function is sufficient for correct identification. If adjustment of the estimation program by change of link function is sufficient, the estimator βP should generate correct rejection frequencies when tested against the

truth. Conclusions of the tests (17) and (18) are used to asses whether or not a specific assumption on the link function, together with a specific estimation program leads to correct parameter identification.

To test these hypotheses, the large samples properties of this estimator are considered. Graham (2017) finds that asymptotically the β follows a normal distribution with mean β0 its true value, and a variance of the form

ˆ H−1J0Hˆ−1, so that approximately: ˆ β ∼ N (β0, 36 n ˆ H−1J0Hˆ−1) (19)

with n = N₂ being the number of dyads. Graham(2017) proves that ˆH can be approximated by extracting the Hessian from a logit estimation program. The J0can be approximated as follows. Define ˆs as the vector of length

n that contains the average of the N −2₂ gradients of the tetrads that each of the n dyads belong to. The indices are set to permute under the conditions S = {i, j} ∩ {k, l} = ∅, i < j, k < l, and the components are subscripted by the indices {i, j} to denote the dyads:

ˆ sij = 1 (n − 2(N − 1) + 1) X S ∇β 1

3(lij,kl(β) + lij,lk(β) + lik,lj(β)) (20) where for example ∇βlij,kl = |Sij,kl|Sij,kl

Wij,kl

(1+eSij,klWij,klβ) for the first of the three permutations of the indices in the case of

the logistic distribution. Finally, the J0 is defined by taking the average of the square of the n components of the

vector in (20), evaluated at the estimated value of β0:

J0= 1 n X i<j ˆ sij( ˆβ)ˆsij( ˆβ)0 (21)

With these specifications, a t-test for each of the parameter differences can be performed. If in the sample of 100 Monte Carlo simulations, the difference between the estimates turns out to be different from 0 for more than the significance level of α = 5%, then it can be concluded that the H0 is rejected in favour for the Ha, suggesting

(16)

4 Results and Analysis

4.1 Network Statistics

The primary analysis is performed on networks which size is set to be N = 40, with 100 Monte Carlo experiments per simulation. The Monte Carlo results and test results for smaller networks of size N = 20 and N = 30 are displayed in Appendix A 11_{. Table [1] shows the network statistics for networks under a standard logistic distribution of}

the U term, and Table [2] shows the same results for U following a standard normal distribution. Networks {A1, A2, A3} describe networks with symmetric uncorrelated heterogeneity, and {B1, B2, B3} describe networks with right-skewed correlated heterogeneity 12.

Table 1: Network Statistics Logistic Distribution

A1 A2 A3 B1 B2 B3 Density 0.3048 0.1443 0.0044 0.3352 0.1713 0.0143 Av. Degree 12.2690 6.3225 1.1395 13.3925 7.3215 1.5050 Std. Degree 3.3791 2.5196 1.0768 4.0462 3.3329 1.3614 Clustering C. 0.4012 0.2360 0.0499 0.4368 0.2812 0.0836 U ∼ Logistic(0,1)

The density of an undirected network is given by _{N (N −3)+2}E−N +1 where E is the number of edges in the networks and N the number of nodes. The average degree is simply the average of the degrees of the nodes in a network, and the standard deviation of the degree is the standard deviation of the degrees of the nodes within the network. The clustering coefficient is defined as the ratio of the number of existing links connecting a certain node’s neighbors to each other over the maximum possible number of such links. This can be formulated as Ci = _b_i_(b2a_i₋₁₎i where bi is

the number of neighbours of node i, and ai is the number of connections between the neighbours (Bondy & Murty,

1976). From the tables it can be seen that networks A1 and B1 always have the highest value for density, whereas networks of kind A3 and B3 exhibit the lowest. This pattern is also observed in the average degree of the networks,

11_{Analysis on bigger size networks are not available due to technical constraints and limited processing capacity of the devices used}

for this study

12_{A1: {−}1 2, − 1 2, 1, 1}, A2: {−1, −1, 1, 1}, A3: {−2, −2, 1, 1}; B1: {− 2 3, − 1 6, 1 4, 3 4}, B2: {− 7 6, − 2 3, 1 4, 3 4}, B3: {− 13 6, − 5 3, 1 4, 3 4}

(17)

and clustering coefficients. From this it can be deduced that within the subgroups A and B, across the 3 different network topographies the first one is the densest while the last one is the sparsest, as intentionally calibrated by the chosen values of α and λ. Also, the standard deviation of the degrees is observed to become smaller for sparser networks, suggesting that in the sparser networks nodes tend to have less variation in the degree level.

Table 2: Network Statistics Normal Distribution

A1 A2 A3 B1 B2 B3 Density 0.2538 0.0689 0.0000 0.2823 0.1043 0.0000 Av. Degree 10.3780 3.5275 0.0480 11.4345 4.8380 0.1920 Std. Degree 3.5318 2.1720 0.1687 4.8025 3.4873 0.4468 Clustering C. 0.5039 0.2319 NA 0.5259 0.3392 NA U ∼ N ormal(0, 1)

Moreover, it is observed that for the case of the normal distribution, the density of the sparse networks reaches 0 and the clustering coefficients are numerically not found. This suggests that the estimation of βN and βP for

networks A3 and B3 might turn out to be problematic. Comparing the statistics of Table [2] to those of Table [1], it can be said that for the case of the normal distribution, all values excluding the standard deviations of the degree are lower than for the logistic distribution. This suggests that in general, if U is normally distributed and the network is described according to (1), it is likely to be sparser and with higher difference in degree than when U is set to be logistically distributed.

Furthermore, a comparison of Table [1] to the analogous table in Graham (2017) infers that even if the size of the network differs by 60 nodes (N = 40 versus N = 100), the densities of the networks and the clustering coefficients are comparable. This suggests that the network simulation procedure proposed in Graham (2017) has the potential to lead to reliable results in smaller networks than the ones studied in his paper. The other statistics however, turn out to be very different. This is no surprise though, as by definition the average degree and the standard deviation of the degree directly depend on quantity of nodes present in the network. Based on the above results, it can be expected that estimation of βL will lead to better results compared to those of βN and βP.

(18)

4.2 Estimation Outputs

Table [3] shows the Monte Carlo results for βL, βN and βP. For each estimate, the first row shows the median of

the vector of estimates, the second the mean, and the third the median of the standard deviation. The choice to display both median and mean stems from observing that for smaller simulated networks, the difference between the two is observed to be substantial even for βLin the densest networks (Appendix A, Table A1). For networks of

size N = 40 and N = 30, this difference is observed only for the sparser networks. For example, for N = 40 for the sparse network A3 the reported median of βL is 1.2051 and the mean is 10.1365. While the median is found to be

Table 3: Monte Carlo Results for N=40

A1 A2 A3 B1 B2 B3 med(βL) 1.0070 0.9894 1.2051 0.9830 1.0297 1.1582 mean(βL) 1.0030 1.0440 10.1365 1.0055 1.0327 8.3402 (8.2473) (1.3010) (0.0023) (8.4646) (1.8419) (0.0054) med(βN) 1.7955 37.3042 ** 1.7356 2.6409 ** mean(βN) 1.8791 20.8419 ** 1.8078 12.8318 ** (0.4230) (0.0000) ** (0.5792) (0.0025) ** med(βP) 0.9558 6.4218 ** 0.9408 1.2965 ** mean(βP) 0.9797 4.1868 ** 0.9396 3.0107 ** (4.8909) (0.000) ** (6.6947) (0.0326) **

very close to the true value β0= 1, the high value of the average suggests that there are very heavy outliers on the

right. This can be seen in the histogram of the estimate (Appendix B, Figure 27). The same instance is observed for sparse network B3, which displays median 1.1582 versus mean 8.3402. For N = 30, this difference becomes even bigger, exhibiting a mean of 14.4243 versus a median of 1.1482 for A3, and for N = 20, these estimates are not even available due to the near-singularity of the adjacency matrix. Specifically, for N = 40 this difference becomes almost negligible even for βN networks A1 and B1, while in the smaller networks this difference is observed to be

(19)

difference between the median and the mean is observed to become smaller as the size of the network is increased for both dense and sparse networks, it can be inferred that for larger networks this difference might possibly get close to 0 even for the sparser cases.

Furthermore, βN and βP are found to display bimodality in their histograms for the case of smaller networks,

for example in A1 and B1 (Appendix B, Figures 2,3, Figures 12,13). The latter findings could be an indication that for N = 40 however, this might not be the case. Taking a look at the histograms of the estimates for networks N = 40 for cases A1 and B1 (Appendix B, Figures 19-21, Figures 29-31), confirms that bimodality is indeed no longer observed for these N = 40 networks: the distribution of all the three estimates is close to a bell shaped distribution around one mean (with the exception of light right skewness for βN). Based on this and the previous

conclusion, it could be expected that for larger networks, the shape of the distribution of the other network types might also get close to a bell shaped one and lose bimodality.

Moreover, Table [3] allows for comparison between the three parameters of interest. It is observed that the estimates that differ the most from the true value are the ones of βN. While for the densest networks A1 and B1

the reported median values of βN − 1.7955 and 1.7356 respectively − are relatively not far from the true value 1,

for all other networks the estimates explode, up to being not obtainable numerically in the sparsest networks A3 and B3. Estimation of βP and βL on the other hand exhibit median and mean both close to the true value in all

other less sparse network configurations, for the exception βP in A2. Taking a look at Figure 1 shows that in the

A2 case, the βP indeed deviates from the true value. To visualise this, figure [1] displays the line graphs of the three

parameter estimates for networks A1 and A2, to compare how these behave in dense versus sparser networks. As expected, in the dense network A1 it can be observed that the βN presents a heavily irregular behaviour often with

large deviations from the true mean, whereas the βL and βP present a much more stable behaviour, remarkably

close around the true value throughout the 100 estimates. Nonetheless, it can be said that in the case of the dense network βP fluctuates around a mean that is approximately a unit away from the true one, while on the other

hand, for the sparser network this estimate explodes ranging between values 0 and 37 throughout the sample. This observation contributes to a possible conclusion that βN is not adequate for estimation of β0, while on the other

(20)

Figure 1: Dense (A1) versus sparse (A2) network

The line graph comparison for size N = 30 of network A1 versus A2 (Appendix B, Figure 4,8) also supports the claim that for a smaller network, the three estimators are less stable around the mean in the case of sparser networks. Specifically, it can be seen that for network A2, βN fluctuates extremely around its mean, peaking up to

value 38 in the A2 networks for both N = 30 and N = 40, versus peaking to 3 only in the dense A1 network of size N = 40 (Figure 1). As for network type B1, the distributions of βL and βP for N = 30 (Appendix B, Figures 11-13)

are also less close to a bell shaped one than for N = 40 networks (Appendix B, Figures 29-31). The improvements resulting from increasing network size, are due to the fact that the result about the distribution of ˆβ in (19) is an asymptotic one, meaning that is it is only approximately true for small and medium sized simulation frameworks. Since the TLE (Graham, 2017) was submitted to Monte Carlo studies of much larger size with a 1000 replications and networks sized powers of N = 100 and N = 200, it can be expected that for the networks considered here the results might be less attractive.

4.3 Test Results

The tests (16), (17) and (18) are performed by estimating both the type 1 error, namely the probability that H0is

rejected if H0 is true, and for comparison by also estimating the total number of rejections by p-value divided by

the simulation sample size. These results are shown in Table [4], with the p-value comparison in brackets. Ideally, both of these values are close to 0 for well specified models (Bain&Engelhardt, 1987).

(21)

It is observed that larger network size yields to better test results. For example, the hypothesis (16) for network A1 which is rejected 0% of the times for network N = 40, for the analogous network of size N = 30 it is rejected at a rate of 5% and at a rate of 44% for networks of size N = 20 (Appendix A, Table A2 and Table A4). Furthermore, it is observed that the rejection rates of βLfor networks A2 and B2 are still relatively low − 26% and 18% respectively

− compared to the analogous results for N = 30 which are 48% and 36%, and for N = 20 which are as high as 89% and 87%.

Table 4: Test Results for N=40

H0 A1 A2 A3 B1 B2 B3 (16) βL=1 0.00 0.26 0.99 0.00* 0.18 0.99 (0.01) 0.30 (0.61) (0.00) (0.19) (0.62) (17) βN=1 1.00 1.00 ** 0.96 1.00 ** (1.00) (1.00) ** (0.90) (1.00) ** (18) βP=1 0.07 0.98 ** 0.02 0.91 ** (0.09) 0.85 ** (0.02) (0.72) **

*do not reject normality, **near singular adjacency matrix

For the largest networks N = 40, testing of hypotheses (16) and (18) turn out to give the most promising results for the estimates of βLand βP in the densest networks A1 and B1, with type 1 error values of 0% in both networks

for βL, and 7% and 2% respectively for βP. On the other hand, testing of hypothesis (17) heavily supports rejection

of the null, meaning that there is enough evidence in the simulated sample that βN is not equal to the true value

β0= 1. This was to be expected since βN is found to exhibit very erratic behavior. This could be an indication

that these results could become even more attractive for larger networks, possibly also so for the sparser cases A3 and B3.

The analysis of network statistics, estimation outputs and test results, across the majority of the aforementioned findings a seem to display a pattern around one network characteristic, namely the density of the network. Firstly, it is observed that the density of a network simulated according to (1) does not depend on the size of the network when U is logistically distributed. Secondly, it was found that the difference between the true value and the median

(22)

of estimates is larger for sparser networks, across all six network designs and all three estimators. Thirdly, it was found that for all three hypotheses (16), (17) and (18) the rejection rate increases as the network gets sparser, yielding even opposing results for type 1 versus 3 networks under the same hypothesis testing procedure (0% versus 99% for testing of (16) for networks A1 versus A3 respectively). This could be evidence that for small networks the TLE procedure performs the best when applied to denser networks, and should be used with caution when applied to networks with density lower than ∼ 0.3. By the same argument, it could be said that the TLE performs best when applied to networks with transitivity (clustering coefficient) between 0.4 and 0.5.

5 Conclusion

The estimation of the homophily parameter in a network of undirected dyadic link formation by tetrad logit estimation as proposed by Graham (2017) can be concluded to be useful and relevant within the field of econometric network formation models. The theoretical basis of the tetrad logit estimation program can be summarised by two main theoretical pillars: i) the micro-economic theory of link surplus, by which an agent is assumed to link with another one only if the total surplus of doing so is positive, and ii) on the econometric theory of error estimation. The first is introduced as a basis for the simulation framework and is used to derive the structure of the TLE. The second takes care of the statistical estimation program and the computational aspect of the estimation of the TLE. As a consequence, it can be said that the tetrad logit estimation procedure on non-simulated observed networks which agents are assumed to link according to positive surplus, relies mostly on the econometric theory on the distribution of an observed network. The econometric nature of the model allows for it to distinguish itself from game theoretic models which rely less on the statistical structure of a given network, and more on the agent individual assumed behaviour.

For this reason, it is important to take a deeper look into the statistics of the TLE. Graham (2017) shows that the estimator proves itself to be legitimate, but only if the distribution of the network is intrinsically logistically distributed. The analysis of this paper is used to challenge and open a discussion around this assumption, since in observed networks it might turn out to restrict the usage of the TLE. The study proposed by this paper aims at understanding whether the TLE can be applied to networks which underlying distribution is found to be close to that of a normal distribution. The conclusions that can be drawn from the Monte Carlo analysis is that the tetrad

(23)

logit estimation program requires to be adjusted when the distribution is not found to follow a logistic one, and that specifically, in case it is found to be normally distributed, for dense networks an appropriate adaptation of the link function is to be sufficient. Thus, when the shock of the network is normally distributed, an alternative option to the TLE could be introduced by the name of tetrad probit estimator. To finalise the foregoing conclusions however further analysis and diagnostics need to be conducted, and other adaptation options of the estimation program to be explored.

In effect, it can be said that within the field of network science, much research has been done around networks with high density, while less attention has been given to sparse networks. As it was found in the results of this paper, many methods that are found to be accurate for dense networks, result to be inefficient or ineffective in networks where links are less frequent (Zhang&Martonosi, 2008). Graham (2017) also found that the other estimator − namely the JMLE − ceased to exist for the biggest portion of the non-dense networks. Techniques beyond regular maximum likelihood have been developed for sparser networks, such as penalized likelihood methods (Shojaie&Michailidis, 2010) or pseudo-likelihood methods (Höfling&Tibshirani, 2009), and could be explored as adaptation options to improve tetrad estimation in the case of sparser networks.

The limitations on the analysis of this paper are mainly caused by the small scale of the sample and small size of the simulated networks. The science of network formation has need to be applicable to very large networks, as observed networks across numerous fields are larger than thirty agents13_{. However, for the subgroup of dense}

networks − especially for those that are symmetric with uncorrelated degree heterogeneity − significant statistical conclusions were found. The results show that increasing the size of the network even by ten nodes at a time, resulted in great improvement in both estimation and testing of the hypotheses. This is a good indicator that the tetrad estimation procedure is to some degree resistant to sample size in certain network configurations.

To conclude, this paper is able to provide to the answer that was raised in the introduction, by providing statistical evidence that the TLE as derived by Graham (2017) is not directly applicable to networks in which the shock is not assumed to be logistically distributed, but that it has promising potential to be easily adaptable. This was found to be especially true for symmetric dense networks with uncorrelated degree heterogeneity, in which adaptation by changing of the link function of its estimation program is found to suffice in the case that the

13_{The choice to focus on small networks is a consequence of the fact that the computation of the variance estimator is heavy, and}

(24)

distribution of the error term is gaussian. Areas for further work include simulation and testing with distributions other than logistic and normal, and comparison with link functions other than the logit and the probit. Most importantly, the study proposed in this paper demands for replication with bigger sized networks.

(25)

Bibliography

Bain, L. J., & Engelhardt, M. (1987). Introduction to probability and mathematical statistics. Brooks/Cole.

Bazán, J. L., Romeo, J. S., & Rodrigues, J. (2014). Bayesian skew-probit regression for binary response data. Brazilian Journal of Probability and Statistics, 28(4), 467-482.

Bondy, J. A., & Murty, U. S. R. (1976). Graph theory with applications (Vol. 290). London: Macmillan.

Borshchev, A ., & Filippov, A. (2004). From System Dynamics and Discrete Event to Practical Agent Based Modeling: Reasons, Techniques, Tools. The 22nd International Conference of the System Dynamics Society

Boucher, V. (2015). Structural Homophily, International Economic Review, 56 (1), 235â264. [827]

Box G., Cox D., (1964). An analysis of transformations.

Chambers, E. A., & D. R. Cox (1967). Discrimination between alternative binary response models, Biometrika, 54, 573â578.

Desmarais, B. A., & Cranmer, S. J. (2012). Statistical inference for valued-edge networks: The generalized expo-nential random graph model. 7(1), 30136.

Erdos, P., & Rényi, A. (1960). On the evolution of random graphs. Publ. Math. Inst. Hung. Acad. Sci, 5(1), 17-60.

Ferguson, & Thomas, S. (2005). UâStatistics, Notes for Statistics 200c. Spring.

Figueredo, G., & Aickelin, U. (2011). Comparing System Dynamics and Agent-Based Simulation for Tumour Growth and its Interactions with Effector Cells. 5-22.

Graham, B. (2017). An econometric model of network formation with degree heterogeneity. Econometrica : Journal of the Econometric Society, an Internat. Society for the Advancement of Economic Theory in Its Relation to Statistics and Mathematics, 85(4), 1033-1063.

Gualdani, C. (2017). Econometric analysis of network formation models. Doctoral Thesis, Doctoral thesis, UCL (University College London).

Gilchrist R., Green M. (1984). Fitting parametric link functions in generalised linear models, Computational Statis-tics Data Analysis, Volume 2, Issue 1. 37-49,

(26)

Guerrero, V. M., & Johnson, R. A. (1982). Use of the Box-Cox transformation with binary response models. Biometrika, 69(2), 309-314.

Heij C., Boer P., Franses, P. H., Kloek T., & Dijk, H. K. (2004). Econometric Methods with Applications in Business and Economics

HĂśfling, H., & Tibshiran, R. (2009). Estimation of Sparse Binary Pairwise Markov Networks using Pseudo-likelihoods pp. 883–906, 2009.

Koenker, R., & Yoon, J. (2009). Parametric links for binary choice models: A FisherianâBayesian colloquy. Journal of Econometrics, 152(2), 120-130.

Leung, M. P. (2015). Two-step estimation of network-formation models with incomplete information. Journal of Econometrics, 188(1), 182-195.

Lu, Y., Jerath, K., & Singh, P. V. (2013). The emergence of opinion leaders in a networked online community: A dyadic model with time dynamics and a heuristic for fast estimation. Management Science, 59(8), 1783-1799.

McPherson, M., Smith-Lovin, L., & Cook, J. M. (2001). Birds of a feather: Homophily in social networks. Annual review of sociology, 27(1), 415-444.

Mele, A. (2017). A Structural Model of Dense Network Formation. Econometrica, 85(3), 825-850.

Newman, M., Watts, D., & Strogatz, S. (2002). Random Graph Models of Social Networks. Proceedings of the National Academy of Sciences of the United States of America, 99(3), 2566-2572.

Pregibon, D. (1980). Goodness of link tests for generalized linear models. Applied statistics, 15-14.

Ridder, G., & Sheng, S. (2015). Estimation of large network formation games. Working papers, UCLA.

Sheng, S. (2012). Identification and estimation of network formation games.

Shojaie, A., Michailidis, G. (2010). Penalized likelihood methods for estimation of sparse high-dimensional directed acyclic graphs. Biometrika, Volume 97, Issue 3, pp.519â538

(27)

Zhang, P. & Martonosi, M. (2008). LOCALE: Collaborative Localization Estimation for Sparse Mobile Sensor Networks. International Conference on Information Processing in Sensor Networks. pp. 195-206.

(28)

APPENDIX

A: Smaller networks

Table A1: Monte Carlo Results for N=20

A1 A2 A3 B1 B2 B3 med(βL) 1.0139 1.1455 ** 1.0118 1.0507 ** mean(βL) 2.1784 6.0214 ** 1.7846 4.4628 ** (0.8124) (0.1110) (**) (1.0418) (0.1739) (**) med(βN) 2.2991 37.3042 ** 1.8674 37.3042 ** mean(βN) 15.4609 33.4459 ** 11.1274 31.1617 ** (0.0165) (NA) (**) (0.0475) (NA) (**) med(βP) 1.0499 6.4218 ** 1.0070 6.4218 ** mean(βP) 2.4382 5.7509 ** 2.1886 4.8642 ** (0.3678) (NA) (**) (0.5651) (NA) (**)

**near singular adjacency matrix

Table A2: Test Results for N=20

H0 A1 A2 A3 B1 B2 B3 (16) βL=1 0.44 0.89 ** 0.44 0.87 ** (0.38) (0.56) (**) (0.39) (0.51) (**) (17) βN=1 0.97 0.97 ** 0.97 0.94 ** (0.97) (0.98) (**) (0.97) (0.91) (**) (18) βP=1 0.53 0.97 ** 0.42 0.94 ** (0.48) (0.89) (**) (0.39) (0.76) (**)

(29)

Table A3: Monte Carlo Results N=30 A1 A2 A3 B1 B2 B3 med(βL) 1.0135 0.9554 1.1482 0.9892 0.9764 1.2988 mean(βL) 1.0200 1.0583 14.4243 1.0044 1.0284 15.5698 (3.1363) (0.6682) (NA) (3.6615) (0.7563) (0.0021) med(βN) 1.9414 37.3042 ** 1.7744 37.3042 ** mean(βN) 3.7185 26.3628 ** 2.1927 26.0222 ** (0.1384) (0.0000) ** (0.2015) (0.0000) (**) med(βP) 0.9818 6.4218 ** 0.9417 6.4218 ** mean(βP) 1.2035 5.0735 ** 1.0130 4.1977 ** (2.0033) (0.0000) ** (2.6851) (NA) (**)

**near singular adjacency matrix

Table A4: Test Results for N=30

H0 A1 A2 A3 B1 B2 B3 (16) βL=1 0.05* 0.48 0.94 0.05* 0.36 0.95 (0.08) (0.34) (0.59) (0.07) (0.28) (0.59) (17) βN=1 0.99 1.00 ** 0.98 0.99 ** (0.99) (1.00) (**) (0.98) (0.98) (**) (18) βP=1 0.28 1.00 ** 0.11* 0.95 ** (0.28) (1.00) (**) (0.11) (0.89) (**)

(30)

B: Figures

Network A1: N30

Figure 1: Histogram BL:A1:30

(31)

Figure 3: Histogram BP:A1:30

(32)

Network A2: N30

(33)

(34)

Network A3: N30

(35)

Network B1: N30

Figure 11: Histogram BL:B1:30

(36)

Figure 13: Histogram BP:B1:30

(37)

Network B2: N30

(38)

Figure 17: Histogram BN:B2:30

(39)

Network A1: N40

(40)

(41)

Network A2: N40

(42)

(43)

Network A3: N40

(44)

Network B1: N40

(45)

(46)

Network B2: N40

(47)

An Econometric Model of Network Formation