# Reconciling econometrics with continuous maximum-entropy network models

Marzio Di Vece<EMAIL_ADDRESS>IMT School for Advanced Studies,
P.zza San Francesco 19, 55100 Lucca (Italy)  Diego Garlaschelli  IMT School
for Advanced Studies, P.zza San Francesco 19, 55100 Lucca (Italy)  Lorentz
Institute for Theoretical Physics (LION), Leiden University, Niels Bohrweg 2,
2333CA Leiden (The Netherlands)  Tiziano Squartini  IMT School for Advanced
Studies, P.zza San Francesco 19, 55100 Lucca (Italy)  Institute for Advanced
Study (IAS), University of Amsterdam, Oude Turfmarkt 145, 1012GC Amsterdam
(The Netherlands)

###### Abstract

In the study of economic networks, econometric approaches interpret the
traditional Gravity Model specification as the expected link weight coming
from a probability distribution whose functional form can be chosen
arbitrarily, while statistical-physics approaches construct maximum-entropy
distributions of weighted graphs, constrained to satisfy a given set of
measurable network properties. In a recent, companion paper, we integrated the
two approaches and applied them to the World Trade Web, i.e. the network of
international trade among world countries. While the companion paper dealt
only with discrete-valued link weights, the present paper extends the
theoretical framework to continuous-valued link weights. In particular, we
construct two broad classes of maximum-entropy models, namely the integrated
and the conditional ones, defined by different criteria to derive and combine
the probabilistic rules for placing links and loading them with weights. In
the integrated models, both rules follow from a single, constrained
optimization of the continuous Kullback-Leibler divergence; in the conditional
models, the two rules are disentangled and the functional form of the weight
distribution follows from a conditional, optimization procedure. After
deriving the general functional form of the two classes, we turn each of them
into a proper family of econometric models via a suitable identification of
the econometric function relating the corresponding, expected link weights to
macroeconomic factors. After testing the two classes of models on World Trade
Web data, we discuss their strengths and weaknesses.

###### pacs:

89.75.Fb; 02.50.Tt; 89.65.Gh

##  I Introduction

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) [1, 2, 3, 4, 5, 6,
7, 12, 8, 9, 10, 11].

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) [13] 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

$\langle
w_{ij}\rangle_{\text{GM}}=f(\omega_{i},\omega_{j},d_{ij}|\underline{\psi})=\tau\omega_{i}^{\beta_{1}}\omega_{j}^{\beta_{2}}d_{ij}^{\gamma}$
(1)

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 1 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 [14] and the
maximization of utility functions constrained to satisfy requirements on the
rate of information acquisition [15].

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 [16]. 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. 1 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. 1
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. 1 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 [1, 2, 3, 4, 5, 6, 7, 12, 8, 9, 10, 11].

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) [17]. 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. 1. 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 [19] and _hurdle_ models [20]. 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 [18];
moreover, only discrete distributions (specifically, either the Poisson or the
negative-binomial one [19]) 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 [20]: 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) [37, 38, 39] applied to network
ensembles [21, 22], prescribing to maximize Shannon entropy [23, 24] 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 [25] 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 [5, 6, 26, 27, 28, 29, 30, 31, 32, 18]. The entire framework
can also accommodate possibly degenerate, discrete-valued, single- or multi-
edges [33].

Although maximum-entropy models have been also studied from an economic
perspective (see [34] 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 [5, 6, 29, 30, 31, 32, 18] or by
introducing network-related statistics into otherwise purely econometric
models [35]. 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 [18], 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) [36] 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 [21] or weighted ERGs [25]
and endows them with macroeconomic factors replacing certain model parameters.

The rest of the paper is organized as follows: in Sec. II, after introducing
the basic quantities, we derive the class of conditional models; in Sec. III
we derive the class of integrated models; in Sec. IV we apply all models to
the analysis of WTW data; in Sec. V we discuss the results and provide our
concluding remarks.

##  II Conditional Models

Discrete maximum-entropy models can be derived by performing a constrained
maximization of Shannon entropy [37, 38, 39]. However, unlike the companion
paper [18], 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

$D_{\text{KL}}(Q||R)=\int_{\mathbb{W}}Q(\mathbf{W})\ln\frac{Q(\mathbf{W})}{R(\mathbf{W})}d\mathbf{W}$
(2)

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 [40] 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 [23] and machine learning [41],
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

$Q(\mathbf{W})=P(\mathbf{A})Q(\mathbf{W}|\mathbf{A}),$ (3)

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

$D_{\text{KL}}(Q||R)=S(Q,R)-S(P)-S(Q_{\bot}|P)$ (4)

where

$S(P)=-\sum_{\mathbf{A}\in\mathbb{A}}P(\mathbf{A})\ln P(\mathbf{A})$ (5)

is the _Shannon entropy_ of the probability distribution describing the binary
projection of the network structure,

$S(Q_{\bot}|P)=-\sum_{\mathbf{A}\in\mathbb{A}}P(\mathbf{A})\int_{\mathbb{W}_{\mathbf{A}}}Q(\mathbf{W}|\mathbf{A})\ln
Q(\mathbf{W}|\mathbf{A})d\mathbf{W}$ (6)

is the _conditional Shannon entropy_ of the probability distribution of the
weighted network structure given the binary projection and

$S(Q,R)=-\sum_{\mathbf{A}\in\mathbb{A}}P(\mathbf{A})\int_{\mathbb{W}_{\mathbf{A}}}Q(\mathbf{W}|\mathbf{A})\ln
R(\mathbf{W})d\mathbf{W}$ (7)

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

$R(\mathbf{W})=T(\mathbf{A})R(\mathbf{W}|\mathbf{A})$ (8)

an expression that allows us to write $S(Q,R)$ as

$\displaystyle S(Q,R)=$
$\displaystyle-\sum_{\mathbf{A}\in\mathbb{A}}P(\mathbf{A})\ln T(\mathbf{A})$
$\displaystyle-\sum_{\mathbf{A}\in\mathbb{A}}P(\mathbf{A})\int_{\mathbb{W}_{\mathbf{A}}}Q(\mathbf{W}|\mathbf{A})\ln
R(\mathbf{W}|\mathbf{A})d\mathbf{W}$ (9)

which, in turn, allows the KL divergence to be rewritten as

$D_{\text{KL}}(Q||R)=-D_{\text{KL}}(P||T)-D_{\text{KL}}(Q_{\bot}||R_{\bot})$
(10)

i.e. as a sum of two terms, one of which involves conditional distributions;
specifically,

$\displaystyle D_{\text{KL}}(P||T)$
$\displaystyle=-\sum_{\mathbf{A}\in\mathbb{A}}P(\mathbf{A})\ln\frac{P(\mathbf{A})}{T(\mathbf{A})},$
(11) $\displaystyle D_{\text{KL}}(Q_{\bot}||R_{\bot})$
$\displaystyle=-\sum_{\mathbf{A}\in\mathbb{A}}P(\mathbf{A})\int_{\mathbb{W}_{\mathbf{A}}}Q(\mathbf{W}|\mathbf{A})\ln\frac{Q(\mathbf{W}|\mathbf{A})}{R(\mathbf{W}|\mathbf{A})}d\mathbf{W}$
(12)

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

$D_{\text{KL}}(Q||R)=-S(P)-S(Q_{\bot}|P)$ (13)

with

$\displaystyle S(P)$
$\displaystyle=-\sum_{\mathbf{A}\in\mathbb{A}}P(\mathbf{A})\ln P(\mathbf{A}),$
(14) $\displaystyle S(Q_{\bot}|P)$
$\displaystyle=-\sum_{\mathbf{A}\in\mathbb{A}}P(\mathbf{A})\int_{\mathbb{W}_{\mathbf{A}}}Q(\mathbf{W}|\mathbf{A})\ln
Q(\mathbf{W}|\mathbf{A})d\mathbf{W}.$ (15)

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.

###  II.1 Choosing the binary constraints

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

$S(P)=-\sum_{\mathbf{A}\in\mathbb{A}}P(\mathbf{A})\ln P(\mathbf{A})$ (16)

leading to a probability mass function reading

$P(\mathbf{A})=\frac{e^{-H(\mathbf{A})}}{\sum_{\mathbb{A}}e^{-H(\mathbf{A})}}$
(17)

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

$P(\mathbf{A})=\prod_{i<j}p_{ij}^{a_{ij}}(1-p_{ij})^{1-a_{ij}}$ (18)

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

$p_{ij}^{\text{UBCM}}=\frac{x_{i}x_{j}}{1+x_{i}x_{j}}$ (19)

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

$p_{ij}^{\text{FM}}=\frac{\delta\omega_{i}\omega_{j}}{1+\delta\omega_{i}\omega_{j}}$
(20)

and requiring the estimation of a global parameter only, i.e. $\delta$. The FM
represents a particular case of the logit model [42], being defined by a
vector of external properties (the ‘fitnesses’) that replace the information
provided by some kind of (otherwise) purely structural properties [43, 5]: in
fact, Eq. 20 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 [26, 6] as accurately as
the UBCM, although requiring much less information.

In what follows, we will consider both the UBCM and the FM specifications.

###  II.2 Choosing the weighted constraints

The constrained maximization of $S(Q_{\bot}|P)$ proceeds by specifying the
following set of weighted constraints

$\displaystyle 1$
$\displaystyle=\int_{\mathbb{W}_{\mathbf{A}}}P(\mathbf{W}|\mathbf{A})d\mathbf{W},\>\forall\>\mathbf{A}\in\mathbb{A}$
(21) $\displaystyle\langle C_{\alpha}\rangle$
$\displaystyle=\sum_{\mathbf{A}\in\mathbb{A}}P(\mathbf{A})\int_{\mathbb{W}_{\mathbf{A}}}Q(\mathbf{W}|\mathbf{A})C_{\alpha}(\mathbf{W})d\mathbf{W},\>\forall\>\alpha$
(22)

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

$Q(\mathbf{W}|\mathbf{A})=\begin{cases}\frac{e^{-H(\mathbf{W})}}{Z_{\mathbf{A}}},&\>\mathbf{W}\in\mathbb{W}_{\mathbf{A}}\\\
0,&\>\mathbf{W}\notin\mathbb{W}_{\mathbf{A}}\end{cases}$ (23)

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

$\displaystyle
H(\mathbf{W})=\sum_{i<j}f(w_{ij}|\beta_{ij},\xi_{ij},\gamma_{ij})$ (24)

with the Lagrange multipliers $\left(\beta_{ij},\xi_{ij},\gamma_{ij}\right)$
satisfying the following requirements:

  * •

$\beta_{ij}\equiv\beta_{0}+\beta_{ij}$, where $\beta_{0}$ is the Lagrange
multiplier associated with the total weight $\sum_{i<j}w_{ij}\equiv W_{1}$ and
$\beta_{ij}$ encodes the dependence on purely econometric quantities;

  * •

$\xi_{ij}$ will be kept either in its dyadic form, to constrain the logarithm
of each weight, or in its global form, $\xi_{ij}\equiv\xi_{0}$, to constrain
the sum of the logarithms of the weights, i.e. $\sum_{i<j}\ln(w_{ij})\equiv
W_{2}$;

  * •

$\gamma_{ij}\equiv\gamma_{0}$ plays the role of the Lagrange multiplier
associated with (a function of) the total variance of the logarithms of the
weights, i.e. $\sum_{i<j}\ln^{2}(w_{ij})\equiv W_{3}$.

####  II.2.1 Conditional exponential model.

Let us start by considering the simplest, conditional model, defined by the
positions $\gamma_{ij}=\xi_{ij}=0$ and inducing the Hamiltonian

$\displaystyle H(\mathbf{W})$
$\displaystyle=\sum_{i<j}f(w_{ij}|\beta_{0}+\beta_{ij})$
$\displaystyle=\sum_{i<j}(\beta_{0}+\beta_{ij})w_{ij};$ (25)

inserting the expression above into Eq. 23 leads to the distribution

$\displaystyle Q(\mathbf{W}|\mathbf{A})$
$\displaystyle=\prod_{i<j}q_{ij}(w_{ij}|a_{ij})$
$\displaystyle=\prod_{i<j}\frac{e^{-(\beta_{0}+\beta_{ij})w_{ij}}}{\zeta_{ij}}$
$\displaystyle=\prod_{i<j}(\beta_{0}+\beta_{ij})e^{-(\beta_{0}+\beta_{ij})w_{ij}}$
(26)

and each node pair-specific distribution induces a (conditional) expected
weight reading

$\langle w_{ij}|a_{ij}=1\rangle=\frac{1}{\beta_{0}+\beta_{ij}}.$ (27)

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.

$\beta_{ij}\equiv z_{ij}^{-1},$ (28)

a position that turns Eq. 27 into

$\langle
w_{ij}|a_{ij}=1\rangle=\frac{1}{\beta_{0}+z_{ij}^{-1}}=\frac{z_{ij}}{1+\beta_{0}z_{ij}};$
(29)

notice that the only structural constraint is, now, represented by the total
weight (see also Appendix A).

####  II.2.2 Conditional gamma model.

Let us, now, consider a different Hamiltonian, constraining each weight, their
total sum and the sum of their logarithms, i.e.

$\displaystyle H(\mathbf{W})$
$\displaystyle=\sum_{i<j}f(w_{ij}|\beta_{0}+\beta_{ij},\xi_{0})$
$\displaystyle=\sum_{i<j}[(\beta_{0}+\beta_{ij})w_{ij}+\xi_{0}\ln(w_{ij})]$
$\displaystyle=\sum_{i<j}\beta_{ij}w_{ij}+\beta_{0}W_{1}+\xi_{0}W_{2};$ (30)

it induces the distribution reading

$\displaystyle Q(\mathbf{W}|\mathbf{A})$
$\displaystyle=\prod_{i<j}q_{ij}(w_{ij}|a_{ij})$
$\displaystyle=\prod_{i<j}\dfrac{e^{-(\beta_{0}+\beta_{ij})w_{ij}-\xi_{0}\ln(w_{ij})}}{\zeta_{ij}}$
$\displaystyle=\prod_{i<j}\dfrac{e^{-(\beta_{0}+\beta_{ij})w_{ij}}w_{ij}^{-\xi_{0}}}{\zeta_{ij}}$
$\displaystyle=\prod_{i<j}\dfrac{(\beta_{0}+\beta_{ij})^{1-\xi_{0}}}{\Gamma(1-\xi_{0})}e^{-(\beta_{0}+\beta_{ij})w_{ij}}w_{ij}^{-\xi_{0}};$
(31)

each node pair-specific distribution is characterized by a (conditional)
expected weight reading

$\langle w_{ij}|a_{ij}=1\rangle=\frac{1-\xi_{0}}{\beta_{0}+\beta_{ij}}$ (32)

and by a (conditional) expected logarithmic weight reading

$\langle\ln(w_{ij})|a_{ij}=1\rangle=\psi(1-\xi_{0})-\ln(\beta_{0}+\beta_{ij})$
(33)

where the function $\psi(x)=\Gamma^{\prime}(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.

$\beta_{ij}\equiv z_{ij}^{-1}$ (34)

a position that, in turn, leads to the expressions

$\langle
w_{ij}|a_{ij}=1\rangle=\frac{1-\xi_{0}}{\beta_{0}+z_{ij}^{-1}}=\frac{(1-\xi_{0})z_{ij}}{1+\beta_{0}z_{ij}}$
(35)

(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

$\langle\ln(w_{ij})|a_{ij}=1\rangle=\psi(1-\xi_{0})-\ln(\beta_{0}+z_{ij}^{-1})$
(36)

(see also Appendix A).

####  II.2.3 Conditional Pareto model.

Constraining a slightly more complex function of the weights, i.e. their
logarithm, leads to the Hamiltonian

$\displaystyle H(\mathbf{W})$ $\displaystyle=\sum_{i<j}f(w_{ij}|\xi_{ij})$
$\displaystyle=\sum_{i<j}\xi_{ij}\ln(w_{ij})$ (37)

which, in turn, induces the distribution

$\displaystyle Q(\mathbf{W}|\mathbf{A})$
$\displaystyle=\prod_{i<j}q_{ij}(w_{ij}|a_{ij})$
$\displaystyle=\prod_{i<j}\dfrac{e^{-\xi_{ij}\ln(w_{ij})}}{\zeta_{ij}}$
$\displaystyle=\prod_{i<j}\dfrac{w_{ij}^{-\xi_{ij}}}{\zeta_{ij}}$
$\displaystyle=\prod_{i<j}\dfrac{(\xi_{ij}-1)}{m_{ij}^{1-\xi_{ij}}}w_{ij}^{-\xi_{ij}}$
(38)

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

$\langle
w_{ij}|a_{ij}=1\rangle=\left(\dfrac{\xi_{ij}-1}{\xi_{ij}-2}\right)m_{ij}.$
(39)

Such a model can be turned into a proper, econometric one by considering the
positions

$\displaystyle\xi_{ij}-2$ $\displaystyle\equiv z_{ij}^{-1},$ $\displaystyle
m_{ij}$ $\displaystyle\equiv w_{min}$ (40)

ensuring that the expected weights are monotonically increasing functions of
the gravity specification and leading to the expression

$\langle w_{ij}|a_{ij}=1\rangle=(1+z_{ij})w_{min}$ (41)

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
[44]; 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.

####  II.2.4 Conditional log-normal model.

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

$\displaystyle H(\mathbf{W})$
$\displaystyle=\sum_{i<j}f(w_{ij}|\gamma_{0},\xi_{ij})$
$\displaystyle=\sum_{i<j}[\xi_{ij}\ln(w_{ij})+\gamma_{0}\ln^{2}(w_{ij})]$
$\displaystyle=\sum_{i<j}\xi_{ij}\ln(w_{ij})+\gamma_{0}W_{3};$ (42)

the Hamiltonian above induces a distribution reading

$\displaystyle Q(\mathbf{W}|\mathbf{A})$
$\displaystyle=\prod_{i<j}q_{ij}(w_{ij}|a_{ij})$
$\displaystyle=\prod_{i<j}\dfrac{e^{-\xi_{ij}\ln(w_{ij})-\gamma_{0}\ln^{2}(w_{ij})}}{\zeta_{ij}}$
$\displaystyle=\prod_{i<j}\frac{e^{-\xi_{ij}\ln(w_{ij})-\gamma_{0}\ln^{2}(w_{ij})}}{\sqrt{\frac{\pi}{\gamma_{0}}}e^{\frac{(\xi_{ij}-1)^{2}}{4\gamma_{0}}}};$
(43)

each node pair-specific distribution is characterized by a (conditional)
expected weight reading

$\langle w_{ij}|a_{ij}=1\rangle=e^{\frac{3-2\xi_{ij}}{4\gamma_{0}}},$ (44)

by a (conditional) expected, logarithmic weight reading

$\langle\ln(w_{ij})|a_{ij}=1\rangle=\frac{1-\xi_{ij}}{2\gamma_{0}}$ (45)

and by a (conditional) logarithmic weight whose squared expectation reads

$\langle\ln^{2}(w_{ij})|a_{ij}=1\rangle=\frac{2\gamma_{0}+(1-\xi_{ij})^{2}}{4\gamma_{0}^{2}}.$
(46)

Such a model can be turned into a proper, econometric one by considering the
position

$1-\xi_{ij}\equiv\ln(z_{ij})$ (47)

ensuring that the expected weights are monotonically increasing functions of
the gravity specification and leading to the expressions

$\displaystyle\langle w_{ij}|a_{ij}=1\rangle$
$\displaystyle=e^{\frac{1+2\ln(z_{ij})}{4\gamma_{0}}},$ (48)
$\displaystyle\langle\ln(w_{ij})|a_{ij}=1\rangle$
$\displaystyle=\frac{\ln(z_{ij})}{2\gamma_{0}},$ (49)
$\displaystyle\langle\ln^{2}(w_{ij})|a_{ij}=1\rangle$
$\displaystyle=\frac{2\gamma_{0}+\ln^{2}(z_{ij})}{4\gamma_{0}^{2}}$ (50)

(see also Appendix A).

##  III Integrated Models

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

$\displaystyle 1$ $\displaystyle=\int_{\mathbb{W}}Q(\mathbf{W})d\mathbf{W},$
(51) $\displaystyle\langle C_{\alpha}\rangle$
$\displaystyle=\int_{\mathbb{W}}Q(\mathbf{W})C_{\alpha}(\mathbf{W})d\mathbf{W},\>\forall\>\alpha$
(52)

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

$Q(\mathbf{W})=\frac{R(\mathbf{W})e^{-H(\mathbf{W})}}{\int_{\mathbb{W}}R(\mathbf{W})e^{-H(\mathbf{W})}d\mathbf{W}}$
(53)

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

$Q(\mathbf{W})=\frac{e^{-H(\mathbf{W})}}{\int_{\mathbb{W}}e^{-H(\mathbf{W})}d\mathbf{W}};$
(54)

notice that the result above could have been also derived by carrying out a
constrained minimization of

$D_{\text{KL}}(Q||R)=\int_{\mathbb{W}}Q(\mathbf{W})\ln
Q(\mathbf{W})d\mathbf{W}\equiv-S(Q)$ (55)

i.e. of (minus) the functional named _differential entropy_ into which the KL
divergence ‘degenerates’ in case completely uninformative priors are
considered.

###  III.1 Choosing the constraints

Let us, now, specify the functional form of the constraints. In what follows,
we will deal with a specific instance of the generic Hamiltonian

$\displaystyle H(\mathbf{W})=\sum_{i<j}f(w_{ij}|\alpha_{ij},\beta_{ij});$ (56)

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.

$\displaystyle H(\mathbf{W})$
$\displaystyle=\sum_{i<j}f(w_{ij}|\alpha_{i}+\alpha_{j},\beta_{0}+\beta_{ij})$
$\displaystyle=\sum_{i<j}[(\alpha_{i}+\alpha_{j})a_{ij}+(\beta_{0}+\beta_{ij})w_{ij}]$
$\displaystyle=\sum_{i}\alpha_{i}k_{i}+\sum_{i<j}(\beta_{0}+\beta_{ij})w_{ij}$
$\displaystyle=\sum_{i}\alpha_{i}k_{i}+\sum_{i<j}\beta_{ij}w_{ij}+\beta_{0}W_{1}$

that leads to the distribution

$\displaystyle Q(\mathbf{W})$ $\displaystyle=\prod_{i<j}q_{ij}(w_{ij})$
$\displaystyle=\prod_{i<j}\dfrac{(x_{i}x_{j})^{a_{ij}}e^{-(\beta_{0}+\beta_{ij})w_{ij}}}{Z_{ij}}$
$\displaystyle=\prod_{i<j}\dfrac{(x_{i}x_{j})^{a_{ij}}e^{-(\beta_{0}+\beta_{ij})w_{ij}}}{1+x_{i}x_{j}(\beta_{0}+\beta_{ij})^{-1}}$
(57)

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

$\displaystyle p_{ij}=1-q_{ij}(0)$
$\displaystyle=\dfrac{x_{i}x_{j}(\beta_{0}+\beta_{ij})^{-1}}{1+x_{i}x_{j}(\beta_{0}+\beta_{ij})^{-1}}$
$\displaystyle=\dfrac{x_{i}x_{j}\zeta_{ij}}{1+x_{i}x_{j}\zeta_{ij}};$ (58)

besides, the corresponding expected weight reads

$\langle w_{ij}\rangle=\frac{p_{ij}}{\beta_{0}+\beta_{ij}}.$ (59)

Equation III.1 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.

$\beta_{ij}\equiv z_{ij}^{-1}$ (60)

a position that turns Eq. 59 into

$\langle w_{ij}\rangle=\frac{p_{ij}}{\beta_{0}+z_{ij}^{-1}}$ (61)

where

$p_{ij}=\frac{x_{i}x_{j}}{x_{i}x_{j}+\beta_{0}+z_{ij}^{-1}}$ (62)

(see also Appendix B).

##  IV Results

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 [45]) and the BACI one (covering 11 years, from 2007 to 2017
[46]). 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}$: as usual, a real number
$u_{ij}$ is drawn from the uniform distribution whose support coincides with
the unit interval, i.e. $U[0,1]$, and compared with $p_{ij}$; if $u_{ij}\leq
p_{ij}$, $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. According to it, it is possible to generate
random variables by considering the complementary cumulative distribution
$F(v_{ij})$, defined as

$F(v_{ij})=\int_{0}^{v_{ij}}q(w_{ij}|a_{ij}=1)dw_{ij};$ (63)

since $F(v_{ij})$ behaves like a random variable, say $v$, uniformly
distributed between 0 and 1, the equation $F(v_{ij})=v$ can be inverted to
obtain $v_{ij}$, i.e. the weight to be assigned to the pair $i<j$. Each
ensemble is constituted by a number of configurations amounting at $10^{4}$;
the error accompanying the estimate of any quantity of interest is quantified
via the confidence intervals (CI) induced by the ensemble distribution of the
former one.

Dataset |  Model |  $f^{k_{i}}$ |  $f^{k_{nn}}$ |  $f^{c_{i}}$
---|---|---|---|---
Gleditsch |  I-Exp |  $1$ |  $1$ |  $1$
Gleditsch |  UBCM |  $1$ |  $1$ |  $1$
Gleditsch |  FM |  $1$ |  $1$ |  $1$

Dataset |  Model |  $f^{k_{i}}$ |  $f^{k_{nn}}$ |  $f^{c_{i}}$
---|---|---|---|---
BACI |  I-Exp |  $1$ |  $0.09$ |  $0.09$
BACI |  UBCM |  $1$ |  $0.09$ |  $0.09$
BACI |  FM |  $0.09$ |  $0.09$ |  $0.09$

Dataset |  Model |  $f^{s_{i}}$ |  $f^{s_{nn}}$ |  $f^{c_{w}}$
---|---|---|---|---
Gleditsch |  I-Exp |  $1$ |  $1$ |  $0.18$
Gleditsch |  C-Exp |  $1$ |  $1$ |  $0.18$
Gleditsch |  C-Pareto |  $0$ |  $0$ |  $0$
Gleditsch |  C-Gamma |  $1$ |  $1$ |  $0.27$
Gleditsch |  C-Lognormal |  $1$ |  $0$ |  $0$

Dataset |  Model |  $f^{s_{i}}$ |  $f^{s_{nn}}$ |  $f^{c_{w}}$
---|---|---|---|---
BACI |  I-Exp |  $1$ |  $1$ |  $1$
BACI |  C-Exp |  $1$ |  $1$ |  $1$
BACI |  C-Pareto |  $0$ |  $0$ |  $0$
BACI |  C-Gamma |  $1$ |  $1$ |  $1$
BACI |  C-Lognormal |  $0$ |  $0$ |  $0$

Table 1: Compatibility between the distributions of the expected values of the
statistics output by the models considered in the present work and their
empirical counterparts for the Gleditsch (left) and the BACI (right) dataset.
A large value of $f^{s}_{m}$ indicates a large percentage of years for which
the distribution of the network statistic predicted by model $m$ is compatible
with the empirical one. While all models seem to perform quite well in
reproducing the binary statistics on the Gleditsch dataset, this is no longer
true when considering the BACI dataset, on which the UBCM outperform the FM -
a result that leads us to prefer the former as the ‘first step-algorithm’ of
our conditional models. For what concerns the set of weighted statistics, the
models constraining $W_{1}$ (i.e. the integrated exponential model, the
conditional exponential model and the conditional gamma model) clearly
outperform the others.

###  IV.1 Model selection via statistical indicators

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

$k_{i}=\sum_{j(\neq i)=1}^{N}a_{ij},\quad\forall\>i$ (64)

which gives information about the tendency of single nodes to link with other
trade partners, the average nearest neighbors degree

$k^{nn}_{i}=\dfrac{\sum_{j(\neq i)=1}^{N}a_{ij}k_{j}}{k_{i}},\quad\forall\>i$
(65)

which gives information about the degree correlations, the clustering
coefficient

$c_{i}=\dfrac{\sum_{j(\neq i)=1}^{N}\sum_{k(\neq
i,j)=1}^{N}a_{ij}a_{jk}a_{ki}}{k_{i}(k_{i}-1)},\quad\forall\>i$ (66)

which counts the percentage of node partners that are also partners
themselves. For what concerns the weighted statistics, we have considered the
strength sequence

(a) Reconstruction accuracy for the binary network statistics.

(b) Reconstruction accuracy for the weighted network statistics.

Figure 1: Reconstruction accuracy for the statistics of interest, calculated
as the percentage of node-specific values of the statistics $s$ falling within
the CI, induced by model $m$, at the significance level of $5\%$; yearly
percentages are, then, averaged. The whiskers represent the $2.5$ and $97.5$
percentiles of each $\text{RA}^{s}_{m}$ distribution, across different years.
Overall, all models perform quite well in reproducing the degrees, with the
only exception of the FM - whence our choice of employing the UBCM as the
‘first step-algorithm’ of our conditional models. For what concerns higher-
order, binary statistics, both the UBCM and the I-Exp perform quite well while
the performance of the FM is much poorer - a result that holds true for both
the Gleditsch (left panels) and the BACI (right panels) dataset. For what
concerns the weighted statistics, only the average nearest neighbors strength
is satisfactorily recovered - however, only by the (integrated and
conditional) exponential models and the conditional gamma one.

$s_{i}=\sum_{j(\neq i)=1}^{N}w_{ij},\quad\forall\>i$ (67)

which gives information about the trade flow of a country, the average nearest
neighbors strength

$s^{nn}_{i}=\dfrac{\sum_{j(\neq i)=1}^{N}a_{ij}s_{j}}{k_{i}},\quad\forall\>i$
(68)

which gives information about the strength correlations, the weighted
clustering coefficient

$c^{w}_{i}=\dfrac{\sum_{j(\neq i)=1}^{N}\sum_{k(\neq
i,j)=1}^{N}w_{ij}w_{jk}w_{ki}}{k_{i}(k_{i}-1)},\quad\forall\>i$ (69)

that weighs the closed triangular patterns node $i$ establishes with other
trade partners.

Figure 2: Reconstruction accuracy $\text{RA}^{w}_{m}$ for the Gleditsch
dataset (left panel) and the BACI dataset (right panel). 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 while
the best-performing model on the BACI dataset is the conditional gamma model.

####  IV.1.1 KS compatibility frequency

Table 1 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 1, 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.

####  IV.1.2 Reconstruction accuracy

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 1 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. 1a), 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. 1b), 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\%$ [48]. 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. 2 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.

####  IV.1.3 Confusion matrix

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 2, 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.

###  IV.2 Model selection via statistical tests

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.

###  IV.3 Model selection via information criteria

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) [49],
reading

$\text{AIC}_{m}=2k-2\mathcal{L}_{m}$ (70)

where $k$ is the number of free parameters of the model and $\mathcal{L}_{m}$
is its log-likelihood, evaluated at its maximum.

Dataset |  $\Delta_{\text{TPR}}$ |  $\Delta_{\text{SPC}}$ |  $\Delta_{\text{PPV}}$ |  $\Delta_{\text{ACC}}$
---|---|---|---|---
Gleditsch 90 |  $0.017$ |  $0.026$ |  $0.017$ |  $0.020$
Gleditsch 91 |  $0.016$ |  $0.020$ |  $0.016$ |  $0.018$
Gleditsch 92 |  $0.015$ |  $0.019$ |  $0.015$ |  $0.017$
Gleditsch 93 |  $0.015$ |  $0.019$ |  $0.015$ |  $0.017$
Gleditsch 94 |  $0.013$ |  $0.018$ |  $0.013$ |  $0.015$
Gleditsch 95 |  $0.013$ |  $0.019$ |  $0.013$ |  $0.016$
Gleditsch 96 |  $0.012$ |  $0.019$ |  $0.012$ |  $0.015$
Gleditsch 97 |  $0.013$ |  $0.022$ |  $0.013$ |  $0.016$
Gleditsch 98 |  $0.013$ |  $0.022$ |  $0.013$ |  $0.016$
Gleditsch 99 |  $0.013$ |  $0.023$ |  $0.014$ |  $0.017$
Gleditsch 00 |  $0.014$ |  $0.023$ |  $0.014$ |  $0.017$

Dataset |  $\Delta_{\text{TPR}}$ |  $\Delta_{\text{SPC}}$ |  $\Delta_{\text{PPV}}$ |  $\Delta_{\text{ACC}}$
---|---|---|---|---
BACI 07 |  $0.007$ |  $0.037$ |  $0.007$ |  $0.012$
BACI 08 |  $0.006$ |  $0.035$ |  $0.006$ |  $0.010$
BACI 09 |  $0.005$ |  $0.031$ |  $0.005$ |  $0.009$
BACI 10 |  $0.005$ |  $0.032$ |  $0.005$ |  $0.009$
BACI 11 |  $0.005$ |  $0.029$ |  $0.006$ |  $0.008$
BACI 12 |  $0.005$ |  $0.033$ |  $0.005$ |  $0.009$
BACI 13 |  $0.005$ |  $0.034$ |  $0.005$ |  $0.009$
BACI 14 |  $0.005$ |  $0.031$ |  $0.005$ |  $0.008$
BACI 15 |  $0.005$ |  $0.028$ |  $0.005$ |  $0.008$
BACI 16 |  $0.003$ |  $0.021$ |  $0.004$ |  $0.006$
BACI 17 |  $0.004$ |  $0.028$ |  $0.004$ |  $0.007$

Table 2: Increments of the four indicators composing the confusion matrix,
i.e. the true positive rate (TPR), the specificity (SPC), the positive
predicted value (PPV) and the accuracy (ACC) when passing from the UBCM to the
integrated exponential model for the Gleditsch (left table) and the BACI
(right table) datasets. All increments are significant at the $1\%$ level,
according to the (non-parametric) Wilcoxon rank-sum test on the ensemble
distributions of the statistics to compare.

The purely binary log-likelihood induced by model $m$ is readily obtained from
Eq. (18) and reads

$\mathcal{L}^{(b)}_{m}=\ln
P(\mathbf{A})=\sum_{i<j}\left[a_{ij}\ln(p_{ij})+(1-a_{ij})\ln(1-p_{ij})\right]$
(71)

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. 3a: 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

$\mathcal{L}^{(f)}_{\text{I-}m}=\ln
Q(\mathbf{W})=\sum_{i<j}\ln(q_{ij}(w_{ij}))$ (72)

for integrated models and

$\displaystyle\mathcal{L}^{(f)}_{\text{C-}m}=$ $\displaystyle\ln
P(\mathbf{A})+\ln Q(\mathbf{W}|\mathbf{A})$ $\displaystyle=$
$\displaystyle\sum_{i<j}[a_{ij}\ln(p_{ij})+(1-a_{ij})\ln(1-p_{ij})$
$\displaystyle+\ln(q_{ij}(w_{ij}|a_{ij}))]$ (73)

for conditional models (see Fig. 3b), 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. 2b).

##  V Discussion and conclusions

In a companion paper [18] 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 [18] 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 [50, 28] (also called Bose-Fermi
[47]) distribution; examples of continuous, conditional models are provided by
the $\text{CReM}_{\text{A}}$ and the $\text{CReM}_{\text{B}}$ [48]. 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.

(a) AIC values for the binary log-likelihood (normalized by the yearly
maximum, across models).

(b) AIC values for the full log-likelihood (normalized by the yearly maximum,
across models).

Figure 3: AIC values for the binary and the ‘full’ log-likelihood, normalized
by the yearly maximum, across models, for better visualization, for the
Gleditsch (left panels) and the BACI (right panel) datasets: a lower AIC value
is associated to a better performance. Our plots clearly show that
constraining the degree sequence increases a model performance in reproducing
a network topology: moreover, the integrated exponential model steadily
outperforms the UBCM, signalling that the information gain due to inclusion of
economic variables does not affect the parsimony of its description. For what
concerns the capability of our models in reproducing weighted properties, the
conditional log-normal and gamma models compete, outperforming the other 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_{i}k_{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.

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.

##  VI The DyGyS Python package

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 [18]. The package is available at the
following URL: https://github.com/MarsMDK/DyGyS.

## Acknowledgements

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’. 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 Italian
‘Programma di Attività Integrata’ (PAI) project ‘Prosociality, Cognition and
Peer Effects’ (Pro.Co.P.E.), funded by IMT School for Advanced Studies.

## Appendix A - Conditional models

Any member of the class of conditional models is described by the expression

$Q(\mathbf{W}|\mathbf{A})=\dfrac{e^{-H(\mathbf{W})}}{\int_{\mathbb{W}_{\mathbb{A}}}e^{-H(\mathbf{W})}d\mathbf{W}}$
(74)

the single instances being characterized by different expressions of the
Hamiltonian. Since, however, each Hamiltonian considered here is a sum over
the node pairs, the result

$\displaystyle Q(\mathbf{W}|\mathbf{A})$
$\displaystyle=\dfrac{e^{-\sum_{i<j}H_{ij}(w_{ij})}}{\int_{\mathbb{W}_{\mathbb{A}}}e^{-\sum_{i<j}H_{ij}(w_{ij})}d\mathbf{W}}$
$\displaystyle=\prod_{i<j}\dfrac{e^{-H_{ij}(w_{ij})}}{\left[\int_{m_{ij}}^{+\infty}e^{-H_{ij}(w_{ij})}dw_{ij}\right]^{a_{ij}}}$
$\displaystyle=\prod_{i<j}\dfrac{e^{-H_{ij}(w_{ij})}}{\zeta_{ij}}$ (75)

holds irrespectively from the specific functional form of $H_{ij}(w_{ij})$.
Notice that $m_{ij}$ is the minimum, pair-specific weight allowed by the
model. The identifications

$\displaystyle H_{ij}(w_{ij})\equiv
f(w_{ij}|\underline{\theta},z_{ij}(\underline{\psi}))$ (76)

and

$\mathcal{L}\equiv\ln Q(\mathbf{W}|\mathbf{A})$ (77)

lead to estimate parameters by solving the (coupled) systems of purely
structural equations
$\frac{\partial\mathcal{L}}{\partial\underline{\theta}}=\underline{0}$ and
econometric-like equations
$\frac{\partial\mathcal{L}}{\partial\underline{\tau}}=\underline{0}$, i.e.

$\begin{cases}\sum_{i<j|a_{ij}=1}\left[\frac{\partial
f(w_{ij}|\underline{\theta},z_{ij}(\underline{\psi}))}{\partial\underline{\theta}}+\frac{1}{\zeta_{ij}}\left(\frac{\partial\zeta_{ij}}{\partial\underline{\theta}}\right)\right]=\underline{0},\\\
\sum_{i<j|a_{ij}=1}\left[\frac{\partial
f(w_{ij}|\underline{\theta},z_{ij}(\underline{\psi}))}{\partial
z_{ij}}+\frac{1}{\zeta_{ij}}\left(\frac{\partial\zeta_{ij}}{\partial
z_{ij}}\right)\right]\frac{\partial
z_{ij}}{\partial\underline{\psi}}=\underline{0}.\end{cases}$ (78)

Notice that the estimation of the parameters carried out by maximizing the
conditional likelihood, and letting only the positive weights to be accounted
for, is perfectly consistent with the theory of hurdle models [51, 52]:
although alternative estimation procedures can be devised (see, for example,
[48]), in the present paper, we will stick to the proper, econometric one -
which has been already employed in our companion paper [18], to estimate the
parameters of conditional, discrete-valued models.

###  VI.1 Conditional exponential model

The conditional exponential model is defined by the expression

$H_{ij}(w_{ij})=(\beta_{0}+\beta_{ij})w_{ij}$ (79)

that induces the following node pair-specific partition function

$\zeta_{ij}=\int_{0}^{+\infty}e^{-(\beta_{0}+\beta_{ij})w_{ij}}dw_{ij}=\frac{1}{\beta_{0}+\beta_{ij}}.$
(80)

After the econometric reparametrization, according to which $\beta_{ij}\equiv
z_{ij}^{-1}$, the log-likelihood function of the conditional exponential model
reads

$\mathcal{L}=\sum_{\begin{subarray}{c}i<j\\\
(a_{ij}=1)\end{subarray}}[-(\beta_{0}+z_{ij}^{-1})w_{ij}-\ln(\zeta_{ij})];$
(81)

hence, its maximization leads to the system of equations

$\begin{cases}\sum_{i<j|a_{ij}=1}[\langle w_{ij}|a_{ij}=1\rangle-w_{ij}]&=0\\\
\sum_{i<j|a_{ij}=1}[\langle w_{ij}|a_{ij}=1\rangle-
w_{ij}]\partial_{\underline{\alpha}}(z_{ij}^{-1})&=0\end{cases}$ (82)

where $\langle w_{ij}|a_{ij}=1\rangle=\frac{z_{ij}}{1+\beta_{0}z_{ij}}$.
Notice that we have a condition on the parameters, reading
$\beta_{0}+\beta_{ij}>0$.

###  VI.2 Conditional gamma model

The conditional gamma model is defined by the expression

$H_{ij}(w_{ij})=(\beta_{0}+\beta_{ij})w_{ij}+\xi_{0}\ln(w_{ij})$ (83)

and induces the following node pair-specific partition function

$\displaystyle\zeta_{ij}$
$\displaystyle=\int_{0}^{\infty}e^{-(\beta_{0}+\beta_{ij})w_{ij}}w_{ij}^{-\xi_{0}}dw_{ij}$
$\displaystyle=\frac{\Gamma(1-\xi_{0})}{(\beta_{0}+\beta_{ij})^{1-\xi_{0}}}.$
(84)

After the econometric reparametrization, according to which $\beta_{ij}\equiv
z_{ij}^{-1}$, the log-likelihood function of the conditional gamma model reads

$\displaystyle\mathcal{L}$ $\displaystyle=\sum_{\begin{subarray}{c}i<j\\\
(a_{ij}=1)\end{subarray}}[-(\beta_{0}+z_{ij}^{-1})w_{ij}-\xi_{0}\ln(w_{ij})-\ln(\zeta_{ij})];$
(85)

hence, its maximization leads to the system of equations

$\begin{cases}\sum_{i<j|a_{ij}=1}[\langle w_{ij}|a_{ij}=1\rangle-w_{ij}]&=0\\\
\sum_{i<j|a_{ij}=1}[\langle\ln(w_{ij})|a_{ij}=1\rangle-\ln(w_{ij})]&=0\\\
\sum_{i<j|a_{ij}=1}[\langle w_{ij}|a_{ij}=1\rangle-
w_{ij}]\partial_{\underline{\alpha}}(z_{ij}^{-1})&=0\end{cases}$ (86)

where $\langle
w_{ij}|a_{ij}=1\rangle=\frac{z_{ij}(1-\xi_{0})}{1+\beta_{0}z_{ij}}$ and
$\langle\ln(w_{ij})|a_{ij}=1\rangle=\psi(1-\xi_{0})-\ln(\beta_{0}+z_{ij}^{-1})$.
Notice that we have conditions on the parameters, reading
$\beta_{0}+\beta_{ij}>0$ and $\xi_{0}<1$.

###  VI.3 Conditional Pareto model

The conditional Pareto model is defined by the expression

$H_{ij}(w_{ij})=\xi_{ij}\ln(w_{ij})$ (87)

that induces the following node pair-specific partition function

$\displaystyle\zeta_{ij}$
$\displaystyle=\int_{m_{ij}}^{+\infty}e^{-\xi_{ij}\ln(w_{ij})}dw_{ij}$
$\displaystyle=\int_{m_{ij}}^{+\infty}w_{ij}^{-\xi_{ij}}dw_{ij}$
$\displaystyle=\frac{m_{ij}^{1-\xi_{ij}}}{\xi_{ij}-1}.$ (88)

After the econometric reparametrization, according to which $\xi_{ij}-2\equiv
z_{ij}^{-1}$ and $m_{ij}\equiv w_{min}$, the log-likelihood function of the
conditional Pareto model reads

$\displaystyle\mathcal{L}=\sum_{\begin{subarray}{c}i<j\\\
(a_{ij}=1)\end{subarray}}[-(2+z_{ij}^{-1})\ln(w_{ij})-\ln(\zeta_{ij})];$ (89)

hence, its maximization leads to the system of equations

$\sum_{i<j|a_{ij}=1}[\langle\ln(w_{ij})|a_{ij}=1\rangle-\ln(w_{ij})]\partial_{\underline{\alpha}}(z_{ij}^{-1})=0$
(90)

where
$\langle\ln(w_{ij})|a_{ij}=1\rangle=\ln(w_{min})+\frac{z_{ij}}{1+z_{ij}}$.
Notice that we have conditions on the parameters, reading $w_{min}>0$ and
$\xi_{ij}>2$.

###  VI.4 Conditional log-normal model

The conditional log-normal model is defined by the expression

$\displaystyle H_{ij}(w_{ij})=\xi_{ij}\ln(w_{ij})+\gamma_{0}\ln^{2}(w_{ij})$
(91)

that induces the following node pair-specific partition function

$\displaystyle\zeta_{ij}$
$\displaystyle=\int_{0}^{+\infty}e^{-\xi_{ij}\ln(w_{ij})-\gamma_{0}\ln^{2}(w_{ij})}dw_{ij}$
$\displaystyle=\int_{-\infty}^{+\infty}e^{(1-\xi_{ij})t_{ij}}e^{-\gamma_{0}t_{ij}^{2}}dt_{ij}$
$\displaystyle=\sqrt{\frac{\pi}{\gamma_{0}}}e^{\frac{(\xi_{ij}-1)^{2}}{4\gamma_{0}}}$
(92)

a result that is readily obtained by putting $t_{ij}=\ln(w_{ij})$ and
exploiting the relationship
$\int_{-\infty}^{+\infty}e^{-ax^{2}+bx+c}dx=\sqrt{\frac{\pi}{a}}e^{\frac{b^{2}}{4a}+c}$.

After the econometric reparametrization, according to which
$1-\xi_{ij}\equiv\ln(z_{ij})$, the log-likelihood function of the conditional
log-normal model reads

$\displaystyle\mathcal{L}=\sum_{\begin{subarray}{c}i<j\\\
(a_{ij}=1)\end{subarray}}[(\ln(z_{ij})+1)\ln(w_{ij})-\gamma_{0}\ln^{2}(w_{ij})-\ln(\zeta_{ij})];$
(93)

hence, its maximization leads to the system of equations

$\begin{cases}\sum_{i<j|a_{ij}=1}[\langle\ln^{2}(w_{ij})|a_{ij}=1\rangle-\ln^{2}(w_{ij})]&=0\\\
\sum_{i<j|a_{ij}=1}[\langle\ln(w_{ij})|a_{ij}=1\rangle-\ln(w_{ij})]\partial_{\underline{\alpha}}\ln(z_{ij})&=0\end{cases}$
(94)

where
$\langle\ln^{2}(w_{ij})|a_{ij}=1\rangle=\frac{2\gamma_{0}+\ln^{2}(z_{ij})}{4\gamma_{0}^{2}}$
and $\langle\ln(w_{ij})|a_{ij}=1\rangle=\frac{\ln(z_{ij})}{2\gamma_{0}}$.
Notice that we have a condition on the parameters, reading $\gamma_{0}>0$.

## Appendix B - Integrated models

Any member of the class of integrated models is described by the expression

$Q(\mathbf{W})=\dfrac{e^{-H(\mathbf{W})}}{\int_{\mathbb{W}}e^{-H(\mathbf{W})}d\mathbf{W}}$
(95)

the single instances being characterized by different expressions of the
Hamiltonian. Since, however, each Hamiltonian considered here is a sum over
the node pairs, the result

$\displaystyle Q(\mathbf{W})$
$\displaystyle=\dfrac{e^{-\sum_{i<j}H_{ij}(w_{ij})}}{\int_{\mathbb{W}}e^{-\sum_{i,j}H_{ij}(w_{ij})}d\mathbf{W}}$
$\displaystyle=\prod_{i<j}\dfrac{e^{-H_{ij}(w_{ij})}}{\sum_{a_{ij}=0,1}\int_{\Theta[w_{ij}]=a_{ij}}e^{-H_{ij}(w_{ij})}dw_{ij}}$
$\displaystyle=\prod_{i<j}\dfrac{e^{-H_{ij}(w_{ij})}}{Z_{ij}}$ (96)

holds irrespectively from the specific functional form of $H_{ij}(w_{ij})$.
The identifications

$\displaystyle H_{ij}(w_{ij})\equiv
f(w_{ij}|\underline{\theta},z_{ij}(\underline{\psi}))$ (97)

and

$\mathcal{L}\equiv\ln Q(\mathbf{W})$ (98)

lead to estimate parameters by solving the (coupled) systems of purely
structural equations
$\frac{\partial\mathcal{L}}{\partial\underline{\theta}}=\underline{0}$ and
econometric-like equations
$\frac{\partial\mathcal{L}}{\partial\underline{\psi}}=\underline{0}$, i.e.

$\begin{cases}\sum_{i<j}\left[\frac{\partial
f(w_{ij}|\underline{\theta},z_{ij}(\underline{\psi}))}{\partial\underline{\theta}}+\frac{1}{Z_{ij}}\left(\frac{\partial
Z_{ij}}{\partial\underline{\theta}}\right)\right]=\underline{0},\\\
\sum_{i<j}\left[\frac{\partial
f(w_{ij}|\underline{\theta},z_{ij}(\underline{\psi}))}{\partial
z_{ij}}+\frac{1}{Z_{ij}}\left(\frac{\partial Z_{ij}}{\partial
z_{ij}}\right)\right]\frac{\partial
z_{ij}}{\partial\underline{\psi}}=\underline{0}.\end{cases}$ (99)

The integrated exponential model we have considered in the present paper is
defined by the expression

$H_{ij}(w_{ij})=(\alpha_{i}+\alpha_{j})a_{ij}+(\beta_{0}+\beta_{ij})w_{ij}$
(100)

that induces the following node pair-specific partition function

$\displaystyle Z_{ij}$
$\displaystyle=\sum_{a_{ij}=0}^{1}\int_{\Theta[w_{ij}]=a_{ij}}e^{-(\alpha_{i}+\alpha_{j})a_{ij}-(\beta_{0}+\beta_{ij})w_{ij}}dw_{ij}$
$\displaystyle=1+e^{-(\alpha_{i}+\alpha_{j})}\int_{0}^{+\infty}e^{-(\beta_{0}+\beta_{ij})w_{ij}}dw_{ij}$
$\displaystyle=1+\frac{e^{-(\alpha_{i}+\alpha_{j})}}{\beta_{0}+\beta_{ij}}.$
(101)

After the econometric reparametrization, according to which $\beta_{ij}\equiv
z_{ij}^{-1}$, the log-likelihood function of the exponential model reads

$\mathcal{L}=\sum_{i<j}[-(\alpha_{i}+\alpha_{j})a_{ij}-(\beta_{0}+z_{ij}^{-1})w_{ij}-\ln(Z_{ij})];$
(102)

hence, its maximization leads to the system of equations

$\begin{cases}\langle k_{i}\rangle-k_{i}&=0,\>\forall\>i\\\ \langle
W\rangle-W&=0\\\ \sum_{i<j}[\langle w_{ij}\rangle-
w_{ij}]\partial_{\underline{\alpha}}(z_{ij}^{-1})&=0\\\ \end{cases}$ (103)

where $k_{i}=\sum_{j(\neq i)}a_{ij}$ is the empirical degree of node $i$,
$w_{ij}$ is the empirical, pair-specific weight and $W=\sum_{i<j}w_{ij}$ is
the empirical, total weight; $\langle k_{i}\rangle=\sum_{j(\neq i)}p_{ij}$,
$\langle w_{ij}\rangle$ and $\langle W\rangle=\sum_{i<j}\langle w_{ij}\rangle$
are their expected counterparts.

Notice that we have a condition on the parameters, reading
$\beta_{0}+\beta_{ij}>0$.

## Appendix C - Turning structural models into econometric models

So far, we have derived two classes of models, by explicitly solving the
constrained maximization of a number of functionals derived from the KL
divergence. As the functional form of the probability distributions belonging
to the two classes (solely) depends on the enforced constraints, such models
‘are born’ as purely structural ones.

In order to turn them into candidate models to be employed for econometric
purposes, we need to properly transform (some of) the Lagrange multipliers
into functions of the econometric quantities of relevance for the problem at
hand. In this respect, the theory of GLMs provides helpful suggestions about
how to proceed; besides, one can figure out some (sets of) basic requirements
such a transformation should satisfy:

  * •

the transformation should turn the expected values $\langle w_{ij}\rangle$ and
$\langle w_{ij}|a_{ij}=1\rangle$ into positive, monotonically increasing
functions of $z_{ij}$;

  * •

the transformation should not violate the mathematical requirements to have
well-defined (first and second) distribution moments.

In what follows, we will focus on the conditional models.

_Conditional exponential model._ It is characterized by the expression

$\langle w_{ij}|a_{ij}=1\rangle=\frac{1}{\beta_{0}+\beta_{ij}}$ (104)

that can be turned into an econometric one by posing

$\langle w_{ij}|a_{ij}=1\rangle=\frac{1}{\beta_{0}+\beta_{ij}}\equiv
g(z_{ij})$ (105)

according to the prescription informing the so-called generalized linear
models (GLMs). Before specifying the functional form of $g$, let us consider
that the dyadic parameter $\beta_{ij}$ must be decreasing in $z_{ij}$ \- a
requirement that can be justified upon identifying $\beta_{ij}$ as the ‘shadow
price’ that countries $i$ and $j$ have to pay to trade a unity of goods [53];
analogously, $\beta_{0}$ can be interpreted as modelling a global tax that
everyone has to pay to exchange goods - independently of its trade ‘capacity’.
These considerations lead us to impose $\beta_{ij}\equiv z_{ij}^{-1}$, a
choice inducing the expression

$\langle
w_{ij}|a_{ij}=1\rangle=\frac{1}{\beta_{0}+z_{ij}^{-1}}=\frac{z_{ij}}{1+\beta_{0}z_{ij}}$
(106)

which violates none of the requirements listed at the beginning of the
section.

_Conditional gamma model._ It is characterized by the expressions

$\displaystyle\langle w_{ij}|a_{ij}=1\rangle$
$\displaystyle=\frac{1-\xi_{0}}{\beta_{0}+\beta_{ij}},$ (107)
$\displaystyle\langle\ln(w_{ij})|a_{ij}=1\rangle$
$\displaystyle=\psi(1-\xi_{0})-\ln(\beta_{0}+\beta_{ij})$ (108)

(where $\psi(x)=\Gamma^{\prime}(x)/\Gamma(x)$ is the digamma function) that
can be turned into econometric ones by posing $\beta_{ij}\equiv z_{ij}^{-1}$,
according to considerations which are analogous to those driving the
econometric reparametrization of the conditional, exponential model. This
choice induces the expressions

$\displaystyle\langle w_{ij}|a_{ij}=1\rangle$
$\displaystyle=\frac{1-\xi_{0}}{\beta_{0}+z_{ij}^{-1}}=\frac{(1-\xi_{0})z_{ij}}{1+\beta_{0}z_{ij}},$
(109) $\displaystyle\langle\ln(w_{ij})|a_{ij}=1\rangle$
$\displaystyle=\psi(1-\xi_{0})-\ln(\beta_{0}+z_{ij}^{-1});$ (110)

notice that the conditional, exponential model is recovered in case
$\xi_{0}=0$ (i.e. when the constraint on the sum of the logarithms of weights
is switched-off).

_Conditional Pareto model._ It is characterized by the expression

$\langle
w_{ij}|a_{ij}=1\rangle=\left(\dfrac{\xi_{ij}-1}{\xi_{ij}-2}\right)m_{ij}$
(111)

that can be turned into an econometric one by posing

$\langle
w_{ij}|a_{ij}=1\rangle=\left(\dfrac{\xi_{ij}-1}{\xi_{ij}-2}\right)m_{ij}\equiv
g(z_{ij})$ (112)

according to the prescription informing the GLMs. Upon considering that 1) the
(conditional) expected value is well defined only if $\xi_{ij}-2>0$ and that
2) a linear relationship between the former and $z_{ij}$ would be desirable, a
suitable reparametrization may read $\xi_{ij}-2\equiv z_{ij}^{-1}$ and
$m_{ij}\equiv w_{min}$, in turn leading to

$\langle w_{ij}|a_{ij}=1\rangle=(1+z_{ij})w_{min}$ (113)

which violates none of the requirements listed at the beginning of the
section.

_Conditional log-normal model._ It is characterized by the expressions

$\displaystyle\langle\ln(w_{ij})|a_{ij}=1\rangle$
$\displaystyle=\frac{1-\xi_{ij}}{2\gamma_{0}},$ (114)
$\displaystyle\langle\ln^{2}(w_{ij})|a_{ij}=1\rangle$
$\displaystyle=\frac{2\gamma_{0}+(1-\xi_{ij})^{2}}{4\gamma_{0}^{2}}.$ (115)

Upon considering that the logarithm of weights can admit negative values, i.e.
when $w_{ij}\in(0,1)$, and that there are no theoretical restrictions on the
sign of $\langle\ln(w_{ij})|a_{ij}=1\rangle$, a suitable reparametrization may
read $1-\xi_{ij}\equiv\ln(z_{ij})$, in turn leading to

$\displaystyle\langle\ln(w_{ij})|a_{ij}=1\rangle$
$\displaystyle=\frac{\ln(z_{ij})}{2\gamma_{0}},$ (116)
$\displaystyle\langle\ln^{2}(w_{ij})|a_{ij}=1\rangle$
$\displaystyle=\frac{2\gamma_{0}+\ln^{2}(z_{ij})}{4\gamma_{0}^{2}}$ (117)

which violate none of the requirements listed at the beginning of the section.

## References

  * [1] M. Barigozzi, G. Fagiolo and D. Garlaschelli, Multinetwork of international trade: A commodity-specific analysis, Phys. Rev. E 81, 046104 (2010)
  * [2] A. Fronczak, Structural Hamiltonian of the international trade network, Acta Phys. Pol. B 5, 31–46 (2012)
  * [3] A. Fronczak and P. Fronczak, Statistical mechanics of the international trade network, Phys. Rev. E 85, 056113 (2012).
  * [4] M. Á. Serrano and M. Boguña, Topology of the world trade web, Phys. Rev. E 68, 015101(R) (2003).
  * [5] D. Garlaschelli and M. I. Loffredo, Fitness-Dependent Topological Properties of the World Trade Web, Phys. Rev. Lett. 93, 188701 (2004).
  * [6] D. Garlaschelli and M. I. Loffredo, Structure and evolution of the world trade network, Phys. A: Stat. Mech. Appl. 355, 138–144 (2005).
  * [7] G. Fagiolo, J. Reyes and S. Schiavo, The evolution of the world trade web: a weighted-network analysis., J. Evol. Econ. 20:479–514 (2010).
  * [8] G. Fagiolo, J. Reyes and S. Schiavo, On the topological properties of the world trade web: A weighted network analysis, Phys. A: Stat. Mech. Appl. 387 (15), 3868-3873 (2009).
  * [9] F. Schweitzer, G. Fagiolo, D. Sornette and F. Vega-Redondo, Economic networks: the new challenges, Science 325(5939):422-5 (2009).
  * [10] S. Vitali, J. B. Glattfelder and S. Battiston, The Network of Global Corporate Control, PLoS ONE , 6(10), e25995 (2011).
  * [11] S. Schiavo, J. Reyes, G. Fagiolo, International trade and financial integration: a weighted network analysis, Quant. Finance 10(4), 389-399 (2010).
  * [12] G. Fagiolo, J. Reyes and S. Schiavo, World-trade web: Topological properties, dynamics, and evolution, Phys. Rev. E 79 036115 (2009).
  * [13] J. Tinbergen, Shaping the World Economy: Suggestions for an International Economic Policy, New York: The Twentieth Century Fund (1962).
  * [14] A. G. Wilson, The Use of Entropy Maximising Models, in the Theory of Trip Distribution, Mode Split and Route Split, J Trans Econ Pol, Vol. 3(1): pp. 108-126, (1969),
  * [15] H Wang, X. Yan and J. Wu, Free utility model for explaining the social gravity law, J. Stat. Mech. 033418 (2021).
  * [16] G. Fagiolo, The international-trade network: gravity equations and topological properties, J. Econ. Interact. Coord. 5:1–25 (2010).
  * [17] E. Helpman, M. Melitz, and Y. Rubinstein, Estimating Trade Flows: Trading Partners and Trading Volumes, Q. J. Econ. 123[2], 441-487 (2008).
  * [18] M. Di Vece, D. Garlaschelli and T. Squartini, Gravity models of networks: integrating maximum-entropy and econometric approaches, Phys. Rev. Research 4, 033105 (2022).
  * [19] M. Burger, F. van Oort and G. Linders, On the Specification of the Gravity Model of Trade: Zeros, Excess Zeros and Zero-Inflated Estimation, Spat. Econ. Anal. 4:2, 167-190, (2009).
  * [20] J. S. Long, Regression models for categorical and limited dependent variables, Thousand Oaks, CA: Sage, (1997).
  * [21] T. Squartini and D. Garlaschelli, Maximum-Entropy Networks: Pattern Detection, Network Reconstruction and Graph Combinatorics, Springer International Publishing (2017).
  * [22] G. Cimini, T. Squartini, F. Saracco, D. Garlaschelli, A. Gabrielli and G. Caldarelli, The statistical physics of real-world networks, Nat. Rev. Phys. 1, 58–71 (2019).
  * [23] C. E. Shannon, A mathematical theory of communication,The Bell system technical journal 27(3), 379-423 (1948).
  * [24] T.M. Cover and J.A. Thomas, Elements of information theory (John Wiley & Sons, 1999).
  * [25] S. Wasserman and K. Faust, Social Network Analysis: Methods and Applications (Cambridge University Press, 1994).
  * [26] T. Squartini, G. Fagiolo and D. Garlaschelli, Randomizing world trade. I. A binary network analysis, Phys. Rev. E 84, 046117 (2011).
  * [27] T. Squartini, G. Fagiolo and D. Garlaschelli, Randomizing world trade. II. A weighted network analysis, Phys. Rev. E 84, 046118 (2011).
  * [28] R. Mastrandrea, T. Squartini, G. Fagiolo and D. Garlaschelli, Reconstructing the world trade multiplex: the role of intensive and extensive biases, Physical Review E 90(6), 062804 (2014).
  * [29] T. Squartini and D. Garlaschelli, Jan Tinbergen’s legacy for economic networks: from the gravity model to quantum statistics, in Econophysics of Agent-Based Models, pp. 161-186 (Springer, Cham, 2014).
  * [30] A. Almog, T. Squartini and D. Garlaschelli, A GDP-driven model for the binary and weighted structure of the International Trade Network, New Journal of Physics 17(1), 013009 (2015).
  * [31] A. Almog, T. Squartini and D. Garlaschelli, The double role of GDP in shaping the structure of the international trade network, International Journal of Computational Economics and Econometrics 7(4), 381-398 (2017).
  * [32] A. Almog, R. Bird and D. Garlaschelli, Enhanced gravity model of trade: reconciling macroeconomic and network models, Frontiers in Physics 7, 55 (2019).
  * [33] O. Sagarra, C. J. Pérez Vicente, and A. Díaz-Guilera, Role of adjacency-matrix degeneracy in maximum-entropy-weighted network models, Phys. Rev. E 92, 052816 (2015).
  * [34] L. Bargigli, Statistical ensembles for economic networks, J. Stat. Phys. 155, 810–825 (2014)
  * [35] P R Herman, Modeling Complex Network Patterns in International Trade, Rev World Econ 158, 127–179 (2022).
  * [36] J. A. Nelder and R. W. M. Wedderburn, Generalized Linear Models, J. R. Stat. Soc. Ser. A (General) Vol. 135, No. 3 (1972), pp. 370-384 (15 pages)
  * [37] E. T. Jaynes, Information Theory and Statistical Mechanics, Phys. Rev. 106, 620 (1957)
  * [38] E. T. Jaynes, Information Theory and Statistical Mechanics II, Phys. Rev. 108, 171 (1957)
  * [39] E. T. Jaynes, On the rationale of maximum-entropy methods, Proceedings of the IEEE 70:9, pp. 939-952 (1982)
  * [40] S. Kullback, R. A. Leibler, On Information and Sufficiency, Ann. Math. Statist. 22(1): 79-86 (1951). DOI: 10.1214/aoms/1177729694
  * [41] I. J. Goodfellow, J. Pouget-Abadie, M. Mirza, B. Xu, D. Warde-Farley, S. Ozairy, A. Courville, Y. Bengioz, Generative adversarial nets, Adv. Neural Inf. Process. Syst. 3(11) (2014).
  * [42] S. H. Walker and D. B. Duncan, Estimation of the probability of an event as a function of several independent variables, Biometrika Vol. 54, No. 1/2 (June 1967), pp. 167-179 (13 pages).
  * [43] G. Caldarelli, A. Capocci, P. De Los Rios, and M. A. Muñoz, Scale-Free Networks from Varying Vertex Intrinsic Fitness, Phys. Rev. Lett. 89:25, 258702
  * [44] M. Visser, Zipf’s law, power laws and maximum entropy, New J. Phys. 15 (2013).
  * [45] K. Gleditsch, Expanded Trade and GDP Data, J. Confl. Resolut., 46(5), 712-724 (2002)
  * [46] K. Head and T. Mayer, Gravity Equations: Toolkit, Cookbook, Workhorse., Handbook of International Economics, Vol.4, Elsevier (2014)
  * [47] D. Garlaschelli and M.I. Loffredo, Generalized Bose-Fermi statistics and structural correlations in weighted networks, Physical Review Letters 102(3), 038701 (2009).
  * [48] F. Parisi, T. Squartini and D. Garlaschelli, A faster horse on a safer trail: generalized inference for the efficient reconstruction of weighted networks, New J. Phys.22 053053 (2020).
  * [49] H. Akaike, Information theory and an extension of the maximum likelihood principle, in: Petrov, B.N. and Csaki, F., Eds., International Symposium on Information Theory, 267-281, (1973).
  * [50] R. Mastrandrea, T. Squartini, G. Fagiolo and D. Garlaschelli, Enhanced reconstruction of weighted networks from strengths and degrees, New J. Phys. 16, 043022 (2014).
  * [51] J. G. Cragg, Some Statistical Models for Limited Dependent Variables with Application to the Demand for Durable Goods, Econometrica 39(5): 829–844 (1971).
  * [52] J. Mullahy, Specification and testing of some modified count data models, J Econom 33(3): 341–365 (1986).
  * [53] L. Bargigli, S. Viaggiu and A. Lionetto, A Statistical Test of Walrasian Equilibrium by Means of Complex Networks Theory, J.Stat. Phys. 165, 351–370 (2016).