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.
72,837 characters · 21 sections · 50 citation commands
Reconciling econometrics with continuous maximum-entropy network models
\email{[email removed]}
\pacs{89.75.Fb; 02.50.Tt; 89.65.Gh}
Over the last couple of decades, the growth of network science has impacted several disciplines by establishing new, empirical facts about many, real-world systems, in terms of the structural, network properties that are typically found in those systems. In the context of trade economics, the growing availability of data about import/export relationships among world countries has prompted researchers to explore and model the architecture of the international trade network, or World Trade Web (WTW) Barigozzi2010,Fronczak2012a,Fronczak2012b,Serrano2003,Garlaschelli2004,Garlaschelli2005,Fagiolo2010,Fagiolo2008a,Fagiolo2008b,Schweitzer2009,Vitali2011,Schiavo2010.
This approach has complemented, and in many ways enriched, the traditional econometric exercise of modelling individual trade flows, i.e. relating the volume of individual trade exchanges to the most relevant covariates (generally, macroeconomic factors) they may depend on. The earliest example of an econometric model for international trade is the celebrated Gravity Model (GM) Tinbergen1962 that predicts that the expected value $\langle w_{ij}\rangle_\text{GM}$ of the trade volume $w_{ij}$ from country $i$ to country $j$ can be expressed via the econometric function
where $\omega_i\equiv{\text{GDP}_i}/{\overline{\text{GDP}}}$ is the GDP of country $i$ divided by the arithmetic mean of the GDPs of all countries, $d_{ij}$ is the geographic distance between (generally the capitals of) countries $i$ and $j$ and $\underline{\psi}\equiv(\tau,\beta_1,\beta_2,\gamma)$ is a vector of parameters (notice that $\beta_1\equiv\beta_2$ when the direction of the exchanges is disregarded, as we will, throughout the paper, and $\tau$ takes care of dimensional units). Equation (ref) has a long tradition in successfully explaining the existing (i.e. positive) trade volumes between pairs of countries. Outside economics, the gravity equation has also been extensively employed in studies concerning transportation, migration Wilson1969 and the maximization of utility functions constrained to satisfy requirements on the rate of information acquisition Wang2021.
Although accurate in reproducing the positive trade volumes, the traditional GM cannot replicate structural network properties, unless the topology is completely fixed via a separate approach Fagiolo2010b. Indeed, if $\mathbf{W}$ denotes the weighted adjacency matrix of the WTW, where the entry $w_{ij}$ represents the trade volume from country $i$ to country $j$, Eq. (ref) predicts that the expected matrix $\langle\textbf{W}\rangle_\text{GM}$ has no off-diagonal zeroes, i.e. an expected positive trade relationship exists between all countries. This means that, when interpreted as an expected value of a regression with small, symmetric, zero-mean (e.g. Gaussian) noise, Eq. (ref) predicts a fully connected network: if $a_{ij}$ denotes the generic entry of the binary adjacency matrix $\mathbf{A}=\Theta[\mathbf{W}]$ (equal to $a_{ij}=1$ if a positive trade volume from country $i$ to country $j$ is present, i.e. $w_{ij}>0$, and equal to $a_{ij}=0$ if a zero trade is, instead, observed, i.e. $w_{ij}=0$), then Eq. (ref) predicts almost surely $a_{ij}=1$, $\forall\:i\neq j$. This result is in obvious contrast with empirical data, which show that the WTW has a rich topological architecture, characterized by a broad degree distribution, (dis)assortative and clustering patterns and other properties Barigozzi2010,Fronczak2012a,Fronczak2012b,Serrano2003,Garlaschelli2004,Garlaschelli2005,Fagiolo2010,Fagiolo2008a,Fagiolo2008b,Schweitzer2009,Vitali2011,Schiavo2010.
In order to overcome such a limitation, the plain Gravity Model needs to be `dressed' with a probability distribution $Q(\mathbf{W})$ that produces $\langle\textbf{W}\rangle_\text{GM}$ as the expected value while, at the same time, accounts for null outcomes as well (i.e. those entries reading $w_{ij}=0$ and representing missing links in the network) Helpman2008. Clearly, the support $\mathbb{W}$ of the probability $Q(\mathbf{W})$ should not include matrices with negative numbers. From the Sixties on, the GM has indeed been interpreted as the expected value of a probability distribution whose functional form needs to be determined.
Trade econometrics models the tendency of countries to establish trade relationships relating it to accepted, macroeconomic determinants (the so-called `factors') such as GDP and geographic distance, as in the expression of Eq. (ref). Econometricians have considered increasingly flexible distributions, the most recent versions of them being capable of disentangling the estimation of the presence of a trade exchange from the estimation of the traded amount. This has led to the definition of two distinct classes of models, i.e. zero-inflated models Burger2009 and hurdle models Long1997. Zero-inflated (ZI) models have been introduced to model the following two scenarios: the possibility of exchanging a zero amount of goods even after having established a trade partnership (e.g. because of a limited trading capacity); the possibility of establishing a very small amount of trade (in fact, so small to be compatible with statistical noise and, as such, removed). A general drawback of employing ZI models is that of predicting sparser-than-observed network structures Marzio2022; moreover, only discrete distributions (specifically, either the Poisson or the negative-binomial one Burger2009) have been considered, so far, to carry out the proper weight-estimation step. Hurdle models, introduced to overcome the limitations affecting zero-inflated models, can predict zeros only at the first step Long1997: in any case, the presence of links is established by employing either a logit or a probit estimation step.
Network science has tackled the aforementioned inference problem using techniques rooted in statistical physics. The most prominent examples descend from the Maximum-Entropy Principle (MEP) Jaynes1957a,Jaynes1957b,Jaynes1982 applied to network ensembles DSBook,Cimini2019, prescribing to maximize Shannon entropy Shannon,Cover in a constrained fashion to obtain the maximally unbiased distribution of networks compatible with a chosen set of structural constraints. This approach is formally equivalent to the construction of so-called Exponential Random Graphs (ERGs) for social network analysis Wasserman1994 but differs in the typical choice of the constraints: in particular, when the enforced constraints are local, such as the degree (number of links) and the strength (total link weight) of each node, maximum-entropy network models have been shown to successfully replicate both the topology and the weights of many economic and financial networks, including the WTW Garlaschelli2004,Garlaschelli2005,Squartini2011a,Squartini2011b,Mastrandrea2014b,Squartini2014,Almog2015,Almog2017,Almog2019,Marzio2022. The entire framework can also accommodate possibly degenerate, discrete-valued, single- or multi-edges Sagarra2015.
Although maximum-entropy models have been also studied from an economic perspective (see Bargigli2014 for a discussion of the economic relevance of the constraints defining the Poisson and the geometric network models), it is only recently that progress has been made to reconcile the above two approaches, allowing for economic factors parametrizing the maximum-entropy probability distribution producing links and weights Garlaschelli2004,Garlaschelli2005,Squartini2014,Almog2015,Almog2017,Almog2019,Marzio2022 or by introducing network-related statistics into otherwise purely econometric models Herman2022. On one hand, the novel framework enriches the methods developed by network scientists with an econometric interpretation; on the other, it enlarges the list of candidate distributions usable for econometric purposes.
With this contribution, we refine the theoretical picture provided in a companion paper Marzio2022, introducing models to infer the topology and the weights of undirected networks defined by continuous-valued data. In order to do so, we present a theoretical, physics-inspired framework capable of accommodating both integrated and conditional, continuous models, our goal being threefold: 1) testing the performance of both classes of models on the WTW in order to understand which one is best suited for the task; 2) offering a principled derivation of currently available, conditional, econometric models; 3) enlarging the list of continuous-valued distributions to be used for econometric purposes. From an econometric point of view, our work moves along the methodological guidelines defining the class of Generalized Linear Models (GLMs) Nelder1972 while enriching it with distributions defined by both econometric and structural parameters. From a statistical physics point of view, our work expands the class of maximum-entropy network models DSBook or weighted ERGs Wasserman1994 and endows them with macroeconomic factors replacing certain model parameters.\\
The rest of the paper is organized as follows: in Sec. (ref), after introducing the basic quantities, we derive the class of conditional models; in Sec. (ref) we derive the class of integrated models; in Sec. (ref) we apply all models to the analysis of WTW data; in Sec. (ref) we discuss the results and provide our concluding remarks.
Discrete maximum-entropy models can be derived by performing a constrained maximization of Shannon entropy Jaynes1957a,Jaynes1957b,Jaynes1982. However, unlike the companion paper Marzio2022, our focus, here, is on continuous probability distributions. In such a case, mathematical problems are known to affect the definition of Shannon entropy and the resulting inference procedure. To restore the framework, one has to introduce the Kullback-Leibler (KL) divergence $D_\text{KL}(Q||R)$ of a distribution $Q$ from a prior distribution $R$ and re-interpret the maximization of the entropy of $Q$ as the minimization of $D_\text{KL}(Q||R)$ from a given prior distribution $R$. In formulas, the KL divergence is defined as
where $\mathbf{W}$ is one of the possible values of a continuous random variable (in our setting, an entire network with continuous-valued link weights), $\mathbb{W}$ is the set of possible values that $\mathbf{W}$ can take, $Q(\mathbf{W})$ is the (multivariate) probability density function to be estimated and $R(\mathbf{W})$ plays the role of prior distribution, the divergence of $Q(\mathbf{W})$ from which must be minimized. Such an optimization scheme embodies the so-called Minimum Discrimination Information Principle (MDIP), originally proposed by Kullback and Leibler Kullback1951 and implementing the idea that, given a prior distribution $R(\mathbf{W})$ and new information that becomes available, an updated distribution $Q(\mathbf{W})$ should be chosen in order to make its discrimination from $R(\mathbf{W})$ as hard as possible. In other words, the MDIP demands that new data produce an information gain that is as small as possible. The use of the KL divergence is widespread in the fields of information theory Shannon and machine learning Goodfellow2014, e.g. as a loss function within the Generative Adversarial Network (GAN) scheme (the aim of the `generating' neural network being that of producing samples that cannot be distinguished from those constituting the training set by the `discriminating' neural network).
In order to introduce the class of conditional models, we write the posterior distribution $Q(\mathbf{W})$ as
where $\mathbf{A}$ denotes the adjacency matrix for the binary projection of the weighted network $\mathbf{W}$. The above equation allows us to split the KL divergence into the following sum of three terms
where
is the Shannon entropy of the probability distribution describing the binary projection of the network structure,
is the conditional Shannon entropy of the probability distribution of the weighted network structure given the binary projection and
is the cross entropy quantifying the amount of information required to identify a weighted network sampled from the distribution $Q(\mathbf{W})$ by employing the distribution $R(\mathbf{W})$. When continuous models are considered, $S(Q_\bot|P)$ is defined by a first sum running over all the binary configurations within the ensemble $\mathbb{A}$ and an integral over all the weighted configurations that are compatible with each, specific, binary structure - embodied by the adjacency matrix $\mathbf{A}$, i.e. such that $\mathbb{W}_\mathbf{A}=\{\mathbf{W}:\Theta[\mathbf{W}]=\mathbf{A}\}$).
The expression for $S(Q,R)$ can be further manipulated as follows. Upon separating the prior distribution itself into a purely binary part and a conditional, weighted one, one can write
an expression that allows us to write $S(Q,R)$ as
which, in turn, allows the KL divergence to be rewritten as
i.e. as a sum of two terms, one of which involves conditional distributions; specifically,
with $T(\mathbf{A})$ representing the binary prior and $R(\mathbf{W}|\mathbf{A})$ representing the conditional, weighted one. In what follows, we will deal with completely uninformative priors: this amounts at considering the somehow `simplified' expression
with
The (independent) constrained optimization of $S(P)$ and $S(Q_\bot|P)$ represents the starting point for deriving the members of the class of conditional models.
The functional form controlling for the binary part of conditional models can be derived by carrying out a constrained maximization of the binary Shannon entropy
leading to a probability mass function reading
where the functional form of the Hamiltonian reads $H(\mathbf{A})=\sum_{i<j}\alpha_{ij} a_{ij}$. This choice induces a factorization of the probability mass function $P(\mathbf{A})$, which becomes
i.e. a product of a number of Bernoulli-like probability mass functions with $p_{ij}=\frac{x_{ij}}{1+x_{ij}}$, $\forall\:i<j$, where $x_{ij}=e^{-\alpha_{ij}}$ is the Lagrange multiplier controlling for the generic entry of the adjacency matrix $\mathbf{A}$.
In what follows, we will consider the specification of the tensor-like Hamiltonian introduced above reading $\alpha_{ij}=\alpha_i+\alpha_j$, $\forall\:i<j$, a choice inducing the Undirected Binary Configuration Model (UBCM), characterized by the following, pair-specific probability coefficient
and ensuring the entire degree sequence of the network at hand to be reproduced.
The econometric reparametrization of the UBCM can be achieved by posing $x_i\equiv\sqrt{\delta}\omega_i$, $\forall\:i$, a choice inducing the so-called fitness model (FM), characterized by the pair-specific probability coefficient
and requiring the estimation of a global parameter only, i.e. $\delta$. The FM represents a particular case of the logit model Walker1967, being defined by a vector of external properties (the `fitnesses') that replace the information provided by some kind of (otherwise) purely structural properties Caldarelli2003,Garlaschelli2004: in fact, Eq. (ref) can be equivalently rewritten as $\text{logit}\left[p_{ij}^\text{FM}\right]\equiv e^{\underline{X}\cdot\underline{\phi}}$ with $\underline{X}\equiv[1,\ln(\omega_i\omega_j)]$ and $\underline{\phi}\equiv[\ln\delta,1]$. The global constant $\delta$ can be determined by imposing the total number of links as the only constraint. Remarkably, the fitness model has been proven to reproduce the (binary) properties of a wide spectrum of real-world systems Squartini2011a, Garlaschelli2005 as accurately as the UBCM, although requiring much less information.
In what follows, we will consider both the UBCM and the FM specifications.
The constrained maximization of $S(Q_\bot|P)$ proceeds by specifying the following set of weighted constraints
the first condition ensuring the normalization of the probability distribution and the vector $\{C_\alpha(\mathbf{W})\}$ representing the `proper' set of weighted constraints (weights are, now, treated as continuous random variables, i.e. $w_{ij}\in\mathbb{R}^+_0$, $\forall\:i<j$). They induce the distribution reading
where $H(\mathbf{W})=\sum_\alpha\psi_\alpha C_\alpha$ is the so-called Hamiltonian, listing the constrained quantities, and $Z_\mathbf{A}=\int_{\mathbb{W}_\mathbf{A}}e^{-H(\mathbf{W})}d\mathbf{W}$ is the partition function, conditional on the `fixed topology' $\mathbf{A}$.
The explicit functional form of $Q(\mathbf{W}|\mathbf{A})$ can be obtained only once the functional form of the constraints has been specified. In what follows, we will deal with the Hamiltonian reading
with the Lagrange multipliers $\left(\beta_{ij},\xi_{ij},\gamma_{ij}\right)$ satisfying the following requirements:
Let us start by considering the simplest, conditional model, defined by the positions $\gamma_{ij}=\xi_{ij}=0$ and inducing the Hamiltonian
inserting the expression above into Eq. (ref) leads to the distribution
and each node pair-specific distribution induces a (conditional) expected weight reading
From a purely topological point of view, constraining each weight and their total sum is redundant. However, this is no longer true when turning the conditional, exponential model into a proper econometric one. Its econometric reparametrization should be consistent with the literature on trade, stating that the weights are monotonically increasing functions of the gravity specification, i.e. $\langle w_{ij}\rangle_\text{GM}=e^{\rho+\beta\cdot\ln(\omega_i\omega_j)+\gamma\cdot\ln(d_{ij})}\equiv z_{ij}$, $\forall\:i<j$ and with $e^\rho\equiv\tau$; for this reason, the link function usually associated with the exponential distribution prescribes to identify the linear predictor with the inverse of the purely econometric parameter of the model, i.e.
a position that turns Eq. (ref) into
notice that the only structural constraint is, now, represented by the total weight (see also Appendix A).
Let us, now, consider a different Hamiltonian, constraining each weight, their total sum and the sum of their logarithms, i.e.
it induces the distribution reading
each node pair-specific distribution is characterized by a (conditional) expected weight reading
and by a (conditional) expected logarithmic weight reading
where the function $\psi(x)=\Gamma'(x)/\Gamma(x)$ is the so-called digamma function.\\
Such a model can be turned into a proper, econometric one by considering the inference scheme of the gamma model with inverse response, which allows us to identify the linear predictor with the inverse of the purely econometric parameter of the model, i.e.
a position that, in turn, leads to the expressions
(allowing the conditional, exponential model to be recovered in case $\xi_{0}=0$, i.e. when the constraint on the sum of the logarithms of the weights is switched-off) and
(see also Appendix A).
Constraining a slightly more complex function of the weights, i.e. their logarithm, leads to the Hamiltonian
which, in turn, induces the distribution
where $m_{ij}$ is the minimum, node pair-specific weight allowed by the model. Each node pair-specific distribution is characterized by a (conditional) expected weight reading
Such a model can be turned into a proper, econometric one by considering the positions
ensuring that the expected weights are monotonically increasing functions of the gravity specification and leading to the expression
where $w_{min}$ is the empirical, minimum weight (see also Appendix A).\\
Let us explicitly notice that the derivation of the gamma and Pareto distributions within the maximum-entropy framework has been already studied in Visser2013; here, however, we aim at making a step further, by individuating a suitable redefinition of these models parameters capable of turning them into proper, econometric ones.
Adding a global constraint on (a function of) the total variance of the logarithms of the weights to the Hamiltonian defining the Pareto model leads to the expression
the Hamiltonian above induces a distribution reading
each node pair-specific distribution is characterized by a (conditional) expected weight reading
by a (conditional) expected, logarithmic weight reading
and by a (conditional) logarithmic weight whose squared expectation reads
Such a model can be turned into a proper, econometric one by considering the position
ensuring that the expected weights are monotonically increasing functions of the gravity specification and leading to the expressions
(see also Appendix A).
MDIP can be also implemented in a straightforward way, by carrying out a constrained optimization of $D_\text{KL}(Q||R)$. In this second case, the following set of constraints
can be specified, with obvious meaning of the symbols. Differentiating the corresponding Lagrangean functional with respect to $Q(\mathbf{W})$ and equating the result to zero leads to
where $H(\mathbf{W})=\sum_\alpha\psi_\alpha C_\alpha$ is, again, the Hamiltonian and $Z=\int_\mathbb{W}e^{-H(\mathbf{W})}d\mathbf{W}$ is the `integrated' partition function.
The explicit functional form of $Q(\mathbf{W})$ can be obtained only once the functional form of both the prior distribution and the constraints has been specified as well. In what follows, we will deal with completely uninformative priors, a choice that amounts at considering the simplified expression
notice that the result above could have been also derived by carrying out a constrained minimization of
i.e. of (minus) the functional named differential entropy into which the KL divergence `degenerates' in case completely uninformative priors are considered.
Let us, now, specify the functional form of the constraints. In what follows, we will deal with a specific instance of the generic Hamiltonian
in particular, one could pose $\alpha_{ij}\equiv \alpha_{0}$, a choice that would lead to constrain the total number of links, or $\alpha_{ij}\equiv\alpha_i+\alpha_j$, a choice that would lead to constrain the whole degree sequence. If not specified otherwise, in what follows we will employ the second functional form and pose $\beta_{ij}\equiv\beta_0+\beta_{ij}$, where $\beta_0$ is the Lagrange multiplier associated with the total weight and $\beta_{ij}$ encodes the dependence on purely econometric quantities. Our choices induce the Hamiltonian of the so-called integrated exponential model, i.e.
that leads to the distribution
where $x_i\equiv e^{-\alpha_i}$. The generic node pair-specific distribution induces a probability for nodes $i$ and $j$ to be connected reading
besides, the corresponding expected weight reads
Equation (ref) clarifies why the models considered in the present section are classified as `integrated': each node pair-specific probability of connection is a function of the parameters controlling for both topological and weighted properties. Models of the kind are, thus, capable of `integrating' information concerning a network structure with information concerning its weights, hence employing them in a joint fashion to define both inference steps.
The recipe for the econometric reparametrization of the integrated exponential model can read as the one of its conditional counterpart, i.e.
a position that turns Eq. (ref) into
where
(see also Appendix B).
The effectiveness of the two classes of models considered in the present paper to reproduce the topological properties of the World Trade Web has been tested on two different datasets, i.e. the Gleditsch one (covering 11 years, from 1990 to 2000 Gled2002) and the BACI one (covering 11 years, from 2007 to 2017 Baci2014). To carry out our analyses, we have built the ensemble induced by each model as follows. First, the presence of a link connecting any two nodes $i$ and $j$ is established with probability $p_{ij}$. Numerically, a real number $u$ is drawn from the uniform distribution $U[0,1]$ with unit support and compared with $p_{ij}$: if $u\leq p_{ij}$, then $i$ and $j$ are linked, otherwise they are not. Once the presence of a link is established, it is loaded with a weight by employing the inverse transform sampling technique: another random variable $\eta$, uniformly distributed between 0 and 1, is set equal to the value of the complementary cumulative distribution
and, inverting the equation $F(v_{ij})=\eta$, one obtains the value of the random variable $v_{ij}$ to be assigned as link weight to the pair $i<j$. Each ensemble is sampled repeatedly to obtain $10^4$ configurations. The error accompanying the estimate of any quantity of interest is quantified via the confidence intervals (CI) induced by the ensemble distribution of the quantity itself.
Let us consider two measures of goodness-of-fit, i.e. the reconstruction accuracy $\text{RA}^s_m$ and the Kolmogorov-Smirnov (KS) compatibility frequency $f^s_m$: while $\text{RA}^s_m$ is defined as the percentage of node-specific values of the statistic $s$ falling within the CIs, induced by model $m$, at the significance level of $5\%$, $f^s_m$ measures the percentage of times the distribution of a given, expected statistics $s$, under model $m$, is compatible with the empirical one, according to the two-sample Kolmogorov-Smirnov test; for example, a value $f^s_m=0.8$ would indicate that the distribution of the expected values of the statistics $s$, under model $m$, is found to be compatible with the empirical one on the $80\%$ of the years in the dataset, at the significance level of $5\%$.
The network statistics for which the values $\text{RA}^s_m$ and $f^s_m$ have been computed are the degree sequence
(which gives information about the tendency of node $i$ to connect to other trade partners), the average nearest neighbors degree
(which gives information about the degree correlations), the clustering coefficient
(which counts the percentage of node $i$'s partners that are also partners themselves). For what concerns the weighted statistics, we have considered the strength sequence
(which gives information about the trade flow of a country), the average nearest neighbors strength
(which gives information about the strength correlations), the weighted clustering coefficient
(that weighs the closed triangular patterns that node $i$ establishes with other trade partners).
Table (ref) lists the values of $f^s_m$ for both binary and weighted network statistics. For what concerns the binary statistics, we report the performance of three different models, i.e. the UBCM, the FM and the integrated exponential model (denoted as I-Exp).
For what concerns the Gleditsch dataset, compatibility is observed for every year; for what concerns the BACI dataset, instead, this is no longer true: in fact, the FM outputs predictions that are not compatible with the empirical values for a large number of years and irrespectively from the considered quantity; the UBCM and the I-Exp (i.e. the models constraining the degrees), instead, output predictions whose compatibility depends on the considered quantity: higher-order statistics are the ones for which the two aforementioned models `fail' to the larger extent. Overall, these results lead us to prefer the UBCM as the `first step-algorithm' of our conditional models.
Let us, now, comment on the performance of our models in reproducing weighted statistics. As it can be appreciated upon looking at Table (ref), the only models outputting predictions whose distributions are compatible with the empirical analogues are the integrated exponential one, the conditional exponential one and the conditional gamma one. On the other hand, employing only logarithmic constraints (as for the conditional Pareto model and the conditional log-normal model) does not help improving the accuracy of the description of the system at hand.
So far, we have inspected the compatibility of the distributions of the empirical values of each network statistics with the ones of their expected values under each of our models. Let us, now, quantify the extent to which each model is able to recover node-wise information by computing the $\text{RA}^s_m$ values. Figure (ref) shows the temporal average of the latter ones (i.e. across the years covered by our datasets), with the whiskers representing their variation, i.e. an indication of the stability of each model performance.
For what concerns the binary statistics (see Fig. (ref)a), both the UBCM and the I-Exp perform quite well in reproducing them; on the other hand, the performance of the FM is much poorer. For what concerns the weighted statistics (see Fig. (ref)b), only the average nearest neighbors strength is satisfactorily recovered by the (integrated and conditional) exponential models and the conditional gamma one. Still, they are found to perform poorly on the other statistics, i.e. the strength, that is only recovered in distribution on both datasets, and the weighted clustering coefficient, that is only recovered in distribution on the BACI dataset. For what concerns the lognormal model, it performs better than competitors in reproducing the strength and the weighted clustering coefficient on the Gledistch dataset but worse on the BACI dataset, causing its behavior to be dataset-dependent.
Finally, let us inspect the reconstruction accuracy of our models for what concerns our networks link weights. Specifically, let us define $\text{RA}^w_m$, i.e. the percentage of empirical weights falling within the CI, induced by model $m$, at the significance level of $5\%$ Parisi2020. Here, we have proceeded numerically, i.e. by considering the $2.5$ and the $97.5$ percentiles induced by the ensemble distribution of each node pair-specific weight. As Fig. (ref) shows, all models perform quite well in reproducing the weights (across all years, on both datasets) with the only exception of the conditional Pareto model. Overall, the best-performing model on the Gleditsch dataset is the integrated exponential one while the best-performing model on the BACI dataset is the conditional gamma model.
The UBCM and the integrated exponential model perform similarly in reproducing the binary statistics, on both datasets. Let us, now, compare them in reproducing the four indicators composing the so-called confusion matrix, i.e. the true positive rate $\langle\text{TPR}\rangle=\langle\text{TP}\rangle/L=\sum_{i<j}a_{ij}p_{ij}/L$ (measuring the percentage of links correctly recovered by a given reconstruction method), the specificity $\langle\text{SPC}\rangle=\langle\text{TN}\rangle/(N(N-1)/2-L)=\sum_{i<j}(1-a_{ij})(1-p_{ij})/(N(N-1)/2-L)$ (measuring the percentage of zeros correctly recovered by a given reconstruction method), the positive predictive value $\langle\text{PPV}\rangle=\langle\text{TP}\rangle/\langle L\rangle=\sum_{i<j}a_{ij}p_{ij}/\langle L\rangle$ (measuring the percentage of links correctly recovered by a given reconstruction method with respect to the total number of links predicted by it) and the accuracy $\langle\text{ACC}\rangle=(\langle\text{TP}\rangle+\langle\text{TN}\rangle)/N(N-1)/2$ (measuring the overall performance of a given reconstruction method in correctly placing both links and zeros).
The results are reported in Table (ref), that shows the increments of the four indicators, defined as $\Delta_X=\langle X\rangle_\text{I-Exp}-\langle X\rangle_\text{UBCM}$ with $X=\text{TPR},\:\text{SPC},\:\text{PPV},\:\text{ACC}$. Notice that each entry of the table is positive, a result signalling that the integrated exponential model steadily performs better than the UBCM. This is further confirmed by the (non-parametric) Wilcoxon rank-sum test on the ensemble distributions of the statistics to compare: all increments are significant, at the $1\%$ level.
Let us now rigorously test if constraining the entire degree sequence $\{k_i\}_{i=1}^N$ leads to a significantly better description of our data than that obtainable by just constraining the total number of links $L$.
Upon solving the model constraining the entire degree sequence and the one constraining the total number of links, we are able to construct a vector reading $(X_i^m, Y_i^m)$ where $X_i^m$ is either $\text{RA}^s_m$ or $f^s_m$ for the $i$-th statistics under the `$L$-constrained version' of model $m$; on the other hand, $Y_i^m$ is either $\text{RA}^s_m$ or $f^s_m$ for the $i$-th statistics under the `$k$-constrained version' of model $m$ - naturally, both values have been considered for the same year, keeping the same set of weighted constraints. Pairing statistics as described above allows us to employ the (non-parametric) Wilcoxon signed-rank test for testing the hypotheses $\text{RA}^s_k\leq\text{RA}^s_L$ and $f^s_k\leq f^s_L$, i.e. that the models just constraining $L$ perform better, in reproducing the statistics $s$, than those constraining the entire degree sequence.
Our results let us conclude that, for both datasets, constraining the degree sequence leads to a significant improvement, at the level of $5\%$, of the reconstruction accuracy of the average nearest neighbors degree, the clustering coefficient, the strengths and the average nearest neighbors strength; on the other hand, constraining the degree sequence does not lead to any significant improvement of the reconstruction accuracy of the weighted clustering coefficient. For what concerns the KS compatibility frequency, a significant improvement, at the level of $5\%$, is observed in the description accuracy of the average nearest neighbors degree, the clustering coefficient and the average nearest neighbors strength.
Let us, now, compare the performance of our models in a more general fashion. To this aim, let us consider the Akaike Information Criterion (AIC) Akaike1973, reading
where $k$ is the number of free parameters of the model and $\mathcal{L}_m$ is its log-likelihood, evaluated at its maximum.
The purely binary log-likelihood induced by model $m$ is readily obtained from Eq. ((ref)) and reads
where $a_{ij}$ is the generic entry of the empirical adjacency matrix and $p_{ij}$ is the model-dependent probability that node $i$ and node $j$ establish a connection. The `binary' AIC values (normalized by the yearly maximum, across models, for better visualization) are reported in Fig. (ref)a: the integrated exponential model outperforms the others, across all years, for both datasets. This result suggests that the information gained by including economic factors into the connection probabilities predicted by it does not affect the parsimony of its description, allowing it to perform better than the UBCM.
When, instead, the `full' log-likelihood is considered, reading
for integrated models and
for conditional models (see Fig. (ref)b), the conditional log-normal and gamma models compete, outperforming the other ones - although the performance of the first one in predicting the network statistics of interest, on the BACI dataset, was less remarkable than that of the competing models (see Fig. (ref)b).
We now complement the analysis of model performance, given in terms of realized likelihood, with an investigation of model `sensitivity', given in terms of the variability of the likelihood across network configurations sampled from the model. To this end, for each conditional model we build the so-called Shannon-Fisher plane Vignat2003, which is a technique that has acquired some popularity in the study of time-series. For instance it has been employed to understand ordinal patterns Rosso2012, quantify the degree of stochasticity Ravetti2014, classify financial stock markets Wang2018 and build indicators of economic efficiency Fernandes2021.
Within our context, we can use the Shannon-Fisher technique to project a given model onto a plane by assigning two coordinates to each connected dyad, i.e. to each pair of nodes $(i,j)$ with $a_{ij}=1$, where $a_{ij}$ is taken from the empirical adjacency matrix of the network. The $y$ coordinate in the plane is the Shannon entropy
which quantifies the degree of uncertainty encoded in the link weight. Note that, since the above entropy is constructed from a continuous pdf, it can attain negative values. This is a well-known problem that can be regularized by introducing the Kullback-Leibler divergence with respect to a continuous uniform pdf, however the result will only consist in an overall shift and rescaling of the $y$ coordinate that are inessential for our discussion below.
The $x$ coordinate in the plane is the Fisher Information Measure (FIM), defined as
and quantifying the (average) change in probability induced by small changes in the value of the link weight. Notice that the presence of the derivative requires that $q_{ij}(w|a_{ij}=1)$ is continuous throughout the domain of integration, and this is why we consider only conditional models with $a_{ij}=1$ so that there is no `jump' in the unconditional $q_{ij}(w)$ from $w=0$ to $w>0$. The expression $F_{ij}=\langle(H'_{ij})^2\rangle$ captures the `sensitivity' of the dyadic probability distribution with respect to small changes in the corresponding random variable. Note that this sensitivity is not captured by Shannon entropy, which is indifferent to any reordering of the values of the random variable, provided each value retains its probability.
In principle, two dyads with the same Shannon entropy can exhibit very different values of the FIM. This difference is captured by the Shannon-Fisher plane in terms of different positions along the $x$-axis. In general, since different connected dyads are described by a probability distribution with different parameters, scattering all connected dyads in the plane provides an overall representation of the model identified by $q_{ij}(w|a_{ij}=1)$. Different models are described by different probability distributions and hence have different projections in the Shannon-Fisher plane. In Appendix D we compute the explicit values of $S_{ij}$ and $F_{ij}$ for all the conditional models considered. Using these calculations we obtain the results shown in Fig. (ref) for an illustrative pair of datasets. We see that both the conditional exponential model and the conditional log-normal model follow a decreasing pattern. However, the conditional exponential model is, on average, characterized by a smaller FIM, i.e. smaller `sensitivity' to variations of the related random variable. On the other hand, the conditional Pareto model collapses onto a single point in the Shannon-Fisher plane while the conditional gamma model is characterized by a diverging FIM (because of the divergence of the first two negative moments, see Appendix D).
It is interesting to notice that, if we consider the sum of the $y$ values of all the connected dyads (a sort of `area under the curve') for a given model, we obtain the Shannon entropy for the entire weighted network, conditional on the empirical binary adjacency matrix $\mathbf{A}$:
(note that the dyadic entropy of $q(w_{ij}|a_{ij}=0)$ is zero, because if $a_{ij}=0$ then $w_{ij}=0$ deterministically). The above expression also coincides with minus the average likelihood of weighted network configurations sampled from the model (given the empirical binary structure), hence providing an average (inverse) `goodness of fit' of the weighted model. Similarly, summing the $x$ values of all the connected dyads gives an overall value of the FIM, hence the average change in likelihood of different weighted configurations sampled from the model. The results shown in Fig. (ref) therefore indicate that while different models (except the Pareto) are characterized by similar values of the overall entropy and goodness of fit, the conditional exponential model has minimum overall FIM, thereby producing the most stable outcome (in terms of likelihood of realized configurations) when used to sample weighted networks.
In a companion paper Marzio2022 the performance of discrete econometric models in reproducing the structural patterns of the WTW was compared with that of discrete, maximum-entropy ones. The analysis carried out there led to identify the zero-inflated Poisson model as the one performing best among the econometric models; still, it was also found to be largely disfavoured by information criteria such as AIC and BIC. This dilemma has been solved upon looking at a different class of statistical models, i.e. the physics-inspired ones: the latter have been found to outperform the purely econometric ones for reconstruction purposes, the reason lying in the higher accuracy achieved by them in estimating the topological structure of networks.
With this contribution, we extend the work carried out in Marzio2022 by, first, introducing models to infer the topology and the weights of (undirected, weighted) networks defined by continuous-valued data and, then, turning them into proper, econometric ones. In order to do so, we present a theoretical, physics-inspired framework based upon the constrained minimization of the KL divergence - hence, implementing the Minimum Discrimination Information Principle, that generalizes the Maximum-Entropy Principle - and capable of accommodating both integrated and conditional (continuous) models.
The main difference between the models belonging to these classes lies in the way the estimation of the topology is carried out; while conditional models disentangle the purely binary step from the (conditional) weighted one, integrated models do not, letting both topological and weighted constraints determine all relevant, structural features of a network. An example of integrated model is provided by the Enhanced Configuration Model (ECM), defined by constraints such as the degree and the strength sequences and described by a mixed Bernoulli-geometric Mastrandrea2014a,Mastrandrea2014b (also called Bose-Fermi Garlaschelli2009) distribution; examples of continuous, conditional models are provided by the $\text{CReM}_\text{A}$ and the $\text{CReM}_\text{B}$ Parisi2020. From a more econometric perspective, hurdle models are conditional in nature while zero-inflated models can be thought as integrated, the estimation steps being carried out by selecting a distribution out of a basket of available ones.
Our analysis leads to several conclusions: 1) constraining the entire degree sequence leads to a statistically significant improvement in the reconstruction accuracy of the WTW. In particular, the integrated exponential model, described by the Hamiltonian $H(\mathbf{W})=\sum_i\alpha_ik_i+\sum_{i<j}\beta_{ij}w_{ij}+\beta_0 W_1$, provides a very accurate, structural reconstruction while being favoured by information criteria: although it is defined by $N+1$ purely topological constraints, AIC reveals them as `irreducible', i.e. necessary to provide a satisfactory explanation of the network generating process; 2) when considering weighted quantities, the conditional gamma model is the one performing best (although it competes with the integrated exponential one in reproducing properties such as the weights, on some of the temporal snapshots covered by our datasets), according to information criteria. To be noticed, however, that if strengths are not explicitly constrained - jointly with the degrees - maximum-entropy models recover them only `in distribution' while failing to reproduce their exact values. The same consideration holds true for the weighted clustering coefficient.
Coming to comparing the models belonging to the classes considered in the present work, the two, best-performing ones are the integrated exponential model and the conditional gamma model, i.e. the ones constraining the total weight (although the conditional exponential model constrains the total weight as well, it is outperformed by the conditional gamma one within the class of conditional model): hence, $W_1$ seems to constitute a somehow fundamental quantity to be necessarily accounted for in order to achieve a good reconstruction accuracy. From an economic point of view, the parameter $\beta_0$ constraining the total weight can be interpreted as a sort of `shadow price' to be paid by everyone to exchange goods.
Additional information is provided by our analysis of the Shannon-Fisher plane, which combines Shannon entropy, i.e. the (inverse) likelihood of a model, with the Fisher Information Measure, i.e. the average variability of the likelihood itself across different sampled configurations. The conditional exponential model turns out the be the least variable in likelihood, hence the most stable. It is worth noticing at this point that our maximum-entropy approach is formulated for canonical ensembles, i.e. for `soft constraints', which implies that different realizations of the network have fluctuating values of the weighted sufficient statistics. These fluctuations are the origin of the FIM. By contrast, if we were to formulate microcanonical models with `hard constraints', then the sufficient statistics would not fluctuate and the overall FIM would be zero. Therefore the Shannon-Fisher plane shows that, among the canonical models considered here, the conditional exponential is the closest to the `least soft' extreme, while the conditional gamma is at the opposite `softest' extreme where the FIM diverges. As a question left for future research, it would be interesting to relate the behaviour of the FIM to the phenomenon of inequivalence of canonical and microcanonical ensembles of networks squartini2015breaking.
Overall, we believe the framework proposed in this contribution to have the potential of reconciling the approach adopted by network scientists for reconstructing economic networks, and focusing on the purely structural aspects of a network formation, with the approach characterizing econometrics, tailored to inform these same models with macro-economic quantities - in all cases considered here, purely bilateral ones such as the GDPs and the geographic distances. From an operative point of view, our (classes of) models combine the pros of both approaches: the importance of purely structural information (highlighted by physics-inspired models) can be accounted for by constraining the entire degree sequence; on top of that, a second step is needed to estimate a network weighted structure. Although the information provided by the total weight cannot be discarded without affecting the overall performance of a model, such an estimation can rests upon econometric considerations driving the reparametrization of otherwise purely structural models.
As an additional result, we release a Python package named `DyGyS - DYadic GravitY regression models with Soft constraints' and containing routines to implement all models considered in the present work as well as those considered in the companion paper Marzio2022. The package is available at the following URL: \href{https://github.com/MarsMDK/DyGyS}{https://github.com/MarsMDK/DyGyS}.
This work is supported by the European Union – Horizon 2020 Program under the scheme “INFRAIA-01-2018-2019 – Integrating Activities for Advanced Communities”, Grant Agreement n.871042, “SoBigData++: European Integrated Infrastructure for Social Mining and Big Data Analytics”(\href{http://www.sobigdata.eu}{http://www.sobigdata.eu}). This work has been also supported by the project `Network analysis of economic and financial resilience', Italian DM n. 289, 25-03-2021 (PRO3 Scuole) CUP D67G22000130001. DG acknowledges support from the Dutch Econophysics Foundation (Stichting Econophysics, Leiden, the Netherlands) and the Netherlands Organization for Scientific Research (NWO/OCW). MDV and DG acknowledge support from the `Programma di Attivit\`a Integrata' (PAI) project `Prosociality, Cognition and Peer Effects' (Pro.Co.P.E.), funded by IMT School for Advanced Studies.