A stochastic generative model of the World Trade Network

The World Trade Network (WTN) is a network of exchange flows among countries whose topological and statistical properties are a valuable source of information. Degree and strength (weighted degree) are key magnitudes to understand its structure and generative mechanisms. In this work, we describe a stochastic model that yields synthetic networks that closely mimic the properties of annual empirical data. The model combines two popular mechanisms of network generation: preferential attachment and multiplicative process. Agreement between empirical and synthetic networks is checked using the available series from 1962 to 2017.

International trade is the circulatory system of global economy. The flow of goods and services may be modelled as a bipartite network with two classes on nodes: exporters and importers. These nodes are tied by the existence and intensity of exchanges.
The study of the statistical and topological properties of the global trade as a network has been an active research topic in recent years. In these studies, degree and strength are two key connectivity properties of each node. Degree is the number of links to nodes of the opposite class, while strength is the sum of weights of those links. For WTN nodes, degree is the number of countries they trade with. Each link is weighted by the trade volume between both nodes. Strength accounts for the total ammount of trade of a particular node and equates the sum of weights ot its links. Probability distributions of degree and strength shed light on the network properties 1 . If the network is scale-free, these distributions follow a power-law 2 , also known as Pareto's law in Economics. Trying to fit empirical data to a power-law has been a common procedure in network analysis, despite the fact that this behavior is not universal 3 .
The log-normal distribution is a second possible choice in many fields, but specially in Economics 4 . The discussion about the best model has run for a long time, because both may fit the long tail of some empirical series 5 . The hypothesis of log-normality dates back to the end of the XIX th Century, but it was the Dutch astronomer Cornelius Kapteyn who first described its pervasive presence and proposed a generative model 6 . He used the growth of plants and animals to explain the multiplicative process that lies behind the emergence of log-normality. Three decades later, Robert Gibrat wrote about the law of proportional effect, that stated that the average and variance of growth rates for many industrial processes are nearly constant 7 and by virtue of the Central Limit Theorem the distribution tends towards log-normality. Fitting empirical data to a log-normal distribution to test Gibrat's hypothesis was a recurrent research topic, in special working with the size of firms 8 , despite the fact that "the upper tail resembles the Pareto's distribution" in the own words of the classical work of Ijiri and Simon 9 . Stanley et al. found that the distribution of logarithmic annual growth rates was better described by an exponential 10 .
The discussion about the validity of Gibrat's law to describe empirical distributions appears in different fields of Economics, with examples supporting 11,12 and rejecting it 13,14 on statistical grounds. Depending on the particular series, the study may be inconclusive, as in the particular case of Italian hospitality industry 15 .
The dilemma between power-law and log-normal distribution to describe empirical trade volume series 13 has also been a recurrent topic of the World Trade Network analysis. WTN degree and strength distributions do not show a power-law pattern (see Supplementary Information, Fig. S1) but over the skewed right tails. So, Garlaschelli and Loffredo showed 16 that WTN cumulative degree distribution is not a power-law as claimed in a precedent work 17 , but that WTN is scale-free using the cumulative GDP distribution as a hidden variable 18  www.nature.com/scientificreports www.nature.com/scientificreports/ Benedictis et al., hold the view that the in-degree and out-degree right tail distributions follow a power-law over a wide range of values 22 . On the other hand, log-normality seems a good choice to describe degree and strength distributions for Fagiolo et al. 23 and Barigozzi et al. 24 .
Empirical degree or strength data follow neither pure log-normal distributions, nor sharp power-law over all the range of values. A possible explanation is that the statistical properties of the WTN matrixes could be driven by two processes of different nature acting at different time scales, the first driving a fast network build-up, and the second one over a much longer period.
We have designed our model inspired by two precedents, Ijiri and Simon principles to explain log-normality in the size distribution of firms 9 and the generative approach to a well-known class of bipartite networks in a quite different field, Ecology.
Ijiri and Simon proposed a stochastic model to emulate the size growth of firms using two basic principles: the probability that a new business chance is captured by an existing company is proportional to its size and the probability of arrival of a new company is constant over time. This model was well-suited for a unipartite network that acts on a long-term scale.
In Ecology, mutualistic networks work on the principle of mutual benefit for both classes of nodes, like plants and pollinators. In a different way of human trade they also exchange goods (pollen, nectar) and services (fecundation, protection). For this kind of networks, there has been a long discussion on the best generative model to explain why their degree and strength distributions follow a power-law or a truncated power-law (see a review by Guimaraes et al. 25 ). Field data sampling in Ecology is a hard task, and empirical data may consist of just one campaign. Generative models in this walk of research try to reproduce the main statistical parameters of the recorded adjacency matrix, from a pure phenomenological point of view. In this work, we adopt the same approach to work with the historical world trade series, with the advantage that the series spans 56 years.
For our study the WTN is a set with two different guilds of nodes: Exporters and Importers. There exists a directed link from country A as exporter to country B as importer (Exporter A → Importer B ), when Exporter A sells merchandises to Importer B . The weight of the link is the volume they trade during that year. Any country may be both exporter and importer, so the opposite flow from B to A happens between a different pair of nodes (Exporter B → Importer A ) and with a different weight. Exporter nodes may link to any Importer, but never among them (and viceversa), and the weight of the relationship is the trade volume flowing from the Exporter to the Importer. Under these circumstances, the WTN is a weighted bipartite network 26 . The adjacency matrix has as many rows as exporters and as many columns as importers. There is one empirical matrix per year that stores the flow from Exporter i to Importer j ∀ i, j. The number of nodes and links changes from year to year, due to geopolitical events and the non-stationary nature of trade itself.
Our approach takes into account two observations: dominant countries, both as exporters or importers, attract with high probability new trade chances and global yearly degree and volume trade distributions are approximately log-normal. In this paper, we describe a generative stochastic synthetic model of trade (Synthrade) that mixes two processes, preferential attachment and proportional effect, to mimic the properties of the weighted matrixes of WTN historical data series. The first one works during a very short period of simulation time and yields a scale-free heterogeneous distribution of trade volume for exporter and importer guilds. The second, acts on a much larger time scale and by its multiplicative effect makes the yearly trade distributions log-normal.
The term growth is pervasive in Economics and Network Science, but not always with the same meaning. In this paper it makes reference to both possible acceptations, it is mandatory to explain how we use it. When Gibrat writes about firms' growth, it is clear that he refers to the yearly increase of some proxy variable of its size, like income or number of employees. If a researcher analyzes the world trade volume evolution over the last century, growth means a rate of change of that magnitude. On the other hand, a network is usually a dynamic entity that experiments growth, by the addition of new nodes, links or the increase of the weight of links. When dealing with the synthetic World Trade Network model, network growth is not the same than growth rate of trade volume. In fact, the model does not deal with the second meaning, but just with network growth.

Results
Model. The model of the international trade network is created from scratch and computed year by year.
Nodes arrive during a very short period of simulation time, at a decreasing pace, and they attach to nodes of the opposite guild with a probability that is proportional to the volume traded by those countries.
An important difference with Ijiri and Simon approach is that we do not try to model trade volume evolution over a long span of time but just the statistical properties of each annual empirical data, using as parameters just three data of the original set as we explain in this section.
A weighted bipartite network t W ( ) , represents the trade, where the value of each cell i, j is the aggregated volume between Exporter i and Importer j at simulation instant t. The strength of Exporter i is defined as The total trade volume of the synthetic network at instat t is the There is another matrix t P ( ) , that accounts for the probability of trade among countries. At each simulation step t, t P ( ) ( 1) The model simulates the build-up of the network as two simultaneous processes with different time scales: the aggregation of new nodes (we call it node aggregation regime) and the increase of trade volumes (weight aggregation regime). At t = 0, the network is just a tiny seed of connected nodes.
The first process drives the aggregation of nodes of both guilds under the hypothesis of neutrality. That means that the birth probability depends only on the number of nodes of the same guild. Each new node joins the network with a probability that is inverse proportional to the number of nodes in its guild, N γ . A new node may www.nature.com/scientificreports www.nature.com/scientificreports/ appear with probability P birth γ(t), where γ stands for Exporters or Importers. This value must decrease as the number of nodes raises, and for the sake of simplicity we use the following form: where λ γ is an intrinsic birth rate parameter of guild γ and N γ (t) is the number of nodes at simulation step t (Fig. 1a). Under these assumptions, the time interval between the birth of two consecutive nodes is ∆ = λ γ γ t N t ( ) . If the network starts with a seed of size σ × σ nodes, the formation time is the sum of an arithmetic series: . We set the condition that birth probabilities equal to 1 when the formation process starts. At that time there exists a connected seed of σ × σ nodes (3 of each class in our practical implementation).
So, the parameter τ γ depends on the final number of nodes of the network. It changes from year to year for the historical data series, and we perform each simulation experiment using the empirical number of nodes.
The number of nodes is not the same for both guilds due to the applied filtering (see Data Preprocessing section), although figures are quite similar. The formation times of each guild would be different.
In our model we choose for both guilds the maximum value of τ E and τ I . This choice has very little impact in the numerical simulation but sets a sharp boundary between both regimes.
During the first, that we call node aggregation regime, the number of nodes grows up. Each new node joins the network with a preferential mechanism. Let t W ( ) be the weighted matrix at simulation instant t. When an importer [exporter] node appears the probability of attachment to any of the exporter [importer] nodes is , as defined above. If the Bernouilli trial is successful, a link with weight 1 will connect both nodes (Fig. 1b). We call this quantum of trade a token, and is the minimum possible value of any non-empty cell of the network. The newborn node may tie to more than one existing node of the opposite class, although it is unlikely according to the procedure. Importer number 5 appears because of the birth process driven by parameter λ I . It may attach to the existing exporter nodes with a probability proportional to the global trade that each exporter handles at that instant. For this example, the set of Bernouilli experiments yields a new link Exporter 2 → Importer 5 (green square). Its weight is the minimum possible trade volume that we call token. (c) The weight aggregation regime happens from the initial step to the end of the simulation, we have represented it just after the aggregation of Importer 5 for simplicity. A new token attaches to the matrix attracted by a probability that is proportional to the product of the marginal probabilities of Exporter i and Importer j . Probability is nonzero even for empty cells, this the mechanism that allows the birth of new links. Anyway, most tokens are attracted towards main trade associations, such as the pair Exporter 1 and Importer 1 (red square). (2019) 9:18539 | https://doi.org/10.1038/s41598-019-54979-1 www.nature.com/scientificreports www.nature.com/scientificreports/ The weight aggregation regime starts at t = 0 and explains network weight growth. Each token may be seen as a trade chance. At each simulation step a set of new trade chances adds to the t W ( ) EI i j , matrix, running Bernouilli trials for each cell with probabilities t P ( ) ( 1) In case of success one token adds to that cell. By definition, the sum of the elements of t P ( ) is one, so, in average the global trade raises one token per simulation step. Note that an empty cell has a small probability of capturing a trade chance, and by this mechanism the number of links keep on growing even after the network formation instant (from now on, t F , see Fig. 1c). Most trade chances go to existing links, and the time interval between the appearance of new links becomes longer. Eventually, the synthetic matrix reaches the number of links of the empirical matrix and simulation stops at that moment, that we call t T . This aggregation process lasts much longer than the node aggregation regime  t t ( ) T F . Finally, as a result of both processes, we obtain the synthetic trade matrix, whose statistical properties may be compared to those of the empirical one (see Fig. 2 and Supplementary Fig. S2c,d). Note that for each year, only three parameters of the empirical network set the final configuration of the synthetic matrix: number of nodes of each guild and number of links of the network.

network formation dynamics. As
is normalized, the average number of successes of the N Exp (t)N Imp (t) Bernouilli trials is just 1, and the number of tokens equates the number of simulation steps. This result is straightforward after the formation instant, when the size of the matrix is constant (see Supplementary  Information, Fig. S2).
During the node aggregation regime the relationship is more complex, as newborn nodes add new links. This figure is of the order of magnitude of N Exp (t F ), less than a 10% of the number of tokens at the formation instant , while the token falls in a non-empty cell with probability P full (t) = 1 − P empty (t).
Probabilities are not constant over the simulation spain. Figure 3a shows some interesting properties to simplify the analysis. The maximum value of t P ( ) EI i j , corresponds to the strongest trade flow of the matrix. This Exporter/Importer pair attracts the majority of tokens and the probability remains almost constant over time. As trade volume tends to attach to existing links, the sum of probabilities of the set of empty cells decreases, an so the minimum value of the matrix. Two interesting values are those of the probability of the cell where the last link of the simulation appeared before any given simulation instant and the median of the probability matrix. The plot shows a sharp contrast before and after t F . During the node aggregation regime, each attached node comes with new empty cells that will have low probability values in contrast with those of populated cells. After t T , matrix size is fixed, no new empty cells may appear, they only may evolve to populated cells. The probability of the last filled cell remains almost constant for the last two thirds of the simulation time, as well. The probability matrix t P ( ) EI i j , plays in Synthrade the same role of the growth distribution in Ijiri and Simon's model, as ict acts as the network growth driver. It is strongly right-skewed, but the median remains almost constant, and this property yields log-normality in synthetic trade degree and strentgth. As a consequence of this evolution, time between arrival of new links does not grow linearly. Assume that the global strength value is m, and the number of links l, when the first link appears after t F . As the variation of P empty (t) is so slow, the addition of new links may be viewed as a Bernouilli process. The mean time between the appearance of link l and link l + 1 is 1/P empty (t). During that interval, the non-empty cells attract a quantity m of tokens, plus the token that falls in the empty cell, so t l+1 − t l = m + 1. To add a second link, the network adds m + 1 tokens to the non-empty cells, plus the token of the new link, so t l+2 − t l+1 = m + 2. The time to add L links is given by the formula: So, for a linear growth of L links, the strength obeys a quadratic law, and the simulation time is O(n 2 ) (see Supplementary Information, Fig. S2a). During the node aggregation regime, the matrix size is not constant and the hypothesis of almost stationary probabilities is not longer valid. Anyway, the number of links and nodes of both guilds follow quadratic patterns as well (see Supplementary Fig. S2b). www.nature.com/scientificreports www.nature.com/scientificreports/ Strength and degree distributions. The plots of empirical networks suggested that strength and degree follow a power law s ∝ k β (see Supplementary Information, Fig. S4). To estimate the slope β, we use the cumulative distribution 27 and fit a log-log regression line to get β + 1.
As Barrat et al. showed 28 , in the lack of correlations between the weight of the links and the degree of the nodes, the strength of a node is simply proportional to its degree, yielding an exponent 1, consequently, strength and degree provide the same information on the system. In our study, the values of β are always above 1 for empirical and synthetic networks ( Fig. 4 and Supplementary Information, Fig. S3). The strength of the most connected nodes tends to have a value higher than the corresponding to a random process of weights. In other terms, main exporters (importers) attract more volume than if trade flow were just proportional to the number of countries they trade with.
The distribution of the empirical weight matrices tends towards a log-normal, and the synthetic ones closely follow their shapes. To test this hypothesis we perform a goodness of fit analysis of the logarithmic values of the empirical strength distributions of both Exporters and Importers, with the Lilliefors test. Setting a p-value of 0.05, the hypothesis of non-normality is rejected for 8 out of 56 importer distributions and none of exporter distributions (Fig. 5a,b, see Supplementary Information, Table S1).
To measure the shape similarity between the normalized empirical strength distributions and their synthetic counterparts we set the following procedure. For each year, we build 30 synthetic networks, running 30 independent experiments with the same three model parameters values (number of importers and exporters and number  www.nature.com/scientificreports www.nature.com/scientificreports/ of links) of the empirical network. Then, we perform a Kolmogorov-Smirnov distance test between the empirical matrix and each one of the synthetic matrixes. According to these results the model mimics quite well the exporter distributions. Boxplots in Fig. 5c,d are the set of K-S test values for the marginal exporter and importers distributions of each year. The horizontal green line, is set at p = 0.05. For this test, the null hypothesis is that the empirical and the synthetic distributions are different. The hypothesis may be rejected for almost all exporter distributions. On the other hand, the null hypothesis of dissimilarity can't be rejected for importers during the 1990's, when connectance dropped in the aftermath of the dissolution of the Soviet Union. Figure 6 shows, as a typical example, a visual comparison of degree and strength distributions for year 1986.
Marginal probability distributions P E (t) and P I (t) are strongly right-skewed at t T (Fig. 3b,c). Their shape tends towards log-normality as the simulation enters the weight aggregation regime and network size remains constant. The explanation is that empty cells have more choices over time of being populated, and from that instant on, they join the multiplicative process. Despite the fact that their strength is small, the distribution is more even because of the increase of the fraction of non-empty cells. The same pattern is observed for the probability matrix P EI .
Log-normality of the strength and degree distributions is a result of the multiplicative process 29 .
for the sake of simplicity. One simulation step later After m simulation steps the strength will be: Taking natural logarithms, the expression becomes: are independent, as each token addition experiment only depends on the values of P EI at simulation step t. As the average value of is also quite smaller than 1. The standard deviation of a sum of independent random variables is the sum of the standard deviations, and so, the value of Var log W m ( ( ( ))) E j is greater than 1 and finite. It holds also that +  Var log X k Var log W m ( (1 ( ))) ( ( ( ))) E j E j . The Lindeberg-Levy conditions are met and the distribution of W E approaches a log-normal. www.nature.com/scientificreports www.nature.com/scientificreports/

Discussion
In this work we have proposed a generative model (that we call Synthrade) to build synthetic bipartite networks that closely resemble the statistical properties of size distributions of the historical series. The bipartite hypothesis is a valid alternative to the time-evolving directed network model 16 and more expressive than unidirected representations 17,26 .
Although the World Trade Network is a very complex system that emerges from the aggregation of multiple choices at different levels (individuals, companies, regulatory bodies, and so on) Synthrade reproduces the statistical properties of the recorded annual trade data working on two principles: preferential attachment and network strength growth by product of probabilities. Trade flows are driven by the relative strength of exporter and importer nodes, any new trade chance attaches to a cell of the bipartite trade matrix with a probability that is the product of the marginal probabilities of both edges.
Network strength grows by probability, though, this fact does not explain initial network formation. The preferential attachment mechanism that acts before t F is the basis of the observed strength heterogeneity in empirical data. Synthrade is naive in the sense that it loses memory from year to year and trade agreements and contracts do not work like that. The explanation on how this relates to real world trade could be the following. Node aggregation is a very fast process compared to the total simulation time. Nodes that arrive early enjoy higher probability of capturing new trade chances. In addition, aggregation speed decreases over simulation time, so first comer nodes will be the main trade attractors, while minor players will enter the network later. The parameter λ governs that speed. Simulation time should not be regarded as a proxy of real time. Simulation steps only represent the arrival of a new trade token. With this idea in mind, the node aggregation regime is not some small period of each January, but an almost instantaneous and convenient mechanism to generate heterogeneity.
Furthermore, the relationship degree-strength of the synthetic networks follows a power-law, as it happens in the empirical series. Any incoming node connects to major sellers (or buyers) as a first choice and links appear as time goes by during the network build-up period, but at a diminishing pace. Once the network has reached its maximum size, incoming trade chances are the only mechanism to create new links. As noted by Garlaschelli et al. 16 , the WTN may be explained by the network fitness model 18 using GDP as the hidden variable. In Synthrade, the force that fosters the birth of new links through mutual benefit is proportional to the product of the strengths of both nodes. As new trade opportunities tend to join the main flows Exporter-Importer, dominant nodes attract more and more volume. The exponent of the degree-strength distribution shows that they are over-trading, while smaller countries strive to play their role in this global network 28 .
As trade chances (tokens) have a small and constant weight, their power to modify the underlying probability matrix weakens over time. Weight and probability matrixes tend towards quasi-stationary states and the long-term effect of this behaviour is a multiplicative growth process that pushes the strength and degree distributions towards a log-normal.
Synthrade is a remarkable simple model that does not take into account many constraints, as distance, trade agreements, tariffs or competitive advantages 30 . In particular, it ignores factors such as distance or some GDP proxy, unlike the family of gravity bilateral trade models 21,31 . It works on a global scale trying to reproduce the network degree sequences that are highly informative instead 1 . On the other hand, it does not require to know in advance that a trade link exist between a pair of countries. A recent advance by Almog, Squartini and Garlaschelli has reconciled the accuracy of gravity models to estimate the intensity of trade flows with its inability to discover network topology, using GDP as a fitness parameter 32 . Synthrade does not include macroeconomic information, just the interaction matrixes data because they are enough to reach its goal of reproducing trade volume topological properties.
This set of simplifications opens the way to further improvements. The node aggregation regime could benefit of the addition of rules of affinity or prohibitions 33 . Similarly, the probability matrix might include some modulation function to model geographical distance 34,35 . Finally, it would be interesting to check if degree and strength follow similar patterns when trade is restricted to some particular good or category of goods, and if so, the bipartite WTN shows fractal properties 36 .

Methods
Data preprocessing. We work with the historical data series from 1962 to 2017, available at The Observatory of Economic Complexity 37 . For each year there is one file that holds per line the ISO code of the exporter country, the ISO code of the importer, a category key of the traded merchandise and the yearly amount in dollars. For this paper we have worked with global trade, just adding the monetary value of all exchanged goods between each pair Exporter/Importer. With this information we build the empirical weight matrix (Fig. 2a).
Quality issues arise from two different product classification systems (SITC and HC), two different statistical sources (from 1962 to 2000 38 and from 2001 to 2017 39 ), two different harmonization approaches 38,40 and geopolitical events 41,42 during such a long span. In addition, trade networks are strongly meshed (average connectance is 0.49) because of the existence of very weak links. To overcome these problems we apply two mild mitigation strategies.
There is a step decrease of connectance (ratio of existing links over total possible links) after 1983, because of the inclusion of 38 new country codes. Almost all of them are island countries, their data were aggregated in supranational entities before that year. As their contribution to the global trade volume is tiny, we drop those nodes for our analysis from 1983 to 2017.
After 1991, the reshape of the political maps of Eastern Europe and Central Asia produced a second fall of connectance event. It is impossible to follow the same approach because in this case the resulting political entities are as relevant as the Russian Federation. We keep former and new countries as separated entities for the global analysis (see Supplementary Information, Fig. S6).
www.nature.com/scientificreports www.nature.com/scientificreports/ To avoid the strongly meshed condition, a common solution is removing those links that account for a small fraction of the global trade 22,43 . We just drop those that fall within the lower 0.1% of total trade distribution for each year. The effect of the filtering is irrelevant for weight, as 99.9% of trade volume is kept, but not for degree 44,45 , because many tiny interactions are lost. As we are studying the aggregated figures of trade, we think this information loss is not critical for our analysis, but it is important to inform the reader about this issue.
After these preprocessing steps we create the weighted interaction matrix. It has one row per exporter country and one column per importer. We fill each cell with the global trade from exporter to importer or set it zero if there were no recorded trade or the volume falls below the filtering threshold. From this matrix, we obtain the number of exporters an importers and the number of links, the parameters that govern the synthetic network dynamics.
R packages. The statistical analyzes and their respective graphs have been carried out with different R package: Lilliefors and Kolmogorov-Smirnov tests with nortest 46 and densities estimation are plotted with the ggplot2 47 (geom_density() with bandwith adjustment = 2).

code and Data Availability
Code and data are available at https://zenodo.org/badge/latestdoi/120930457. Reproducibility instructions are detailed in the README.md file.