Louis Duvivier
email: louis.duvivier@insa-lyon.fr

Graph space: using both geometric and probabilistic structure to evaluate statistical graph models!

Keywords: statistical graph model, graph ensemble, entropy, statistical hypothesis testing, edit distance

Statistical graph models aim at representing graphs as random realization among a set of possible graphs. To evaluate the quality of a model M with respect to an observed network G, most statistical model selection methods rely on the probability that G was generated by M, which is computed based on the entropy of the associated microcanonical ensemble. In this paper, we introduce another possible definition of the quality of fit of a model based on the edit distance expected value (EDEV). We show that adding a geometric structure to the microcanonical ensemble induces an alternative perspective which may lead to select models which could potentially generate more different graphs, but whose structure is closer to the observed network. Finally we introduce a statistical hypothesis testing methodology based on this distance to evaluate the relevance of a candidate model with respect to an observed graph.

Introduction

The study of large and complex systems in domains as diverse as physics, biology, computer science or social science brings forward interaction networks which are neither regular, nor random. The interaction structure follows no clear pattern, as in a grid for example, but it often presents remarkable properties such as a heterogeneous degree distribution (scale-free) [START_REF] Barabási | Emergence of scaling in random networks[END_REF], small average distances (the small-world phenomenon) [START_REF] Watts | Collective dynamics of 'small-world'networks[END_REF], a high level of transitivity and community structure [START_REF] Newman | The structure and function of complex networks[END_REF], which are not present in random graphs. This suggests that interactions are the result of a random process following some constraints which shape the resulting network.

Various statistical models have been proposed to describe such constrained random graphs: the configuration model, which preserves degree distribution; the stochastic blockmodel, which preserves local densities; spatial models, which focus on distance between nodes [START_REF] Simini | A universal model for mobility and migration patterns[END_REF]. A review of statistical graph models can be found in [START_REF] Goldenberg | A survey of statistical network models[END_REF] and an introduction to their formalism in [START_REF] Cimini | The Statistical Physics of Real-World Networks[END_REF].

This diversity of parametric models raises the issue of model selection: considering one specific graph, which model is the most relevant? And with which set of parameters? Considering models as ensembles makes it possible to leverage results from statistical physics and information theory to perform model selection [START_REF] Ding | Model Selection Techniques -An Overview[END_REF], in particular using maximum likelihood estimators [START_REF] Stoica | Model-order selection: a review of information criterion rules[END_REF] or equivalently the minimum description length [START_REF] Grünwald | A tutorial introduction to the minimum description length principle[END_REF]. In this framework, a model is considered to be good for a given graph if it generates it with high probability, i.e. if the associated ensemble has a small entropy. Therefore, many papers have focused on computing the entropy of various graph ensemble [START_REF] Bianconi | Entropy of network ensembles[END_REF], [START_REF] Peixoto | Entropy of stochastic blockmodel ensembles[END_REF], [START_REF] Zingg | What is the Entropy of a Social Organization?[END_REF]. A whole methodology for community detection based on those principles has been developed in the case of stochastic blockmodels [START_REF] Peixoto | Bayesian stochastic blockmodeling[END_REF]. However, the definition of the likelihood of a model based on the entropy of its microcanonical ensemble implies that by itself this likelihood cannot be interpreted, it can only be compared with the likelihood of another model for the same graph: It is fundamentally a relative measure. More generally, this definition of the likelihood makes it possible to identify the model that generates the observed graph with the highest probability, but it is not always the best possible measure of how relevant a model is with respect to an observed graph. Indeed, graphs are also geometrical objects, in the sense that one can define distances between them. Such a distance induces a structure on a model's ensemble. Much work has been devoted to quantifying how similar two graphs are [START_REF] Wills | Metrics for Graph Comparison: A Practitioner's Guide[END_REF], especially from a topological point of view [START_REF] Monnig | The resistance perturbation distance: A metric for the analysis of dynamic networks[END_REF], [START_REF] Koutra | Deltacon: A principled massive-graph similarity function[END_REF]. These distances between two graphs can be generalized to evaluate the quality of a model by computing a distance between an observed graph and the graph ensemble associated to a model. For example, the widely used measure for community detection known as modularity [START_REF] Newman | Modularity and community structure in networks[END_REF] evaluates the quality of a partition by comparing the edge weight in the observed graph with the expected edge weight of the graph in the configuration model, as developed in section 2.2. In this case, the problem is to evaluate the statistical significance of the results, in order not to mistake noise for structure [START_REF] Guimera | Modularity from fluctuations in random graphs and complex networks[END_REF].

In this paper, after reviewing in section 2.1 and 2.2 existing techniques to measure the relevance of a model with respect to a graph, we introduce in section 3 the edit distance expected value, a measure which takes into account both the geometric and the probabilistic structure of the graph ensemble. Finally, we show how this measure can be used to evaluate a model relevance with respect to a given graph in section 4.

Microcanonical ensemble

The main point of statistical modeling is to describe the structure of a graph G based on a set of global properties of that graph P M [START_REF] Cimini | The Statistical Physics of Real-World Networks[END_REF]. Some examples of common models and their associated properties are given in table 1. Apart from those properties, the graph is considered to be random. In practice, this means that G is considered to have been chosen from the set:

Ω M = {H | P M (H) = P M (G)}
This set is called the microcanonical ensemble, in reference to statistical physics in which it was first introduced to represent all the possible states of a system corresponding to a global property. It can be defined with directed or undirected graphs, weighted or not, with or without self-loops. In the rest of the article, we will consider labelled directed multigraphs with self-loops. Although they are not the most widely used in practice, it makes probability derivations easier, especially for the configuration model. Information theory shows that, in order to ensure that no additional information bias the results, the probability distribution P M on Ω M has to be the one that maximises Shannon's entropy

S = -∑ H∈Ω M

P(H)log(P(H))

With no additional constraint, this optimal distribution is simply the uniform one:

∀H ∈ Ω M , P M [H] = 1 |Ω M | (2.1)
whose entropy is log(|Ω M |). Thus, computing the probability to choose G among all possible graphs in Ω M boils down to counting the number of graphs it contains. This has been done for different models in [START_REF] Peixoto | Entropy of stochastic blockmodel ensembles[END_REF].

The microcanonical ensemble and its entropy provide a common formulation for various statistical models. This is useful to perform model selection. Indeed, entropy is directly related to the likelihood of a given model. If we observe a graph G and consider a set of models M = {M 1 , . . . , M p }, we can find which model G has most likely been sampled from by maximising its likelihood

M * = argmax M i ∈M P[M i |G]
This maximisation can be done using Bayes theorem

P[M i |G] = P[G|M i ] × P[M i ] P[G] (2.2)
As P[G] does not depend on M i , maximising the likelihood is equivalent to maximising the numerator of equation (2.2). The first term of the product P[G|M i ] corresponds to the probability to generate G with the model M i , and according to equation (2.1), maximising it is equivalent to minimising the entropy of the associated microcanonical ensemble log(|Ω M i |). The second term P[M i ] is a prior distribution defined on the set of candidate models M . Its role is to account for the fact that a model with enough parameters can be made arbitrarily close to any given graph G, up to the point where Ω M = {G}. Such a model would generate G with probability P M [G] = 1, but it would be overfitting. The prior distribution prevents this from happening by assigning lower probabilities to models which have more parameters. Typical methods to counterweight models with too many parameters are the Akaike Information Criterion and Bayesian Information Criterion [START_REF] Stoica | Model-order selection: a review of information criterion rules[END_REF]. This idea was also developed and applied to the case of stochastic blockmodels in [START_REF] Peixoto | Bayesian stochastic blockmodeling[END_REF]. All these methods rely on a definition of the likelihood of a model according to the uniform distribution over the microcanonical ensemble, and thus on the hypothesis that the most relevant model M for a graph G is the one which generates G with the highest possible probability. It considers the microcanonical ensemble as a geometrically unstructured set, with no notion of distance between graphs, and thus it does not discriminate between two models M 1 and M 2 such that |Ω M1 | = |Ω M 2 | but where Ω M 1 contains graph similar to G while M 2 does not. This definition of the goodness of fit of a model is not always the best one. In particular, one can be interested in selecting the model that most probably generates graph similar to G, and not necessarily G itself.

For example, if we consider the three graphs on figure 1 with G 1 as a reference, both G 2 and G 3 are different from G 1 , but the topology of G 1 and G 2 is almost the same. Therefore, one could be interested in selecting a model which produces mostly G 2 -like graphs rather than one which produces G 3 -like ones. This cannot be done using the mere minimization of entropy since replacing G 2 by G 3 in the microcanonical ensemble does not change its entropy.

FIG. 1. Three graphs with 8 nodes and 10 edges. While G 2 is clearly more similar to G 1 than G 3 , both G 2 and G 3 are different from G 1 . Therefore, from the point of view of entropy maximisation, replacing G 2 by G 3 in the microcanonical ensemble does not affect the likelihood of a model.

Distance to the barycenter

On the other hand, one may measure how typical G is with respect to Ω M by comparing it with an appropriate representative of this ensemble, for example its barycenter:

G M = ∑ H∈Ω M

P(H)H

If we denote W G the weight matrix of graph G, it can easily be derived that The famous modularity function to evaluate the quality of a node partition B = (b 1 , . . . , b p ) on a graph G = (V, E) with weight matrix W G is defined as the difference between the number of edges

∀(i, j) ∈ V 2 ,W G M (i, j) = E[W H (i, j)] (2.
Model M Parameters G M (i, j) m G M Erdös-Renyi n, m m n 2 m Configuration model (k in i , k out i ) i∈V k out i k in j m m Stochastic block model (m r,s ) r,s m c i ,c j |c i ||c j | m Gravity model (k i ) i∈V , f k i k j f (d(i, j)) m Radiation model (k in i , k out i ) i∈V k out i k in i k in j (k in i +s i j )(k in i +k in j +s i j )
m inside each cluster and the expected number for a random graph with the same degree distribution (i.e. following the configuration model). It can be understood as a comparison with the barycenter G M of the corresponding configuration model.

Q(G, B) = 1 2m p ∑ i=1 ∑ u,v∈b i W G (u, v) - k out u k in v m = 1 2m p ∑ i=1 ∑ u,v∈b i (W G (u, v) -W G M (i, j)) = d(G, G M )
A problem is that G is compared with a single graph G M which is supposed to account for the whole graph ensemble Ω M . In particular, all information about the dispersion around the barycenter is lost, which undermines any attempt to interpret statistically the results.

Graph space and the edit distance expected value

As we have seen, existing techniques to compare a graph G and a model M exploit in different ways the ensemble Ω M . Entropy based techniques described in section 2.1 focus on its cardinality, but they neglect the topological similarities of graphs inside the ensemble. On the other hand, as described in section 2.2, an objective function such as the modularity accounts for these similarities, but it does so with a single graph which is supposed to represent the whole set. Reality is more complex: Ω M is a set of graphs with a probability distribution, and it can be further structured with a metric, making it a graph space. Both aspects, probabilistic and geometric, should be taken into account in order to understand the structure of Ω M , and the plausibility that a graph G was generated by the associated model M.

Many different measures exist to compute a similarity score between two graphs G and H [START_REF] Wills | Metrics for Graph Comparison: A Practitioner's Guide[END_REF]. One of the simplest is the edit distance. For two graphs on the same vertex sets G 1 = (V, E 1 ) and G 2 = (V, E 2 ), it counts the number of differences between their respective sets of edges.

ed(G 1 , G 2 ) = ∑ (i, j)∈V 2 |W G 1 (i, j) -W G 2 (i, j)|
As expected from its name, edit distance is a distance between graphs. Indeed, if we consider the weight matrix W G of a graph G as a point in R n 2 , the edit distance corresponds to the L 1 distance and for any model M, Ω M is a subset of R n 2 . The dimension prevents any direct drawing of it for graphs with more than 2 nodes, but it is possible to obtain some intuition about its shape.

In the following, we will use a normalized version of the edit distance which can be interpreted as the fraction of different edges between G 1 and G 2 .

ned(G 1 , G 2 ) = 1 2m ∑ (i, j)∈V 2 |W G 1 (i, j) -W G 2 (i, j)|
This normalized edit distance is no longer a distance on R n 2 . Yet, for all models M considered here, the number of edges m is constant over the set Ω M . Thus, the normalized edit distance is equivalent to edit distance inside Ω M and it allows to compare more easily results between various models, because whatever the model M, the distance between any two graphs G 1 and G 2 in Ω M is at most 1.

Edit distance to the barycenter

Ω M barycenter has already been introduced in section 2.2, where it was used as a proxy for the whole space. Using normalized edit distance, it is possible to check how much graphs in Ω M are similar to the barycenter G M . We consider six different models: For each model M, we pick a random sample S M of 100 graphs in Ω M and for all G ∈ S M we compute the normalized edit distance to the barycenter ned(G, G M ). Results are shown in figure 2.

The first thing to underline is that whatever the model, ned(G, G M ) is greater than 0.5, which means that most graphs in Ω M have at most half of their edges in common with G M . This observation shows that for those models, the graph space is not concentrated around its barycenter. On the contrary, most graphs in Ω M seem to be at a specific distance from its barycenter, as would happen for a sphere with a radius depending on the model: 0.67 for ER and CFM cst, 0.55 for CFM arith, 0.69 for SBM hom and 0.71 for SBM het.

All models were chosen to have similar entropy, as shown in table 3, yet their characteristic distance to the barycenter vary greatly. Furthermore, we observe that these quantities are not positively correlated: CFM arith, which is the model with the larger entropy is also the one which is the most concentrated around its barycenter. This means that even if this model can generate a higher number of different graphs, the graphs it produces tend to be more similar one to the other than for other models. This is logical as this model preserves a degree distribution, which enforces more constraints on edges' distribution than an Erdös-Renyi or a stochastic block model. This concentration of graphs at a specific distance from the barycenter is a consequence of the dimensionality of the vector space. Let's denote B(G, r) the ball of center G and radius r in (R n 2 , ed). We consider the set

Ω M (r) = {G ∈ Ω M | ned(G, G M ) r} = {G ∈ Ω M | ed(G, G M ) 2mr} = Ω M ∩ B(G M , 2mr) The volume V n (r) of B(G M , 2mr) is proportional to r n 2 , therefore ∀r < 1, V n (r) V n (1) -→ n→∞ 0 (3.1)
The volume of the ball concentrates quickly at its periphery as the dimension increases, and so does the volume of Ω M . The additional constraints on Ω M modify its shape in such a way that graphs too far away from the barycenter are rare, which explains why the concentration does not happen at distance 1 from the barycenter. Still, this phenomenon is strong enough to imply that even graphs generated according to a model M will share only a relatively small fraction of their edges with the barycenter of the model.

Edit distance expected value

The previous observations on the structure of graph spaces show that in order to compare a graph G with a model M, one should consider more than the mere cardinal of Ω M . One possibility to evaluate how similar to G are the graphs in Ω M is to compute the expected value of the normalized edit distance:

EDEV(G, M) = E H∈Ω M   1 2m ∑ (i, j)∈V 2 |W G (i, j) -W H (i, j)|  
To illustrate how EDEV provides further information on the place of G within the graph space, we compare it with entropy for different synthetic graphs. A low value indicates that G is close to other graphs in Ω M , and thus that it is typical of the model, while a high value shows that it is an outlier. As a case study, we consider the Erdös-Renyi model. Let's recall that we consider multigraphs, which implies that we allow for densities rising above 1. The extension of Erdös-Renyi model to multigraphs is straightforward, Ω ER(n,m) contains all multigraphs with n nodes and m edges and each multigraph is generated with the same probability 1 |Ω ER(n,m) | . In practice, we consider models with n = 100 nodes and a number of edges m ranging from 100 to 500000. For each, we consider three graphs:

• G 1 (m), picked uniformly at random inside Ω ER(n,m)
• G 2 (m), a graph made of two equal communities, each with n 2 nodes and m 2 edges, perfectly separated.

• G 3 (m), the graph where all edges are between nodes 0 and 1.

Results are shown on figure 3.

For each value of m, all three graphs belong to the graph ensemble Ω ER(n,m) . We observe that as density increases |Ω ER(n,m) | grows exponentially, which implies that the probability to pick at random G 1 (m), G 2 (m) or G 3 (m) becomes even less probable. Yet, in the case of the random graph G 1 (m) this is counter-intuitive: as density grows and becomes higher than 1, most graph in Ω ER(n,m) become complete graphs with each edge having weight about m n 2 . This is the case of G 1 too with a high probability, so ER(n, m) is very likely to produce graphs similar to G 1 (m), even if it is very unlikely to produce G 1 (m) itself.

On the other hand, edit distance expected value is able to capture this phenomenon. While it is close to 1 for all three types of graphs when density is low because in this situation a random model can hardly predict correctly which edge is present in any graph, it decreases quickly towards 0 when density rises above 0.1 for G 1 (m). For G 2 (m) we have an intermediate situation: edit distance decreases too, but it reaches its minimum around 0.5, indicating that even when it is densely populated, the model is only able to reproduce correctly half of its edges. This is normal as G 2 (m) concentrates them inside the communities, which means on half of all possible node pairs. These observations are actually a particular case of a more general result, which can be stated as:

∑ i, j∈[1,p 1 ] 2 M 1 [i, j] = ∑ k,l∈[1,p 2 ] 2 M 2 [k, l] = m
Let's consider two series of stochastic blockmodels defined as . This is what we observe with G 1 (m). Yet, we also observe on figure 3 that the normalized edit distance converges toward 0 only as density rises above 1. Thus, in practice, the vast majority of graphs are too sparse for this assumption to hold true and most graphs in Ω S(k) are far from G S(k) , as developed in section 3.1

S 1 (k) = (B 1 , k • M 1 ) and S 2 (k) = (B 2 , k • M 2 ), whose barycenters are denoted G 1 (k) and G 2 (k). We have that 1. There exists d ∈ R, ∀k ∈ N, ed(G 1 (k), G 2 (k)) = d 2. Let (G k ) k∈N

Statistical test

As the distance to the barycenter, the expected value of the normalized edit distance is characteristic of a model. For a model M, the values of EDEV(H, M) for graphs H in Ω M are concentrated around a specific value d M . We can use this fact to rule out models which fit badly on an observed graph G.

For example, let's consider the configuration model CFMD(n,

k out i , k in i ) n = 50 ∀i ∈ 0, n -1 , k out i = k in i = i
We use this model to generate a graph G i . CFMD(n, k out i , k in i ) will thus be called the generative model, and G i the observed graph. On this observed graph, we test the stochastic blockmodel SBM i obtained by partitioning its nodes in two blocks: B 0 contains even nodes and B 1 odd nodes (this way we avoid to put all high-degree nodes in the same block) and learning the block adjacency matrix on G i . We call SBM i the candidate model. We generate a sample S i of 100 test graphs with the candidate model SBM i and compare the normalized edit distance EDEV(G i , Ω SBM i ) of the observed graph to the candidate model with EDEV(H, Ω SBM i ) for all test graphs H ∈ S i . This experiment is performed 5 times, and results are shown on figure 4.

of 23

We observe that for all five experiments, the normalized edit distance expected value to the candidate model SBM i for the observed graph G i is around 0.74, while for the test graphs generated by SBM i it is between 0.53 and 0.56. This shows that the normalized edit distance expected value to SBM i is significantly different for the observed graphs, which are generated by CFMD(n, k out i , k in i ), and for the test graphs generated by SBM i . It is thus very unlikely that the observed graph G i was generated by the candidate model SBM i .

We then perform the same experiment the other way round, by considering as generative model a stochastic blockmodel SBM 0 (n, B, M) defined by Finally, as a null case, we consider the situation in which the generative model and the test model are the same. We generate 5 graphs (G i ) i∈ 1,5 with the same stochastic blockmodel SBM 0 (n, B, M) as before, and consider this same model SBM 0 as the candidate model. A sample S of 200 test graphs is randomly picked in Ω SBM 0 and we compare EDEV(G i , Ω SBM 0 ) with EDEV(H, Ω SBM 0 ) for all H ∈ S . Results are shown on figure 6. In this case, we observe that the normalized edit distance expected value for the first 5 graphs is not significantly different from the one of the sampled graphs. We cannot reject the hypothesis that the observed graph G i was generated by SBM 0 .

n = 50 B = 0,

Statistical hypothesis testing

The methodology can be formalized using statistical hypothesis testing. Let's say we have a graph G and a model M (possibly obtained by fitting some parameters on G). We want to test the null hypothesis H 0 : "The observed graph G has been generated by the candidate model M". If we knew the distribution of the edit distance expected value to the model M for graphs generated by M, we could simply compute EDEV(G, M) and evaluate the probability to obtain a value at least as large under the hypothesis that G was generated by M. However, characterizing the distribution P EDEV (M) for various models and deriving the parameters of the distribution from the parametrization of the model would require to further investigate the structure of their associated microcanonical ensembles.

As we do not know this distribution, we use bootstrapping to conduct our test. This approach consists in replacing the unknown distribution by a sample generated using this distribution. We therefore use M to generate a sample of q graphs G 1 , . . . G q . We can then consider G as a sample of size 1 from an unknown distribution M and H O can be reformulated as: "(G i ) i∈ 1,q and G were generated by the same probability distribution". Such an hypothesis can be tested using Fisher's permutation test, as described in [START_REF] Efron | An Introduction to the Bootstrap[END_REF]. We consider for each graph its edit distance expected value to the model M: x i = EDEV(G i , M) and y = EDEV(G, M), and compute the difference between the means of the two samples

θ = x -ȳ = 1 q q ∑ i=1 x i -y
We test H 0 by evaluating, under this hypothesis, the probability p = P H 0 [θ * θ ] to obtain a value of θ * at least as large for two samples of size 1 and q generated by M. The exact probability cannot be computed, as we do not know the probability distribution P EDEV (M), but it can be approximated by considering all the pairs of samples of size 1 and q that can be constructed by picking at random one graph within the q+1 graphs G, G 1 , . . . , G q . We denote P j the pair of samples {(G j ), (G 1 , . . . , G j-1 , G, G j+1 , . . . , G q )}.

For this pair of samples, the difference between the means is denoted

θ * j = 1 q ∑ i = j x i + y -x j
If H is true, the q + 1 pairs of samples P i are a subset of the possible pairs of samples generated by M, and we can make the approximation

P H 0 [θ * θ ] ≈ #{ j | θ * j θ } q + 1 If P H 0 [θ * θ ]
is lower than a preset threshold δ (in the following we will use δ = 0.01), one can conclude that there is enough evidence to reject the hypothesis that G was generated by M (and thus affirm that G was not generated by M). On the other hand, if it is greater than this threshold, there is not enough evidence to reject it. It should be stressed that this last sentence does not mean that one can affirm that G was generated by M, but only that this hypothesis cannot be discarded.

Let's consider a situation in which one wishes to evaluate the relevance of the block structure computed on a graph G, with a weight matrix W G . Whatever the graph, and the partition B = (b 1 , . . . , b p ) of its n nodes, it is always possible to define M ∈ M p (N) as n,B,M) . The objective is to evaluate whether this stochastic blockmodel is a relevant model of G. An even trickier question is to evaluate whether any stochastic blockmodel can be a relevant model. In particular, spatial models can generate graphs with groups of nodes densely connected due to their position rather than to block membership. An example of such a graph is shown in figure 7. It may then be hard to tell whether the blocks found are indeed a legitimate model of the observed graph or should be considered as artefacts, consequences of the underlying spatial structure.

∀i, j ∈ 1, p , M[i, j] = ∑ u∈b i ,v∈b j W G [u, v] such that G ∈ Ω SBM(
To illustrate how statistical hypothesis testing allows to address this issue, we consider eight models: four stochastic blockmodels and four Waxman models for random geometric graphs1 with different sets of parameters. The Waxman model for spatial graphs allows to easily control the strength of the spatial structure, by tuning the speed at which edge probability decays as the distance between nodes rises. The number of nodes is fixed to n = 100 and the parameters are fixed such as to ensure a density d around 0.036. All stochastic blockmodels use a node partition in four blocks of 25 nodes, with a block adjacency matrix of the form:

    m int m ext m ext m ext m ext m int m ext m ext m ext m ext m int m ext m ext m ext m ext m int    
The four stochastic blockmodels are then defined by:

1. M 0 : m int = 90, m ext = 0 2. M 1 : m int = 75, m ext = 5 3. M 2 : m int = 60, m ext = 10 4. M 3 : m int = 45, m ext = 15
This way, the graphs generated using SBM 0 are made of perfectly separated blocks of nodes, while those generated by SBM 3 have blocks with as many internal and external edges.

For the Waxman models, we also consider four parameter sets:

1. M 4 : α = 0.1, β = 1 FIG. 8.
For each model M i in ordinate, and each graph number i in abscissa, we compute the probability p i, j that a graph generated using SBM i, j has a normalized edit distance expected value to SBM i, j further away from the mean value d i, j than G i, j . The probability is plotted in heatmap. The lower the value of α, the stronger the spatial structure.

With each model M i , we generate 8 graphs (G i, j ) j∈ 0,7 . On each of those observed graph, we find the minimum entropy node partition B i, j using graph tools, and fit a candidate stochastic blockmodel SBM i, j on G i, j based on this node partition. We then evaluate the relevance of this stochastic blockmodel using the previously described methodology. We use a confidence level δ of 0.01 and a sample size q = 200. The probabilities p i, j = P H 0 [θ * i, j θ i, j ] obtained are plotted in figure 8.

We observe that for all graphs generated by stochastic blockmodels but three, p i, j > 10 -2 = δ . On the other hand, for all spatial graphs p i, j δ . This means that the hypothesis "G i, j has been generated by the candidate stochastic block model SBM i, j " is rejected for all 32 spatial graphs, and for three of the stochastic blockmodels generated graphs. Let's stress again that a probability p i, j superior to δ does not mean that SBM i, j is the right model for G i, j . It only means that there is not enough statistical evidence to reject it.

These results show that, on 32 spatial graphs generated with various sets of parameters, the statistical hypothesis testing methodology is able to correctly identify that the block structure found is not a relevant model, despite being the solution maximizing the likelihood among the candidate blockmodels. This result is not trivial. Indeed, the fact that this partition is the most likely according to the entropy definition of likelihood means that the block identified are characterized by a specific connection density which cannot be explained by random fluctuations. However, as the likelihood is maximized on a set of models which does not contain the Waxman model used to generate the graph, this estimator is not able to recognize that its structure would be even better explained using a spatial model. What is more, this methodology manages to reject the block structure for all spatial graphs while spuriously rejecting it for only three out of 32 graphs originally generated with a stochastic blockmodel. In other words, there is no false positive, and only three false negatives.

Strictly speaking, these results only allow to rule out one node partition. As the block structure tested were fitted on the observed graphs using minimum entropy, one could argue that ruling out this partition implies that no other node partition can lead to a relevant model. Yet, for most real graphs, there is more than one plausible node partition and minimizing the entropy of a partition is a stochastic process. Therefore, the experiment should be performed more than once to conclude that the observed graph has no block structure.

Conclusion

As a conclusion, we have seen that widely used quality measures for graph models rely either on the number of different graphs they can produce, which neglects the geometric structure of the graph space, or on a direct comparison with the barycenter of those graphs, which discards information about the distribution around this barycenter. Because of these restrictions, they are unable to distinguish between graphs which have a typical structure of a model and graphs which may be generated by this model but as outliers.

In this paper, we show how edit distance can provide additional information on the structure of the graph space which is captured neither by the entropy nor the barycenter. By computing the expected value of the normalized edit distance for a given graph, we obtain a criterion which can be used to evaluate the model quality with respect to this graph. Finally, we incorporate this criterion to a statistical hypothesis testing methodology to perform model selection.

Graph space is a theoretical framework which can be used for any statistical model, and particularly spatial models. It allows to compare them with SBM or configuration model, and perform model selection between models of different nature. What is more, statistical hypothesis testing provides a statistically rigorous methodology to evaluate the relevance of a candidate model to an observed graph.

Apart from its simplicity to compute and interpret, an interesting result about the normalized edit distance expected value is that the same quantity can be used to test many different models. Statistical tests could be performed on any graph property, such as the average path length or the clustering coefficient, to rule out a candidate model. Yet, such tests require to choose for each model and each tested graph a property that could be used as a statistical test. In turn, all these properties can be used to define new models which better reproduce the values measured on the observed graph (as done for example when defining degree-corrected stochastic blockmodels, which incorporate the degree distribution constraint to the block structure). The risk being that by taking more and more properties into account, one eventually overfits the graph. Note, however, that the expected value of the normalized edit distance is somewhat different from the previously mentioned properties. Indeed, as it depends on both the graph and the whole microcanonical ensemble, it is impossible to compute its value on a graph observed beforehand and then to define a model adapted to this specific value.

Indeed, the fact that, for most models, the normalized edit distance expected value distribution quickly clusters around the mean, is a property of the graph space rather than the graph itself. It is a consequence of a geometrical result (the volume of a ball in n dimensions), which highlights the benefits of considering the geometric structure of graph ensembles. However, the edit distance is not the only distance that can be used. Considering other metrics which are more sensitive to the global topology of the network, like the perturbation-resistance distance or spectral distances, could provide additional insight on the structure of the graph space.

Appendix

Barycenter graph weight of various statistical models

We have defined the barycenter of a graph model as

G M = ∑ H∈Ω M P(H) • H Which means that ∀(i, j) ∈ V 2 ,W G M (i, j) = ∑ H∈Ω M P(H) ×W H (i, j) = E[W H (i, j)]
Let's illustrate how this can be computed for some classical models.

ERD ÖS-R ÉNYI MODEL The simplest graph model is the Erdös-Rényi model for random graphs. It's associated microcanonical ensemble can be defined as:

Ω ER(n,m) = {H = (V, E) | |V | = n ∧ |E| = m}
Let's recall that for the sake of simplicity, we chose to consider multigraphs with self loops. Thus, the computation of E[W H (i, j)] is particularly simple. Indeed, if for each pair of node (i, j) ∈ V 2 and each k ∈ [1, m] we define the random variable X i, j,k which is equal to 1 if the k th edge is i → j and to 0 else, then we have that W H (i, j) = ∑ m k=1 X i, j,k . It is a sum of independent Bernouillis' random variable so it follows a binomial law of parameters m and 1 n 2 , and thus

W G ER(n,m) (i, j) = E Ω ER [W H (i, j)] = m n 2 (6.1)
CONFIGURATION MODEL For the configuration model, all graphs in the microcanonical ensemble must have the same degree distribution. Let's consider the directed version.

Ω CFMD = {G | ∀i ∈ V, deg out G (i) = k out i ∧ deg in G (i) = k in i }
To compute the weight of the barycenter graph's edges, we consider that each node i has k out i outgoing stubs and k in i ingoing stubs. Any graph in Ω CFMD is characterized by a configuration of connections of outgoing stubs with ingoing stubs. For every pair of nodes i, j ∈ V 2 and any pair of stub

k ∈ [1, k out i ], l ∈ [1, k in j ],
we define the random variable X i, j,k,l which is equal to 1 if the k th outgoing stub of i is connected to the l th ingoing stub of j, and to 0 otherwise. Then

W H (i, j) = k out i ∑ k=1 k in j ∑ l=1 X i, j,k,l
As each outgoing stub of i has the same probability to be connected to any of the m ingoing stubs ∀i, j, k, l, P[X i, j,k,l = 1] = 1 m Thus, W H (i, j) follows a binomial law of parameters 1 m and k out i × k in j . Finally

W G CFMD (i, j) = E Ω CFMD [W H (i, j)] = k out i × k in j m (6.2)
STOCHASTIC BLOCKMODEL The case of the stochastic blockmodel can be treated in the same way as erdös-rényi. It is defined, considering a partition of the nodes B = (b 1 , . . . , b q ) and a block adjacency matrix M ∈ M q (N) by

Ω SBM = H | ∀b k , b l , ∑ i∈b k ∑ j∈b l W H (i, j) = M(k, l) So, for any pair of nodes i ∈ b k , j ∈ b l , W H (i, j) follows a binomial law of parameters (M(k, l), |b k ||b l |). Thus W G SBM (i, j) = E Ω SBM [W H (i, j)] = M(k, l) |b k ||b l | (6.3)
SPATIAL MODELS References for the gravitational model and the radiation model can be found in [START_REF] Barthélemy | Spatial networks[END_REF] and [START_REF] Simini | A universal model for mobility and migration patterns[END_REF]. In both cases, they are constructed in such a way that edges weight have a given expected value. In the case of the gravitational model, it is

W G grav (i, j) = f(d(i, j)) × k out i × k in j (6.4)
where d(i, j) is the distance from node i to node j, and f is a deterence function.

Finally, in the case of the radiation model, it is and W G k (u, v) ∼ B(k • M i, j , p i, j ) with p i, j = 1 |b i ||b j | . Therefore, according to the Bienaymé-Tchebychev inequality:

W G rad (i, j) = k out i × k in i × k in j (k in i + s i j ) × (k in i + k in j + s i j ) (6.5) with s i j = ∑ u∈C (i, j) k in u and C (i, j) = {u ∈ V | 0 < d(i, u) < d(i, j)}.
P W G k (u, v) -W G S(k) 2km > α n 2 k × M i, j × p i, j × (1 -p i, j ) × n 2 4 × k 2 × m 2 × α M i, j × p i, j × (1 -p i, j ) × n 2 4 × k × m 2 × α -→
W G 1 (k) [i, j] -W G 2 (k) [i, j] = 1 2km ∑ i, j∈[1,n] 2 k • M 1 [b 1 (i), b 1 ( j)] |b 1 (i)||b 1 ( j)| -k • M 2 [b 2 (i), b 2 ( j)] |b 2 (i)||b 2 ( j)| = 1 2m ∑ i, j∈[1,n] 2 M 1 [b 1 (i), b 1 ( j)] |b 1 (i)||b 1 ( j)| - M 2 [b 2 (i), b 2 ( j)] |b 2 (i)||b 2 ( j)|
which is constant with respect to k. In the following we will denote this distance d for the sake of conciseness. We want to show that EDEV(G k , S 2 (k))

1 .

 1 EM: Erdös-Renyi with 50 nodes and 1000 edges 2. CFM cst: configuration model with 50 nodes and a constant degree distribution (k in i = k out i = 20) 3. CFM arith: configuration model with 50 nodes and an arithmetic degree distribution (k in i = k out i = i + 1) 4. SBM hom: stochastic block model with 50 nodes and 5 communities, each having internal density 1.2, and external density 0.2. 5. SBM het: stochastic block model with 50 nodes and 5 communities, with internal density 0.4, 0.8, 1.2, 1.6, 2, and external density 0.2.

FIG. 2 .

 2 FIG.2. Edit distance to the barycenter for 6 different models. For each model M in ordinate, we draw 100 graphs G at random from Ω M and compute for each of them ned(G, G M ). The distribution of results is then plotted as a boxplot.

LEMMA 3 . 1 FIG. 3 .

 313 FIG. 3. Entropy and edit distance expected value against density. For each density, three graphs are generated. G 1 (m) is random, G 2 (m) is made of two random communities and G 3 (m) has its edges concentrated on a single pair of nodes of weight m. On the top plot, the entropy log(|Ω ER(n,m) |) is plotted against the density m n 2 . As all graphs belong to the same graph ensemble Ω ER(n,m) , the three curves are the same. On the bottom plot the edit distance expected value EDEV (G(m), Ω ER(n,m) ) is plotted against density.

  FIG.4. 5 graphs G i were generated with a generative model CFMD(n, k out i , k in i ). Their normalized edit distance expected value with respect to a candidate model SBM i is plotted as diamond. As a comparison point, the distribution of the normalized edit distance expected value for 100 test graphs generated with SBM i is plotted as dots and boxplot.

FIG. 5 .

 5 FIG.5. 5 graphs G i were generated with a generative model SBM(n, B, M). Their normalized edit distance expected value with respect to a test model CFMD i is plotted as diamond. As a comparison point, the distribution of the normalized edit distance expected value for 100 test graphs generated with CFMD i is plotted as dots and boxplot.

FIG. 7 .

 7 FIG. 7. Spatial graph generated using the Waxman random geometric model. 100 nodes are randomly distributed in a [0, 1] × [0, 1] square. They are then connected with a probability depending on the distance d between nodes: p(d) = 10 exp d 0.05L , with L the maximum distance between two nodes. Communities are computed using graph-tool (https://graph-tool.skewed. de/).

2. M 5 :

 5 α = 0.08, β = 1.6 3. M 6 : α = 0.06, β = 2.7 4. M 7 : α = 0.04, β = 8.5

20 of 23 6 . 2 > α n 2 andP> α n 2 Let

 206222 Convergence proofFirst of all, let's prove the following lemma LEMMA 6.1 Let B be a partition of 1, n with p blocks. Let M ∈ M p (N) be a block adjacency matrix. For all k ∈ N, we define the stochastic blockmodel S(k) = (B, k • M), and its barycenter G S(k) . We consider a sequence of random graphs (G k ) k∈N , each drawn from S(k). We have that ed(G k , G S(k) ) above, we want to prove that:∀α > 0, P[ed(G k , G S(k) ) > α] -→ k→∞ 0Let α > 0. Let's denote m = ∑ i, j M i, j the number of edges of graphs in Ω S[START_REF] Barabási | Emergence of scaling in random networks[END_REF] . By definition,ed(G k , G S(k) ) = 1 2km ∑ u,v |W G k (u, v) -W G S(k) (u, v)| Thus, ed(G k , G S(k) ) > α ⇒ ∃(u, v), W G k (u, v) -W G S(k) (u, v) 2km [ed(G k , G S(k) ) > α] ∑ u,v P W G k (u, v) -W G S(k) (u, v)2km's consider two blocks b i and b j in B. We know that ∀u ∈ b i , v ∈ b j ,W G S(k) (u, v) = k • M i, j |b i ||b j |

6 ) 2 M 1 2 M 2

 62122 (G k , G S(k) ) > α]Which proves the lemma.We can now prove the theorem THEOREM 6.1 Let B 1 and B 2 be two partition on 1, n , with p 1 and p 2 blocks respectively. Let M 1 ∈ M p 1 (N) and M 2 ∈ M p 2 (N) be two block adjacency matrices such that∑ i, j∈[1,p 1 ] [i, j] = ∑ k,l∈[1,p 2 ] [k, l] = mLet's consider two series of stochastic blockmodels defined asS 1 (k) = (B 1 , k • M 1 ) and S 2 (k) = (B 2 , k • M 2 ), whose barycenters are denoted G 1 (k) and G 2 (k). We have that 1. There exists d ∈ R, ∀k ∈ N, ed(G 1 (k), G 2 (k)) = d 2. Let (G k ) k∈N be series of random graph each drawn following model S 1 (k). EDEV(G k , S 2 (k))For any pair of nodes i, j, belonging to blocks b(i) and b( j) in B 1 (resp. B 2 ), the weight of the edge i → j in G 1 (k) (resp. G 2 (k)) is given by:W G 1 (k) [i, j] = k • M[b(i), b( j)] |b(i)||b( j)| Therefore, the edit distance between G 1 (k) and G 2 (k) is ed(G 1 (k), G 2 (k)) = 1 2km ∑ i, j∈[1,n] 2

  We start by noticing thatEDEV(G k , S 2 (k))d = E H∈S 2 (k) [ed(G k , H)]d E H∈S 2 (k) [ed(G k , G 1 (k)) + ed(G 1 (k), G 2 (k)) + ed(G 2 (k), H)]d ed(G k , G 1 (k)) + E H∈S 2 (k) [ed(G 2 (k), H)]On the other hand,d -EDEV(G k , S 2 (k)) = E H∈S 2 (k) [ed(G 1 (k), G 2 (k))ed(G k , H)] E H∈S 2 (k) [ed(G 1 (k), G k ) + ed(G k , H) + ed(H, G 2 (k))ed(G k , H)] ed(G k , G 1 (k)) + E H∈S 2 (k) [ed(G 2 (k), H)] Thus |EDEV(G k , S 2 (k)) -d| ed(G k , G 1 (k)) + E H∈Ω S 2 (k) [ed(G 2 (k), H)](6.7)Because G k is generated following S 1 (k), a direct application of lemma 6.1 is thated(G k , G 1 (k)) , if H is generated following S 2 (k), we also have that ed(H, G 2 (k)) obtain that ed(G k , G 1 (k)) + E H∈Ω S 2 (k) [H, G 2 (k)]

Table 1 .

 1 Common statistical graph models and their associated properties.

	Model M	Properties P M
	Erdös-Renyi	number of nodes and number of edges
	Configuration model	degree distribution
	Stochastic block model	block to block density
	Gravity model	node position, strength, and deterrence function
	Radiation model	node position and strength
	2.1 Entropy

Table 2 .

 2 Statistical graph models' barycenter weight.

Table 3 .

 3 Edit distance to the barycenter and entropy.

	Model	Characteristic distance	Entropy
	ER	0.67	2050
	CFM cst	0.67	2500
	CFM arith	0.55	3300
	SBM hom	0.69	1840
	SBM het	0.71	1840

https://networkx.org/documentation/stable/reference/generated/networkx.generators. geometric.waxman\_graph.html

Funding

Work supported by the ACADEMICS grant of the IDEXLYON, project of the Université de Lyon, PIA operated by ANR-16-IDEX-0005 and the BITUNAM grand ANR-18-CE23-0004.


Célestin Coquidé

José Lages
email: jose.lages@univ-fcomte.fr

Dima L Shepelyansky

Dollar-Yuan Battle in the World Trade Network

Keywords: International trade, Commercial flows, Markov chains, Complex networks

From the Bretton Woods agreement in 1944 till the present day, the US dollar has been the dominant currency in the world trade. However, the rise of the Chinese economy led recently to the emergence of trade transactions in Chinese yuan. Here, we analyze mathematically how the structure of the international trade flows would favor a country to trade whether in US dollar or in Chinese yuan. The computation of the trade currency preference is based on the world trade network built from the 2010-2020 UN Comtrade data. The preference of a country to trade in US dollar or Chinese yuan is determined by two multiplicative factors: the relative weight of trade volume exchanged by the country with its direct trade partners, and the relative weight of its trade partners in the global international trade. The performed analysis, based on Ising spin interactions on the world trade network, shows that, from 2010 to present, a transition took place, and the majority of the world countries would have now a preference to trade in Chinese yuan if one only consider the world trade network structure.

Significance Statement -Taking only into account the network structure of the commercial flows, we determine whether a country would prefer to trade in US dollar or in Chinese yuan. The evolution of this trade currency preference is obtained along the last decade: The fraction of countries preferring firmly to trade in Chinese yuan increased from 34% in 2010 to 60% in 2020, and concerns at present all the Asian countries, Australia, the former Soviet Union, and almost all African and Latin American countries. European countries are also on the brink of a transition toward Chinese yuan preference. Our mathematical analysis focuses on the economical aspects of the international trade neglecting other purely geopolitical aspects.

In 1944, the Bretton Woods agreements established a system of payments based on the US dollar (USD). The USD effectively became the world currency, i.e., the standard to which every other currency was pegged [START_REF]Bretton Woods system -Wikipedia, The Free Encyclopedia[END_REF]. Till now, the USD remained the dominant world trade currency and, as an example, the United Nations (UN) reports the world trade transactions between countries in USD [2]. However, the possible end of the dollar dominance is increasingly discussed as the Chinese yuan (CNY) is gradually becoming credible as a reserve currency [START_REF] Raisinghani | Could China's Yuan replace the U.S. dollar as the world's dominant currency?[END_REF][START_REF] Ahmed | Can yuan replace the 'mighty' dollar?, The Business Standard[END_REF][START_REF] Amadeo | How the Yuan Could Become a Global Currency, The balance[END_REF][START_REF] Curran | Dollar's Dominance Is Being Stealthily Eroded[END_REF][START_REF] Liu | Can BRICS De-dollarize the Global Financial System?[END_REF][8]. Moreover, recently important trade transactions were considered to be realized in CNY instead of USD like oil sales from Saudi Arabia to China [START_REF] Said | Saudi Arabia Considers Accepting Yuan Instead of Dollars for Chinese Oil Sales[END_REF]. Thus, the growth of the Chinese economy [START_REF]China Overview: Development news, research[END_REF] opens up the possibility for a country to prefer the CNY to the USD for their trade exchanges. In the present paper, we use a mathematical model mimicking the preference of world countries to trade in USD or CNY. This currency preference model is applied on the world trade network (WTN) built from the 2010-2020 UN Comtrade data [2]. Our analysis, based purely on trade relations without taking account of external political factors or other considerations, clearly shows that, during this decade, a USDto-CNY transition took place with the result that the main part of the world countries would now prefer to trade in CNY instead of USD. The WTN is constructed from the UN Comtrade database [2] which provides the money matrix M encoding the transactions of all the commodities between all the countries. The money matrix element M c ′ c gives the total amount of commodities, expressed in USD of a given year, exported from the country c to the country c ′ . The UN Comtrade database concerns 194 countries for the period 2010-2020. The Markov chain of trade transactions is characterized by two WTN trade matrices S and S * whose elements are S c ′ c = M c ′ c /M * c and S * c ′ c = M cc ′ /M c . Here,

M * c = c ′ M c ′ c (M c = c ′ M cc ′ )
gives the total amount of commodities exported from (imported to) the country c to (from) the rest of the world. Consequently, the matrix element S c ′ c (S * c ′ c ) gives the relative weight of the imports to (exports from) the country c ′ from (to) the country c. These matrix elements can be considered as link weights of a directed network (see e.g. [START_REF] Dorogovtsev | Lectures in Complex Networks[END_REF]). By construction, the sum of each column of the matrices S and S * is equal to 1, i.e., c ′ S c ′ c = 1 and c ′ S * c ′ c = 1, ensuring the fact that the stochastic processes of the trade transactions, both in export and import directions, belong to the class of Markov chains. Let us define, for a given year, the total world trade volume by M = c M c = c M * c . Then, a given country c can be characterized by its import trade probability P c = M c /M and its export trade probability P * c = M * c /M . These probabilities characterize the global capability of a given country to import and to export, respectively. As an example, in 2019, the top 5 countries according to the import trade probability P c are: 1. USA, 2. China, 3. Germany, 4. Japan, [START_REF] Blondel | Fasdt unfolding of communities in large networks[END_REF] with Dugué's algorithm [START_REF] Dugué | Directed Louvain : maximizing modularity in directed networks[END_REF]. Each cluster is labelled using the ISO2 code of the country with the best import and export trade probabilities Pc and P * c . The leaders of the clusters are the USA (blue), China (red), and Germany (green). The cluster Others (gold) gathers other small size clusters. and the top 5 Countries according to the export trade probability P * c are: 1. China, 2. USA, 3. Germany, 4. Japan, 5. France.

In order to determine the trade currency preference (TCP) of the country c, i.e., whether the country c would prefer to trade in USD or in CNY with the other world countries, we assign to the country c an Ising spin σ c = ±1. The value σ c = -1 (+1) indicates that the country c prefers to trade in USD (CNY). Thus, we obtain a network of interacting Ising spins (an Ising spin is attached to each country). The interaction energy of the country (spin) c is

E c = c ′ ̸ =c σ c ′ (S c ′ c + S * c ′ c ) (P c ′ + P * c ′ ) 2 . ( 1
)
which is the sum of all the energies of interaction of the country (spin) c with its direct trade partners (spins) c ′ . Each element of this sum is the product of two factors: the relative trade import-export flows between the country c and a direct trade partner c ′ (terms in the first bracket of (1)), and the global trade capability of the partner c ′ in the WTN (terms in the second bracket of (1)). The term σ c ′ in (1) indicates the preference of the commercial partner c ′ to perform trade in USD or CNY. Consequently, depending on the interactions with its commercial partners, we allow a country c to change its TCP: If E c < 0 then the country c spin becomes σ c = -1 (preference to trade in USD) and if E c > 0 then the country c spin becomes σ c = +1 (preference to trade in CNY). We keep USA and China always trading in USD and CNY respectively. We start with an initial fraction f i of randomly chosen countries which prefer to trade in USD, we assign a -1 value to their spins. Consequently, the complementary fraction 1 -f i of countries prefer to trade in CNY, we assign a +1 value to their spins. A first Monte Carlo shake allows to determine the energy E c1 for a randomly picked spin σ c1 . The spin σ c1 flips or not according to the sign of the newly computed E c1 , and consequently the TCP of the country c 1 either stay the same or possibly changes from one currency to the other (i.e., from USD to CNY or from CNY to USD). Then, with the obtained new configuration of spins, a second shake is performed for another randomly chosen spin σ c2 (c 2 ̸ = c 1 ), and so on. After 192 shakes, the trade currency preference for each country is determined (USA and China spins are always kept fixed). The ensemble of these shakes forms the first time step τ = 1. We observe that after five consecutive time steps (τ > 5) the system converges to a fixed steady-state configuration of spins which stays unchanged for higher values of τ . This procedure is applied to N r = 10 4 initial random spin configurations with a fixed fraction f i of countries with initially a TCP for USD. We follow the evolution with τ of the fraction f (τ ) of countries preferring to trade in USD till the steady-state f f at τ = 10 is obtained. The complementary fraction 1 -f f gives then the fraction of countries preferring to trade in CNY once the steadystate is reached. An example of time evolution f (τ ) is shown in Fig. S1 of the Supplementary Information. We note that a similar approach has been used in [START_REF] Kandiah | PageRank model of opinion formation on social networks[END_REF] for the opinion formation on social networks.

RESULTS

As a preliminary step, let us characterize, for the years 2010 and 2019, the WTN using the Louvain method for cluster detection [START_REF] Blondel | Fasdt unfolding of communities in large networks[END_REF][START_REF] Dugué | Directed Louvain : maximizing modularity in directed networks[END_REF]. The results presented in Fig. 1 show the existence of 3 main clusters formed around the USA, China and Germany. We observe that from 2010 to 2019 the size of the cluster around China extends significantly. Indeed, from 2010 to 2019, the US-cluster loses almost all of the South American countries to the benefit of the CN-cluster, and a dominant part of Africa enters the CN-cluster. Meanwhile, the cluster formed around Germany remains practically unchanged including the EU countries, the countries of the former Soviet Union, and most of the Maghreb. Examples of WTN clustering for other years of the considered decade are shown in Fig. S2 of the Supplementary Information.

However, this preliminary cluster analysis of the WTN, based on the maximization of the modularity [START_REF] Blondel | Fasdt unfolding of communities in large networks[END_REF][START_REF] Dugué | Directed Louvain : maximizing modularity in directed networks[END_REF], does not determine the trade preference of the countries either for USD or CNY. The Monte Carlo procedure,

f i ρ (f i ) f f2 FIG. 2.
Final fraction f f of countries preferring to trade in USD as a function of the initial fraction fi for the years 2010 (left panel) and 2019 (right panel). Left and right panels: for any initial fraction fi, the Monte Carlo procedure converges toward one of two final fractions f f . These two final fractions are f f 1 = 0.16 and f f 2 = 0.66 in 2010 (left panel) and f f 1 = 0.14 and f f 2 = 0.44 in 2019 (right panel). The color of a point (fi,f f ) indicates the portion ρ f f (fi) of the Nr = 10 4 initial configurations, with a corresponding initial fraction fi, which attain the final state with the corresponding final fraction f f . The color ranges from black for ρ f f (fi) = 0 (all the countries preferring to trade in CNY whether than in USD) to bright yellow for ρ f f (fi) = 1 (all the countries preferring to trade in USD whether than in CNY). Middle panel: Portion ρ f f (fi) of the Nr = 10 4 initial configurations, the fraction fi of which initially prefers to trade in USD, which attains the final state with the final fraction f f . The red (blue) curve and the up (down) triangles correspond to the lowest (highest) value f f 1 (f f 2 ) of the two final fractions f f . The full (empty) symbols correspond to the year 2019 (2010).

described in the previous section, allows to obtain the final fraction f f of world countries which prefers to trade in USD, and conversely the final fraction 1 -f f of world countries which prefers to trade in CNY. Fig. 2 shows the final fraction f f as a function of the fraction f i of countries which initially prefer to trade in USD. For each value of f i , we randomly picked N r = 10 4 different configurations of spins. We observe that for any initial fraction f i belonging to the interval [0, 1], only two final fractions f f1 and f f2 can be reached. However, the probability that a given spin configuration reach one or the other final fraction values depends on the initial distribution of the TCPs over the countries. Let us take f f1 < f f2 . Quite naturally, higher (lower) is the initial fraction f i , higher is the probability to obtain the highest (lowest) final value f f2 (f f1 ). The middle panel of Fig. 2 gives, for the years 2010 and 2019, the probabilities ρ f f 1 (f i ) and ρ f f 2 (f i ) to obtain the final fraction f f1 and f f2 as a function of the initial fraction f i . In 2010, see left panel of Fig. 2, each of the final fractions corresponded to a majority of countries with either a USD preference (f f2 = 0.66) or a CNY preference (f f1 = 0.16). This no longer the case in 2019, see right panel of Fig. 2, for which the two final fractions f f1 = 0.14 and f f2 = 0.44 are below 0.5 and give both a CNY preference for the majority of the world countries. In one decade and according to the sole structure of the WTN, we pass from a bipolar USD-CNY trade currency preference to a global domination of the CNY.

Let us define the TCP probability P $ (c) to obtain for FIG. 3. World distribution of the trade currency preference probability for the years 2010 and 2019. Each panel corresponds to a given year and a given fraction fi of countries preferring to initially trade in USD. Left (right) column panels correspond to the year 2010 (2019). The top, middle, and bottom rows correspond to fi = 0.1, 0.5, and 0.9, respectively. Each country is characterized by a probability P $ to obtain a USD trade preference at the end of the Monte Carlo procedure. The CNY trade preference probability of a country is then P ¥ = 1 -P $ . The color ranges from red for P $ = 0 and P ¥ = 1 to blue for P $ = 1 and P ¥ = 0. The average TCP probability P $ has been computed from Nr = 10 4 random initial TCP distributions.

the country c a USD preference at the end of the Monte Carlo procedure. The probability to obtain for the country c a CNY preference is then P ¥ (c) = 1 -P $ (c). The probability P $ is obtained from the application of the Monte Carlo procedure to the N r = 10 4 initial TCP distributions. Fig. 3 shows the TCP probability world distribution. The Fig. 3 left panels illustrate the above described bipolar USD-CNY trade currency preference which exists in 2010: for high (low) f i most of the countries finally prefer USD (CNY). In 2019, Fig. 3 right panels, the CNY dominance is clearly observed. Indeed, even for high f i , most of the countries finally prefer CNY over USD.

The reason of the bistability of the final outcomes f f (see f 1 and f 2 in Fig. 2) can be understood from the analysis of the distribution of the USD or CNY trade preference over the countries. In fact, there are two groups of countries which keep, for any initial fraction f i , a hard preference to trade in USD (the USD group) and in CNY (the CNY group). Otherwise stated, a country of the CNY (USD) group, independently of its initial TCP and of the initial TCPs of the other countries, will always ends up in the CNY (USD) group. A third group of swing states (the swing group) may change their TCP. Amazingly, depending on the initial configuration of countries which prefer to trade in USD or in CNY, the countries belonging to this swing group collectively adopt at the equilibrium either the USD or the CNY as trade currency. This swing states are collectively responsible for the final outcome: if they adopt USD (CNY) at the equilibrium, the final fraction f f will be the highest (lowest) of the two possible fractions, i.e., f 2 (f 1 ). These three groups are shown on the world maps displayed in Fig. 4 for the years 2010 and 2019 (the world maps for the other years of the past decade are shown in Fig. S3 of the Supplementary Information). We clearly see that there is a drastic change from 2010 to 2019: a large number of countries passed from the swing group to the CNY group which has considerably increased in size during the last decade. Indeed, the former Soviet Union countries, almost all South America and Africa belong in 2019 to the CNY group. By contrast, from 2010 to 2019, the size of the USD group reduced only slightly loosing few countries: Venezuela, Nigeria, Chad, South Sudan, Equatorial Guinea, Afghanistan and Federated States of Micronesia switched to the CNY group and, Suriname and Israel to the swing group. In 2019, the swing group is mainly composed by EU countries, the UK and some Mediterranean countries (Turkey, Egypt, Morocco, Algeria, Tunisia and Israel). The lists of the countries belonging to the USD, CNY and swing groups in 2010 and in 2019 are given in Tables S1 to S6 of the Supplementary Information.

The time evolution over the last decade of the size of the three TCP groups is shown in Fig. 5. The left panel of Fig. 5 displays the fraction of countries belonging to The width of a band corresponds to the size of the corresponding group expressed as: (left panel) the ratio between the countries belonging to the group over and the total number of countries, (right panel) the ratio between the total trade volume exchanged by the countries of the group and the total volume exchanged by all the world countries. The trade volumes are expressed in USD of the concerned year.

each group and the right panel displays the fraction of the total volume of import and export exchanged by each group. From Fig. 5 left panel, we observe that the CNY group (red band) steadily grows along the decade from a fraction of 34% of the countries in 2010 to 60% in 2020. This growth of the CNY group is mainly compensated by the depletion of the swing group (green band) whose size drops from 50% of the countries in 2010 to 27% in 2020. Meanwhile, the size of the USD group (blue band) slightly decreased from 16% of the countries in 2010 to 13% in 2020. The trends are the same but less pronounced for the fraction of the trade volume exchanged by the different groups (see Fig. 5 Let us note that we obtain practically the same results if, instead of keeping China always trading in CNY and the USA always trading in USD, we keep China and the other BRICS (Brazil, Russia, India, South Africa) always trading in CNY and the USA and other Anglo-Saxon countries (Canada, UK, Australia, New Zealand) always trading in USD; see e.g., Figs. S4 and S5 in the Supplementary Information which are quite similar to Figs. 2 and4. This result asserts the dominance of China and USA in the world trade network. We have also considered to replace the import trade probabilities P c and the export trade probabilities P * c in (1) by the PageRank and the CheiRank probabilities obtained from the Google matrix of the WTN. These probabilities allowing to measure the capability of a country to import or export products throughout the WTN were used to analyze the international trade [START_REF] Ermann | Google matrix of the world trade network[END_REF]16]. Such a replacement of the probabilities, i.e., of the centrality measures of the WTN, leads again to practically the same results (see e.g., Figs. S6 and S7 in Supplementary Information which are similar to Figs. 2 and4).

CONCLUSION

The question addressed in the current work is the following one: assuming that only the WTN structure matters, what would be the trade currency preference for each country? Our analysis, performed by superimposing a Ising spin network on the WTN, clearly shows that the main part of the world would prefer to trade nowadays in CNY while in 2010 it would have preferred to trade in USD. We observe two final equilibrium states. In 2010, one of them characterizes an USD preference and the other one a CNY preference, whilst in 2019, both of the two final states characterize a CNY preference. Nowadays, according to the WTN structure, for any initial distributions of countries preferring to trade in USD and in CNY, the final state would always favor a world which preferentially trade in CNY. The bistability of the final state is due to a group of swing states which, depending on the initial distribution of the trade currency preferences over the countries, adopt all together a preference for either USD or CNY. Of course, our analysis is based on the mathematical treatment of the trade flows between the world countries and does not take into account any geopolitical relations between the countries. But, it is often claimed that economics determines politics and thus we argue that the obtained results demonstrate a drastic change in the international trade with a transition from dollar dominance to yuan dominance.

[16] C. Coquidé, L. [START_REF] Blondel | Fasdt unfolding of communities in large networks[END_REF] with Dugué's algorithm [START_REF] Dugué | Directed Louvain : maximizing modularity in directed networks[END_REF]. Each cluster is labelled using the ISO2 code of the country with the best import and export trade probabilities Pc and P * c . For all the considered years, the leaders of the clusters are the USA (blue), China (red), and Germany (green). The cluster Others (gold) gathers other small size clusters. S4. Same as Fig. 2 but keeping always USA, Canada, UK, Australia, New Zealand trading in USD and China, Brazil, Russia, India, South Africa in CNY. Final fraction f f of countries preferring to trade in USD as a function of the initial fraction fi for the years 2010 (left panel) and 2019 (right panel). Left and right panels: for any initial fraction fi, the Monte Carlo procedure converges toward one of two final fractions f f . These two final fractions are f f 1 = 0.22 and f f 2 = 0.62 in 2010 (left panel) and f f 1 = 0.18 and f f 2 = 0.45 in 2019 (right panel). The color of a point (fi,f f ) indicates the portion ρ f f (fi) of the Nr = 10 4 initial configurations, with a corresponding initial fraction fi, which attain the final state with the corresponding final fraction f f . The color ranges from black for ρ f f (fi) = 0 (all the countries preferring to trade in CNY whether than in USD) to bright yellow for ρ f f (fi) = 1 (all the countries preferring to trade in USD whether than in CNY). Middle panel: Portion ρ f f (fi) of the Nr = 10 4 initial configurations, the fraction fi of which initially prefers to trade in USD, which attains the final state with the final fraction f f . The red (blue) curve and the up (down) triangles correspond to the lowest (highest) value of the two final fractions f f . The full (empty) symbols correspond to the year 2019 (2010). S6. Same as Fig. 2 but replacing the import trade probabilities Pc and the export trade probabilities P * c in (1) by the PageRank and the CheiRank probabilities obtained from the Google matrix of WTN. Final fraction f f of countries preferring to trade in USD as a function of the initial fraction fi for the years 2010 (left panel) and 2019 (right panel). Left and right panels: for any initial fraction fi, the Monte Carlo procedure converges toward one of two final fractions f f . These two final fractions are f f 1 = 0.16 and f f 2 = 0.67 in 2010 (left panel) and f f 1 = 0.14 and f f 2 = 0.44 in 2019 (right panel). The color of a point (fi,f f ) indicates the portion ρ f f (fi) of the Nr = 10 4 initial configurations, with a corresponding initial fraction fi, which attain the final state with the corresponding final fraction f f . The color ranges from black for ρ f f (fi) = 0 (all the countries preferring to trade in CNY whether than in USD) to bright yellow for ρ f f (fi) = 1 (all the countries preferring to trade in USD whether than in CNY). Middle panel: Portion ρ f f (fi) of the Nr = 10 4 initial configurations, the fraction fi of which initially prefers to trade in USD, which attains the final state with the final fraction f f . The red (blue) curve and the up (down) triangles correspond to the lowest (highest) value of the two final fractions f f . The full (empty) symbols correspond to the year 2019 (2010).

ρ(f i ) f f 1 ρ (f i ) f f 2 FIG.
f i f f ρ (f i ) f f 1 ρ (f i ) f f 2 FIG.

FIG. 4 .

 4 FIG. 4.World distribution of countries belonging to the USD group (blue, hard preference to trade in USD), the CNY group (red, hard preference to trade in CNY) and the swing group (green, the TCP can change between USD and CNY depending on the initial conditions). The world maps are shown for the years 2010 (left panel) and 2019 (right panel).

FIG. 5 .

 5 FIG.5. Time evolution of the sizes of the trade currency preference groups. The 3 bands correspond to the USD group (blue), the CNY group (red), and the swing group (green). The width of a band corresponds to the size of the corresponding group expressed as: (left panel) the ratio between the countries belonging to the group over and the total number of countries, (right panel) the ratio between the total trade volume exchanged by the countries of the group and the total volume exchanged by all the world countries. The trade volumes are expressed in USD of the concerned year.

  FIG. S3. World map distribution of countries belonging to the USD group (blue, hard preference to trade in USD), the CNY group (red, hard preference to trade in CNY) and the swing group (green, the TCP can change between USD and CNY depending on the initial conditions). The world maps are shown for the years2010, 2012, 2014, 2016, 2018, 2019 and 2020.

  FIG.S5. Same as Fig.4but keeping always USA, Canada, UK, Australia, New Zealand trading in USD and China, Brazil, Russia, India, South Africa in CNY. World map distribution of countries belonging to the USD group (blue, hard preference to trade in USD), the CNY group (red, hard preference to trade in CNY) and the swing group (green, the TCP can change between USD and CNY depending on the initial conditions). The world maps are shown for the years 2010 (left panel) and 2019 (right panel).

  FIG. S7. Same as Fig. 4 but replacing the import trade probabilities Pc and the export trade probabilities P * c in (1) by the PageRank and the CheiRank probabilities obtained from the Google matrix of WTN. World map distribution of countries belonging to the USD group (blue, hard preference to trade in USD), the CNY group (red, hard preference to trade in CNY) and the swing group (green, the TCP can change between USD and CNY depending on the initial conditions). The world maps are shown for the years 2010 (left panel) and 2019 (right panel).

  Ermann, J. Lages, and D. Shepelyansky, Influence of petroleum and gas trade on EU economies from the reduced Google matrix analysis of UN COM-TRADE data, Eur. Phys. J. B 92, 71 (2019).

					2010			2019 2012
	US	CN	DE DE	Others Others		US	CN	DE
					2014			2016
	US	CN	DE			US	CN	DE
					2018			2019
	US	CN	DE DE			US	CN	DE
							2020
				US	CN	DE
	FIG. S2. Geographical distribution of clusters in the WTN for the years 2010, 2012, 2014, 2016, 2018, 2019 and 2020 obtained
	using the Louvain modularity method

TABLE S1 .

 S1 List of the 32 countries belonging to the USD group in 2010. The trade currency preference of these countries is always USD at the end of the relaxation process. The countries are sorted by descending value of max (Pc, P * c ), i.e. the maximum value between the relative import volume Pc and the relative export volume P * c , and, in case of tie, by descending value of P * c . The red (green) colored countries switch to the CNY (Swing) group in 2019. The countries are represented by their ISO2 codes.Countries of the USD group in 2010

TABLE S2 .

 S2 List of the 66 countries belonging to the CNY group in 2010. The trade currency preference of these countries is always CNY at the end of the relaxation process. The countries are sorted by descending value of max (Pc, P * c ), i.e. the maximum value between the relative import volume Pc and the relative export volume P * c , and, in case of tie, by descending value of P * c . The blue colored countries switch to the USD group in 2019. The countries are represented by their ISO2 codes. Countries of the CNY group in 2010 TABLE S3. List of the 96 countries belonging to the Swing group in 2010. The trade currency preference of these countries is for all of them either CNY or USD at the end of the relaxation process. The countries are sorted by descending value of max (Pc, P * c ), i.e. the maximum value between the relative import volume Pc and the relative export volume P * c , and, in case of tie, by descending value of P * c . The red (blue) colored countries switch to the CNY (USD) group in 2019. The countries are represented by their ISO2 codes. Countries of the Swing group in 2010 TABLE S4. List of the 28 countries belonging to the USD group in 2019. The trade currency preference of these countries is always USD at the end of the relaxation process. The countries are sorted by descending value of max (Pc, P * c ), i.e. the maximum value between the relative import volume Pc and the relative export volume P * c , and, in case of tie, by descending value of P * c . The red (green) colored countries belonged to the CNY (Swing) group in 2010. The countries are represented by their ISO2 codes. Countries of the USD group in 2019 TABLE S5. List of the 109 countries belonging to the CNY group in 2019. The trade currency preference of these countries is always CNY at the end of the relaxation process. The countries are sorted by descending value of max (Pc, P * c ), i.e. the maximum value between the relative import volume Pc and the relative export volume P * c , and, in case of tie, by descending value of P * c . The blue (green) colored countries belonged to the USD (Swing) group in 2010. The countries are represented by their ISO2 codes. Countries of the CNY group in 2019 TABLE S6. List of the 57 countries belonging to the Swing group in 2019. The trade currency preference of these countries is for all of them either CNY or USD at the end of the relaxation process. The countries are sorted by descending value of max (Pc, P * c ), i.e. the maximum value between the relative import volume Pc and the relative export volume P * c , and, in case of tie, by descending value of P * c . The blue colored countries belonged to the USD group in 2010. The countries are represented by their ISO2 codes. Countries of the Swing group in 2019

ACKNOWLEDGMENTS

We thank L. Ermann for his help to collect data from the UN Comtrade database. We thank the UN Statistics Division to grant us a friendly access to the UN Comtrade database. This research has been partially supported through the grant NANOX N • ANR-17-EURE-0009 (project MTDINA) in the frame of the Programme des Investissements d'Avenir, France. This research has been also supported by the Programme Investissements d'Avenir ANR-15-IDEX-0003.

Supplementary Information