Spectral signatures of structural change in financial networks
Valentina Macchiati Emiliano Marchese Piero Mazzarisi Diego Garlaschelli Tiziano Squartini





Scuola Normale Superiore, P.zza dei Cavalieri 7, 56126 Pisa (Italy)

IMT School for Advanced Studies, P.zza San Francesco 19, 55100 Lucca (Italy)

Università degli Studi di Siena, P.zza S.Francesco 7-8, 53100 Siena (Italy)

IMT School for Advanced Studies, P.zza San Francesco 19, 55100 Lucca (Italy)
Lorentz Institute for Theoretical Physics, University of Leiden, Niels Bohrweg, 2, Leiden, NL-2333 CA, The Netherlands
INdAM-GNAMPA Istituto Nazionale di Alta Matematica `Francesco Severi', P.le Aldo Moro 5, 00185 Rome (Italy)

tiziano.squartini@imtlucca.it
IMT School for Advanced Studies, P.zza San Francesco 19, 55100 Lucca (Italy)
Scuola Normale Superiore, P.zza dei Cavalieri 7, 56126 Pisa (Italy)
INdAM-GNAMPA Istituto Nazionale di Alta Matematica `Francesco Severi', P.le Aldo Moro 5, 00185 Rome (Italy)



§ ABSTRACT 
The level of systemic risk in economic and financial systems is strongly determined by the structure of the underlying networks of interdependent entities that can propagate shocks and stresses. Since changes in network structure imply changes in risk levels, it is important to identify structural transitions potentially leading to system-wide crises. Methods have been proposed to assess whether a real-world network is in a (quasi-)stationary state by checking the consistency of its structural evolution with appropriate maximum-entropy ensembles of graphs. While previous analyses of this kind have focused on dyadic and triadic motifs, hence disregarding higher-order structures, here we consider closed walks of any length. Specifically, we study the ensemble properties of the spectral radius of random graph models calibrated on real-world evolving networks. Our approach is shown to work remarkably well for directed networks, both binary and weighted. As illustrative examples, we consider the Electronic Market for Interbank Deposit (e-MID), the Dutch Interbank Network (DIN) and the International Trade Network (ITN) in their evolution across the 2008 crisis. By monitoring the deviation of the spectral radius from its ensemble expectation, we find that the ITN remains in a (quasi-)equilibrium state throughout the period considered, while both the DIN and e-MID exhibit a clear out-of-equilibrium behaviour. The spectral deviation therefore captures ongoing topological changes, extending over all length scales, to provide a compact proxy of the resilience of economic and financial networks.


89.75.Fb; 89.65.-s; 02.50.Tt

Spectral signatures of structural change in financial networks
    Tiziano Squartini
    September 9, 2024
==============================================================





§ INTRODUCTION


As witnessed by two major recent crises (i.e. the global financial one in 2008 and the Covid-19 pandemic in 2020), having a clear understanding of the intricate structure of economic and financial systems - be they interbank <cit.>, interfirm <cit.> or trade networks <cit.> - is crucial, especially under stress conditions. The interconnectedness of economic and financial agents is, in fact, known to play a major role both during the phase of distress accumulation and after a crisis outbreak in sustaining and reinforcing shock propagation <cit.>. Back in 2008, banks sought to minimise individual risk by diversifying their portfolios: the simultaneous character of such diversification, however, led to an unexpected level of mutual dependency whose net consequence was that of amplifying the effects of individual defaults <cit.>.

A particularly relevant question addresses the (quasi-)stationarity of the temporal evolution of a given, real-world, economic or financial network, i.e. does the system undergo smooth, structural changes controlled by few driving parameters? Should this be the case, the behaviour of the network under analysis would be predictable solely in terms of the dynamics of those parameters; otherwise, the lack of stationarity may lead to abrupt - hence, uncontrollable - regime shifts.

The problem of the (non) stationarity of real-world, economic and financial networks has been addressed by studying whether they can be considered typical members of an evolving, (quasi-)equilibrium ensemble of graphs with given properties <cit.>: while such properties are treated as constraints - hence, assumed to be the `independent variables' undergoing an autonomous evolution - the other network properties are treated as `dependent vareiables' - hence, assumed to vary only as a consequence of the former ones. Broadly speaking, three different situations can occur:



  ∙ The observed network properties are systematically found to agree with what is expected from the evolution of the enforced constraints. In this case, one can conclude that the real-world network is (quasi-)stationary - and its evolution is driven by the dynamics of the constraints;

  ∙ The observed network properties slightly deviate from equilibrium expectations, but the deviating patterns remain coherent. In this case, the network can still be considered (quasi-)stationary - even if its evolution cannot be claimed to be completely driven by the chosen constraints (very likely, with the addition of other appropriate constraints, one would go back to the first situation);

  ∙ The observed network properties significantly deviate from the (quasi-)equilibrium expectations, showing different deviating patterns at different times. In this case, the network can be considered non-stationary.


Analyses of this kind have indeed led to individuate early-warning signals of upcoming, critical events, although the indicators considered so far have just involved dyadic and triadic `debt loops' with different levels of reciprocity <cit.>. The present paper aims to extend the study of early-warnings' emergence by considering closed walks of any length at once. Such a request can be handled by exploiting the theorem stating that a_ij^(n), i.e. the generic entry of the n-th power of the adjacency matrix 𝐀, counts the total number of closed walks of length n connecting node i with node j: all `debt loops' can be, then, compactly accounted for by carrying out a double sum, over n and over the diagonal entries of 𝐀.

From a computational perspective, such a calculation can be greatly sped up by proxying the trace of the adjacency matrix with its principal eigenvalue λ_1, which, then, becomes the only relevant statistics whose z-score needs to be explicitly calculated. Such an appealing simplification, however, comes at a price: the expressions of ⟨λ_1⟩ and Var[λ_1], i.e. of the expected value of λ_1 and its variance, are explicitly known in few cases only, i.e. i) when the random network model is the binary, undirected version of the Erdös-Rényi (ER) model <cit.>; ii) when the random network model is the Chung-Lu (CL) model, either in its binary, undirected version <cit.> or in its binary, directed version <cit.>; iii) if the edges are treated as i.n.i.d. (independent, non-identically distributed) random variables, each one obeying a different Poisson distribution <cit.>; iv) if the considered graphs are infinitely large, locally tree-like and directed <cit.>.

Let us remark that the existing estimations obtained under hypotheses are rarely satisfied by empirical configurations. For instance, the presence of cycles contradicts the assumption of observing locally tree-like structures, and the heterogeneity of the (in- and out-) degree distributions severely limits the applicability of the CL model. On a more general ground, the vast majority of the approaches above requires the knowledge of the (in- and out-) degree sequences, i.e. of a kind of information that data confidentiality issues make often (if not always) unavailable; moreover, none provides estimations of a network spectral properties taking its weighted marginals (i.e. in-strengths and out-strengths) as the sole input.

Motivated by the evidence that general results about the statistical properties of a network principal eigenvalue are currently missing, we propose an approach to their study that is applicable under any random network model. The generality of our approach comes at a price: our results rest upon the validity of several approximations that need to be explicitly verified whenever a particular configuration is studied. Still, although our assumptions may appear quite drastic, our approach works remarkably well for directed networks, be they binary (BDN) or weighted (WDN).

A BDN is described by an adjacency matrix 𝐀 whose generic entry satisfies the relationships a_ij=1 if a link points from node i towards node j and a_ij=0 otherwise. Moreover, a_ij will, in general, differ from a_ji. A WDN is described by an adjacency matrix 𝐖 whose generic entry satisfies the relationships w_ij>0 if a weighted link points from node i towards node j and w_ij=0 otherwise. Moreover, w_ij will, in general, differ from w_ji.



§ DETECTING STRUCTURAL CHANGES


Structural changes can be spotted by comparing the empirical abundance of a quantity of interest with the corresponding expected value, calculated under a properly defined benchmark model[Hereby, the expressions `random network model', `benchmark model' and `null model' will be used interchangeably.]. To this aim, a very useful indicator is represented by the z-score


    z[X]=X-⟨ X⟩/σ[X]

where X is the empirical abundance of the quantity X, ⟨ X⟩ is its expected occurrence under the chosen null model and σ[X]=√(⟨ X^2⟩-⟨ X⟩^2) is the standard deviation of X under the same null model. In words, z[X] quantifies the number of standard deviations by which the empirical abundance of X differs from the expected one after checking for the Gaussianity of X under the null model - often ensured by the fact that X is the sum of several random variables - a result |X|≤2 (|X|≤3) indicates that the empirical abundance of X is compatible with the expected one, at the 5% (1%) level of statistical significance; on the other hand, a value |X|>2 (|X|>3) indicates that the empirical abundance of X is not compatible with the expected one, at the same significance level. In the latter case, a value z[X]>0 (z[X]<0) indicates the tendency of the pattern to be over-represented (under-represented) in the data with respect to the chosen benchmark.



 §.§ Dyadic signature of structural changes


Moving from the observation that


    ∑_j=1^Na_ija_jk=[𝐀^2]_ik

we will pose


    X   =∑_i=1^N∑_j=1^Na_ija_ji=∑_i=1^N[𝐀^2]_ii=Tr[𝐀^2],

noticing that the total number of links having a partner pointing in the opposite direction coincides with the trace of the second power of the adjacency matrix. The position above leads to


    ⟨ X⟩=∑_i=1^N∑_j=1^Np_ijp_ji=∑_i=1^N[𝐏^2]_ii=Tr[𝐏^2]

and to


    σ[X]   =√(Var[∑_i=1^N∑_j=1^Na_ija_ji])
       =√(Var[2·∑_i=1^N∑_j(>i)a_ija_ji])
       =√(4·∑_i=1^N∑_j(>i)Var[a_ija_ji])
       =2·√(∑_i=1^N∑_j(>i)p_ijp_ji(1-p_ijp_ji))

where 𝐏≡{p_ij}_i,j=1^N is the matrix of probability coefficients induced by the chosen null model, and the third passage follows from the evidence that the dyads induce independent random variables (see also Appendix AppAA).



 §.§ Triadic signature of structural changes


Analogously to the dyadic case, let us move from the observation that


    ∑_j=1^N∑_k=1^Na_ija_jka_kl=[𝐀^3]_il

and pose


    X   =∑_i=1^N∑_j=1^N∑_k=1^Na_ija_jka_ki=∑_i=1^N[𝐀^3]_ii=Tr[𝐀^3],

noticing that the total number of triangles is proportional to the trace of the third power of the adjacency matrix. The position above leads to


    ⟨ X⟩=∑_i=1^N∑_j=1^N∑_k=1^Np_ijp_jkp_ki=∑_i=1^N[𝐏^3]_ii=Tr[𝐏^3]

and to


    σ[X]   =√(Var[∑_i=1^N∑_j=1^N∑_k=1^Na_ija_jka_ki])
       =√(Var[3·∑_i=1^N∑_j(>i)∑_k(>j)(a_ija_jka_ki+a_ika_kja_ji)])
       =√(9·Var[∑_i=1^N∑_j(>i)∑_k(>j)(a_ija_jka_ki+a_ika_kja_ji)])
       =3·√(Var[∑_i=1^N∑_j(>i)∑_k(>j)(a_ija_jka_ki+a_ika_kja_ji)])

where 𝐏≡{p_ij}_i,j=1^N is the matrix of probability coefficients induced by the chosen null model. Since triads do not induce independent random variables, the explicit expression of σ[X] is derived in Appendix AppBB. Let us notice that, in case the considered networks are sparse, one can simplify the expression above upon posing


    σ[X]   ≃3·√(∑_i=1^N∑_j(>i)∑_k(>j)Var[a_ija_jka_ki+a_ika_kja_ji])

with


    Var[a_ija_jka_ki+a_ika_kja_ji]   ≃Var[a_ija_jka_ki]+Var[a_ika_kja_ji]

and


    Var[a_ija_jka_ki]   =p_ijp_jkp_ki(1-p_ijp_jkp_ki),
    Var[a_ika_kja_ji]   =p_ikp_kjp_ji(1-p_ikp_kjp_ji).




 §.§ Spectral signature of structural changes


Let us now enlarge the set of patterns to be considered for detecting structural changes by accounting for closed walks of any length.



  §.§.§ The trace of the matrix exponential


Let us start by considering the N× N adjacency matrix 𝐀 of a BDN, with a_ii=0, ∀ i: the following relationship


    𝐈+𝐀+𝐀^2/2!+𝐀^3/3!+…+𝐀^n/n!+…=∑_k=0^∞𝐀^k/k!≡ e^𝐀,

where 𝐀^0≡𝐈, defines the exponential of 𝐀 <cit.>. Let us, now, calculate the trace of such a matrix exponential: since it is invariant under diagonalisation, one obtains that


    Tr[e^𝐀]=∑_k=0^∞Tr[𝐀^k]/k!=∑_k=0^∞Tr[Λ^k]/k!=Tr[e^Λ]

where Λ is the matrix obtained upon diagonalising 𝐀 (see also Appendix AppCC). As the number of walks of length k starting from and ending at the same vertex can be counted by computing the trace of the k-th power of the adjacency matrix, i.e. Tr[𝐀^k]=∑_i=1^N[𝐀^k]_ii, eq. (<ref>) relates the number of walks of any length characterising a binary network 𝐀 with its spectral properties. Such a quantity, named Estrada index, represents a graph invariant quantifying the communicability of a given network, i.e. the `participation' of each node to the walks present in the network itself <cit.>.


Analogously, given the N× N adjacency matrix 𝐖 of a WDN with w_ii=0, ∀ i, the relationships


    𝐈+𝐖+𝐖^2/2!+𝐖^3/3!+…+𝐖^n/n!+…=∑_k=0^∞𝐖^k/k!≡ e^𝐖

and


    Tr[e^𝐖]=∑_k=0^∞Tr[𝐖^k]/k!=∑_k=0^∞Tr[Ω^k]/k!=Tr[e^Ω],

where Ω is the matrix obtained upon diagonalising 𝐖, hold true. As a result, concerning the number of walks of length k starting and ending at the same vertex can be extended to weighted networks, eq. (<ref>) generalises the Estrada index to weighted configurations.


Let us explicitly notice that



  ∙ the absence of self-loops, i.e. Tr[𝐀]=Tr[𝐖]=0, implies that, whenever present, complex eigenvalues must appear in conjugate pairs;

  ∙ eq. (<ref>) implies that Tr[e^𝐀]≥0, i.e. that the trace of the exponential of 𝐀 is a real, non-negative number. Analogously, eq. (<ref>) implies that Tr[e^𝐖]≥0, i.e. that the trace of the exponential of 𝐖 is a real, non-negative number;

  ∙ When computing the number of closed walks of a certain length, edges must be counted repeatedly. For example, the closed walks of length 4 in a binary, directed network are i) the proper cycles like i→ j→ k→ l→ i; ii) the pairs of dyads like i→ j→ k→ j→ i; iii) the single dyads like i→ j→ i. Equations (<ref>) and (<ref>) compactly account for all of them.


The third observation has relevant implications for economic and financial applications: when studying the propagation of a shock, in fact, it is extremely important to account for all possible patterns along which distress can propagate, including the ones leading to multiple reverberations among the same nodes <cit.>. As (combinations of) cycles are supposed to lower the resilience of financial networks by amplifying external shocks <cit.>, eqs. (<ref>) and (<ref>) suggest the trace of the exponential matrix to represent a compact proxy of the stability of the network itself.



  §.§.§ Expected value of the trace of the matrix exponential


Let us now move to analyse the expected value of the quantity Tr[e^𝐀]=Tr[e^Λ], under a properly-defined benchmark model. We will suppose the latter one to be described by an N× N matrix 𝐏 whose generic entry p_ij, with i≠ j, indicates the probability that nodes i and j are connected via a directed link. Following the same steps as above, we find


    Tr[e^𝐏]=∑_k=0^∞Tr[𝐏^k]/k!=∑_k=0^∞Tr[Π^k]/k!=Tr[e^Π]

where Π is the matrix obtained upon diagonalising 𝐏. 

Let us now inspect the relationship between eq. (<ref>) and eq. (<ref>). Since we are considering binary, adjacency matrices, the matrix 𝐏 satisfies the relationship ⟨𝐀⟩=𝐏, a compact notation stating for ⟨ a_ij⟩=p_ij, ∀ i≠ j. To extend this result to higher powers of the adjacency matrix, an explicit expression for the quantity ⟨𝐀^n⟩=f(𝐏), ∀ n is needed. Here, we adopt the recipe defining the so-called delta method <cit.> and prescribing to identify f(𝐏) with 𝐏^n. According to it, the expected value of the number of closed walks of any length satisfies the chain of inequalities



    ⟨Tr[e^𝐀]⟩=∑_k=0^∞⟨Tr[𝐀^k]⟩/k!=∑_k=0^∞Tr[⟨𝐀^k⟩]/k!≥∑_k=0^∞Tr[⟨𝐀⟩^k]/k!=∑_k=0^∞Tr[𝐏^k]/k!=Tr[e^𝐏];


a relationship leading to ⟨Tr[e^Λ]⟩≥Tr[e^Π].

The inequality can be understood upon considering a reciprocated dyad and noticing that relationships like ⟨[𝐀^4]_ii⟩=⟨ a_ija_jia_ija_ji⟩=⟨ a_ija_ji⟩=⟨ a_ij⟩⟨ a_ji⟩=p_ijp_ji≥ p_ij^2p_ji^2=[𝐏^4]_ii hold true; in other words, estimating the number of closed walks of a certain length via the delta method implies overweighing the edges constituting them, whence the mismatch between the correct and the approximated expression. Such a mismatch is absent in case no link is reciprocated: given a square loop, in fact, ⟨[𝐀^4]_ii⟩=⟨ a_ija_jka_kla_li⟩=⟨ a_ij⟩⟨ a_jk⟩⟨ a_kl⟩⟨ a_li⟩=p_ijp_jkp_klp_li=[𝐏^4]_ii. In other words, the larger the number of reciprocal links[Let us remind that L^↔=∑_i=1^N∑_j(≠ i)a_ija_ji.], the less accurate the approximation provided by the delta method. Hereby, we will assume that the symbol ≳ can replace the symbol ≥.


Analogously, upon posing ⟨𝐖⟩=𝐐, the expected value of the quantity Tr[e^𝐖]=Tr[e^Ω] can be approximated as follows



    ⟨Tr[e^𝐖]⟩=∑_k=0^∞⟨Tr[𝐖^k]⟩/k!=∑_k=0^∞Tr[⟨𝐖^k⟩]/k!≥∑_k=0^∞Tr[⟨𝐖⟩^k]/k!=∑_k=0^∞Tr[𝐐^k]/k!=Tr[e^𝐐],


a relationship leading to ⟨Tr[e^Ω]⟩≥Tr[e^Ψ], where Ψ is the matrix obtained upon diagonalising 𝐐.

The inequality can be understood upon considering a weighted, reciprocated dyad and noticing that relationships like ⟨[𝐖^4]_ii⟩=⟨ w_ijw_jiw_ijw_ji⟩=⟨ (w_ijw_ji)^2⟩=⟨ w_ijw_ji⟩^2+Var[w_ijw_ji]=⟨ w_ij⟩^2⟨ w_ji⟩^2+Var[w_ijw_ji]≥⟨ w_ij⟩^2⟨ w_ji⟩^2=[𝐐]^4_ii hold true; as in the binary case, estimating the total weight of closed walks of a certain length via the delta method implies overweighing the edges constituting them. Such a mismatch is absent if no link is reciprocated, as evident upon considering a weighted, square loop. Hereby, we will assume that the symbol ≳ can replace the symbol ≥.



  §.§.§ Expected value of the spectral radius


Let us now recall the statement of the generalised Perron-Frobenius (GPF) theorem <cit.>.


GPF Theorem. Whenever non-negative, irreducible matrices are considered, a unique, real, positive eigenvalue exists whose modulus is maximum and (only) the corresponding left and right eigenvectors have positive components.


Requiring irreducibility, sometimes stated as regularity, implies requiring the existence of a natural number n such that [𝐀^n]_ij>0, ∀ i,j. In other words, when directed networks are considered, requiring irreducibility is equivalent to requiring strongly connectedness. In case such a requirement is not satisfied, the Perron-Frobenius theorem must be weakened as follows.


WPF Theorem. Whenever non-negative matrices are considered, a real, non-negative eigenvalue exists whose modulus is maximum and with associated, non-negative left and right eigenvectors.


The eigenvalue mentioned in any variant of the Perron-Frobenius theorem will be referred to as the principal eigenvalue or spectral radius. The relationship between the matrices 𝐀 and Λ encoded into eq. (<ref>) can be further simplified upon noticing that, in case the spectral radius exists, is unique[A reducible square matrix 𝐌 can be written in a block triangular form <cit.>, each matrix 𝐁_ii on the diagonal being either irreducible or zero. As the spectrum of such a matrix is the union of the spectra of the 𝐁_iis, the GPF Theorem can be applied to each 𝐁_ii: the Perron–Frobenius eigenvalue of 𝐌 is, thus, the largest of the Perron–Frobenius eigenvalues of the 𝐁_iis, hence coinciding with the one of the maximal strongly-connected component of the network under study.] and the spectral gap is (much) larger than zero[Although the condition λ_1-λ_2≫0 can be relaxed, the formulas provided in the present paper hold for this case.], the sum Tr[e^Λ]=∑_i=1^Ne^λ_i is exponentially dominated by the addendum e^λ_1, an observation allowing us to write


    Tr[e^𝐀]≳ e^λ_1;

analogously, 


    Tr[e^𝐖]≳ e^ω_1.


Let us now inspect the relationships between eqs. (<ref>) and (<ref>) and between eqs. (<ref>) and (<ref>). Putting everything together, we obtain



    ⟨ e^λ_1⟩≲⟨Tr[e^𝐀]⟩=∑_k=0^∞⟨Tr[𝐀^k]⟩/k!   =∑_k=0^∞Tr[⟨𝐀^k⟩]/k!≳∑_k=0^∞Tr[𝐏^k]/k!=Tr[e^𝐏]≳ e^π_1,
    ⟨ e^ω_1⟩≲⟨Tr[e^𝐖]⟩=∑_k=0^∞⟨Tr[𝐖^k]⟩/k!   =∑_k=0^∞Tr[⟨𝐖^k⟩]/k!≳∑_k=0^∞Tr[𝐐^k]/k!=Tr[e^𝐐]≳ e^ϕ_1.



The two chains of (in-)equalities above motivate us to explore the possibility of deriving an (approximated) expression for the expected value of the spectral radius. According to the delta method, the expected value of a function, f, of a random variable, x, can be computed by Taylor-expanding f(x) around ⟨ x⟩=μ, taking the expected value of the resulting expression and retaining only the lowest order of the expansion. Such a prescription allows us to write ⟨ e^λ_1⟩≃ e^⟨λ_1⟩ and ⟨ e^ω_1⟩≃ e^⟨ω_1⟩, two positions further leading to the results


    ⟨λ_1⟩≃π_1

and


    ⟨ω_1⟩≃ϕ_1.


Equations (<ref>) and (<ref>) are the main result of our paper, as they establish a (fundamental, although approximated) relationship between the empirical value of the spectral radius of a directed network, be it binary or weighted, and its expected counterpart: in words, the delta method suggests us to identify the latter with the spectral radius of the matrix defining the chosen random network model. Since the calculation of the expected number, or of the expected weight, of walks boils down to calculate the spectral radius of a single matrix, i.e. 𝐏 or 𝐐, eqs. (<ref>) and (<ref>) have deep implications from a purely computational point of view as well: in fact, they prevent the network ensemble induced by 𝐏 or 𝐐 from being explicitly sampled.



  §.§.§ Variance of the spectral radius


Now, let us focus on the variance of the spectral radius calculation. To this aim, we will move from the known expressions of ⟨λ_1⟩, treating them as subject to statistical variability. For instance, let us recall that 


    π_1≃⟨k||k⟩/2L=∑_i=1^Nk_i^2/2L

for binary, undirected networks under the Chung-Lu model, according to which p_ij=k_ik_j/2L, ∀ i,j; upon considering that all quantities defining such an expression are random variables themselves, one is led to write


    Var[λ_1]=Var[∑_i=1^Nk_i^2/2L]

and evaluate such an expression either analytically or numerically. In what follows, we will numerically evaluate the spectral radius variance of our random network models.



  §.§.§ Statistical significance of the spectral radius


Let us now define the quantity to be inspected for spotting the presence of a spectral signature of structural changes: it reads


    z[λ_1]=λ_1-⟨λ_1⟩/σ[λ_1]≃λ_1-π_1/σ[λ_1]

and is nothing but the z-score of the spectral radius λ_1. As already stressed, the statistical meaning of such a quantity is guaranteed by the Gaussianity of the quantity whose z-score is to be calculated. Such a property of the spectral radius is guaranteed by the analytical results obtained in <cit.> and by the numerical checks carried out in Appendix AppDD and depicted in fig. <ref>.



§ RANDOM NETWORK MODELS


Let us now discuss a set of null models to be employed for the subsequent steps of our analysis. To this aim, we will consider some members of the family of Exponential Random Graph Models (ERGMs), i.e. the entropy-based benchmarks that preserve a given set of constraints, otherwise being maximally random. More specifically, we follow the approach introduced in <cit.> and further developed in <cit.>, which prescribes to carry out a constrained maximisation of Shannon entropy


    S=-∑_𝐆P(𝐆)ln P(𝐆),

the sum running over the ensemble 𝔾 of N× N directed networks, be they binary (in which case 𝐆≡𝐀) or weighted (in which case 𝐆≡𝐖).



 §.§ Erdös-Rényi Model


The Erdös-Rényi Model <cit.> is induced by the Hamiltonian 


    H(𝐀)=α L(𝐀),

where L(𝐀)=∑_i=1^N∑_j(≠ i)a_ij represents the total number of directed edges, and α is the Lagrange multiplier associated with such a global constraint. The probability of the generic configuration 𝐀 reads


    P_ER(𝐀)=p^L(𝐀)(1-p)^N(N-1)-L(𝐀)

where p=e^-α/(1+e^-α) is the probability that a link points from node i towards node j.

In order to tune the unknown parameter defining the Erdös-Rényi Model to ensure that ⟨ L⟩_ER=L(𝐀^*), we maximise the likelihood function ℒ_ER=ln P_ER(𝐀^*) with respect to it. Such a recipe leads us to find


    p=L(𝐀^*)/N(N-1),  ∀ i≠ j

with obvious meaning of the symbols.



  §.§.§ Expected value of the spectral radius


Although eq. (<ref>) provides a general recipe for estimating the expected value of the spectral radius of any random network model, a more explicit expression can be derived for the Erdös-Rényi Model. Specifically, let us consider the following equation


    ∑_k=0^∞Tr[𝐏^k]/k!=N+∑_k=2^∞(Np)^k/k!

where 𝐏≡𝐏_ER={p_ij}_i,j=1^N, p_ij≡ p, ∀ i≠ j and each addendum encodes the information about the order of magnitude of the specific contribution to the total number of cycles - to see this explicitly, let us consider that Tr[𝐀^2]=∑_i=1^N[𝐀^2]_ii=∑_i=1^N∑_j(≠ i)a_ija_ji whose expected value reads ⟨Tr[𝐀^2]⟩=∑_i=1^N∑_j(≠ i)⟨ a_ija_ji⟩=∑_i=1^N∑_j(≠ i)p^2≃(Np)^2 and analogously for the higher orders of the expansion. As adding and subtracting 1 and Np leads to


    ∑_k=0^∞Tr[𝐏^k]/k!   =∑_k=0^∞(Np)^k/k!+N(1-p)-1
       =e^Np+N(1-p)-1≳ e^Np≃ e^⟨ k⟩,

eq. (<ref>) can be employed to derive the chain of relationships


    π_1≃ Np≃⟨ k⟩,

stating that the spectral radius, π_1, of the N× N matrix of i.i.d. Bernoulli random variables 𝐏≡𝐏_ER can be accurately approximated by their sum along any row or any column; in network terms, this can be rephrased by saying that the expected value of the spectral radius under the Erdös-Rényi Model coincides with the expected value of the degree of each node.


A second way of identifying π_1 rests upon the following relationship:


    𝐏·1=(N-1)p·1=⟨ k⟩·1;

since 𝐏 obeys the GPF Theorem, the equation above allows us to identify the value of its spectral radius quite straightforwardly by posing


    π_1=(N-1)p=⟨ k⟩≡λ_1^ER.


Such a result is consistent with the one stating that the spectral radius of the deterministic matrix a_ii≡ν, ∀ i=j and a_ij≡μ, ∀ i≠ j is equal to λ_1=(N-1)μ+ν.


A third way of identifying the expected value of λ_1 rests upon the results from <cit.>, i.e.


    λ_1=∑_i=1^N∑_ja_ij/N+σ^2/μ,

where a_ii≡ν, ∀ i=j, ⟨ a_ij⟩=μ and Var[a_ij]=σ^2, ∀ i≠ j. Since, in our case, ν=0, μ=p and σ^2=p(1-p), ∀ i≠ j, such an expression leads to


    λ_1=∑_i=1^N∑_jp/N+σ^2/μ=(N-1)p+(1-p).




  §.§.§ Variance of the spectral radius


Equation (<ref>) offers a straightforward way to calculate the variance of the spectral radius. It is, in fact, enough to evaluate the expression


    Var[λ_1]=∑_i=1^N∑_j(≠ i)Var[a_ij]/N^2   ≃ p(1-p)≡Var[λ_1^ER]

with the symbol ≃ replacing the more correct expression lim_N→∞ N(N-1)p(1-p)/N^2=p(1-p), indicating that Var[λ_1] tends to p(1-p) in the (asymptotic) regime N→∞.



 §.§ Binary Configuration Model


The Binary Configuration Model <cit.> is induced by the Hamiltonian 


    H(𝐀)=∑_i=1^N[α_ik_i(𝐀)+β_ih_i(𝐀)]

where k_i(𝐀)=∑_j(≠ i)a_ij represents the out-degree of node i, i.e. the number of nodes pointed by it and h_i(𝐀)=∑_j(≠ i)a_ji represents the in-degree of node i, i.e. the number of nodes it is pointed by; the vectors {α_i}_i=1^N and {β_i}_i=1^N represent the Lagrange multipliers associated with those above, local constraints. The probability of the generic configuration 𝐀 reads


    P_BCM(𝐀)=∏_i=1^N∏_j(≠ i)p_ij^a_ij(1-p_ij)^1-a_ij

where p_ij=e^-α_i-β_j/(1+e^-α_i-β_j) is the probability that a link points from node i towards node j.

To tune the unknown parameters defining the Binary Configuration Model to ensure that ⟨ k_i⟩_BCM=k_i(𝐀^*), ∀ i and ⟨ h_i⟩_BCM=h_i(𝐀^*), ∀ i, we maximise the likelihood function ℒ_BCM=ln P_BCM(𝐀^*) with respect to them. Such a recipe leads us to solve


    k_i(𝐀^*)   =∑_j(≠ i)e^-α_i-β_j/1+e^-α_i-β_j,   ∀ i
    
    h_i(𝐀^*)   =∑_j(≠ i)e^-α_j-β_i/1+e^-α_j-β_i,   ∀ i

with obvious meaning of the symbols.



  §.§.§ Expected value of the spectral radius


According to eq. (<ref>), π_1 is the spectral radius of the N× N matrix of i.n.i.d. random variables 𝐏≡𝐏_BCM={p_ij}_i,j=1^N, with p_ij=e^-α_i-β_j/(1+e^-α_i-β_j), ∀ i≠ j.


As for the Erdös-Rényi Model, a more explicit expression can also be derived for the Binary Configuration Model. To this aim, let us consider that a way to identify π_1 in case p_ij=k_ih_j/L, ∀ i,j rests upon the relationship


    𝐏=𝐤⊗h/L=|k⟩⟨h|/L,

indicating that the matrix 𝐏 characterising the Binary Configuration Model can be obtained as the direct product of the vector of out-degrees, 𝐤, and the vector of in-degrees, 𝐡. Employing the bra-ket formalism allows the calculations to be carried out quite easily, as


    𝐏|k⟩=|k⟩⟨h|/L|k⟩=⟨h||k⟩/L|k⟩

where ⟨h||k⟩=∑_i=1^Nk_ih_i. Since 𝐏 obeys the GPF Theorem, the equation above allows us to identify the value of its spectral radius[Notice that ⟨h|𝐏=⟨h||k⟩⟨h|/L=⟨h|⟨h||k⟩/L as well.] quite straightforwardly as π_1=⟨h||k⟩/L=∑_i=1^Nk_ih_i/L. The sparse-case approximation of the Binary Configuration Model is, however, defined by the position p_ij=k_ih_j/L, ∀ i≠ j, a piece of evidence leading us to write


    π_1≃⟨h||k⟩/L=∑_i=1^Nk_ih_i/L≡λ_1^CL.




  §.§.§ Variance of the spectral radius


The expression π_1=⟨h||k⟩/L=∑_i=1^Nk_ih_i/L offers a straightforward way to calculate the variance of the spectral radius. Upon considering that all quantities defining such an expression are random variables themselves, one is led to write


    Var[λ_1]=Var[∑_i=1^Nk_ih_i/L]≡Var[λ_1^CL]

and evaluate such an expression either analytically or numerically. In what follows, we will proceed by evaluating it numerically.



 §.§ Reciprocal Configuration Model


The Reciprocal Configuration Model <cit.> is induced by the Hamiltonian 


    H(𝐀)=∑_i=1^N[α_ik_i^→(𝐀)+β_ik_i^←(𝐀)+γ_ik_i^↔(𝐀)]

where k_i^→(𝐀)=∑_j(≠ i)a_ij^→ represents the non-reciprocated out-degree of node i, k_i^←(𝐀)=∑_j(≠ i)a_ij^← represents the non-reciprocated in-degree of node i and k_i^↔(𝐀)=∑_j(≠ i)a_ij^↔ represents the reciprocated degree of node i; the vectors {α_i}_i=1^N, {β_i}_i=1^N and {γ_i}_i=1^N represent the Lagrange multipliers associated with those above, local constraints. The probability of the generic configuration 𝐀 reads


    P_RCM(𝐀)=∏_i=1^N∏_j(>i)(p_ij^→)^a_ij^→(p_ij^←)^a_ij^←(p_ij^↔)^a_ij^↔(p_ij^×)^a_ij^×

where


    p_ij^→=e^-α_i-β_j/1+e^-α_i-β_j+e^-α_j-β_i+e^-γ_i-γ_j

is the probability that a non-reciprocated link points from node i towards node j,


    p_ij^←=e^-α_j-β_i/1+e^-α_i-β_j+e^-α_j-β_i+e^-γ_i-γ_j

is the probability that a non-reciprocated link points from node j towards node i,


    p_ij^↔=e^-γ_i-γ_i/1+e^-α_i-β_j+e^-α_j-β_i+e^-γ_i-γ_j

is the probability that nodes i and j are connected by a reciprocated link and p_ij^×=1-p_ij^→-p_ij^←-p_ij^↔ is the probability that i and j are disconnected.

To tune the unknown parameters defining the Reciprocal Configuration Model to ensure that ⟨ k_i^→⟩_RCM=k_i^→(𝐀^*), ∀ i, ⟨ k_i^←⟩_RCM=k_i^←(𝐀^*), ∀ i and ⟨ k_i^↔⟩_RCM=k_i^↔(𝐀^*), ∀ i, we maximise the likelihood function ℒ_RCM=ln P_RCM(𝐀^*) with respect to them. Such a recipe leads us to solve


    k_i^→(𝐀^*)   =∑_j(≠ i)e^-α_i-β_j/1+e^-α_i-β_j+e^-α_j-β_i+e^-γ_i-γ_j,   ∀ i
    
    k_i^←(𝐀^*)   =∑_j(≠ i)e^-α_j-β_i/1+e^-α_i-β_j+e^-α_j-β_i+e^-γ_i-γ_j,   ∀ i
    
    k_i^↔(𝐀^*)   =∑_j(≠ i)e^-γ_i-γ_i/1+e^-α_i-β_j+e^-α_j-β_i+e^-γ_i-γ_j,   ∀ i

with obvious meaning of the symbols.



  §.§.§ Expected value of the spectral radius


According to eq. (<ref>), π_1 is the spectral radius of the N× N matrix of i.n.i.d. random variables 𝐏≡𝐏_RCM={p_ij}_i,j=1^N, p_ij=p_ij^→+p_ij^↔=(e^-α_i-β_j+e^-γ_i-γ_i)/(1+e^-α_i-β_j+e^-α_j-β_i+e^-γ_i-γ_j), ∀ i≠ j.


As for the Binary Configuration Model, more explicit expressions can also be derived for the Reciprocal Configuration Model. To this aim, let us consider that, in the sparse case, one can write


    𝐏^→|k^→⟩   =|k^→⟩⟨k^←|/L^→|k^→⟩=⟨k^←||k^→⟩/L^→|k^→⟩,
    𝐏^↔|k^↔⟩   =|k^↔⟩⟨k^↔|/2L^↔|k^↔⟩=⟨k^↔||k^↔⟩/2L^↔|k^↔⟩

where ⟨k^←||k^→⟩=∑_i=1^Nk^←_ik^→_i
and ⟨k^↔||k^↔⟩=∑_i=1^Nk^↔_ik^↔_i. Since 𝐏^→
and 𝐏^↔ obey the GPF Theorem, the equations above allow us to identify the values of their spectral radius[An analogous observation to the one in the previous footnote can be made.] quite straightforwardly as


    π_1^→   ≃⟨k^←||k^→⟩/L^→=∑_i=1^Nk^←_ik^→_i/L^→≡λ_1^CL^→,
    π_1^↔   ≃⟨k^↔||k^↔⟩/2L^↔=∑_i=1^Nk^↔_ik^↔_i/2L^↔≡λ_1^CL^↔

(because of the definition of the sparse-case approximation of the Reciprocal Configuration Model, valid ∀ i≠ j).



  §.§.§ Variance of the spectral radius


The expressions above offer a straightforward way to calculate the corresponding variances. In fact, one is led to write


    Var[λ_1^→]   =Var[∑_i=1^Nk^←_ik^→_i/L^→]≡Var[λ_1^CL^→],
    Var[λ_1^↔]   =Var[∑_i=1^Nk^↔_ik^↔_i/2L^↔]≡Var[λ_1^CL^↔]

and evaluate such expressions either analytically or numerically. In what follows, we will proceed by evaluating them numerically.



 §.§ Global Reciprocity Model


The Global Reciprocity Model <cit.> is a special case of the Reciprocal Configuration Model, induced by the Hamiltonian


    H(𝐀)=∑_i=1^N[α_ik_i(𝐀)+β_ih_i(𝐀)]+γ L^↔(𝐀)

where L^↔(𝐀)=∑_i=1^N∑_j(≠ i)a_ij^↔ represents the total number of reciprocated links; the parameters {α_i}_i=1^N, {β_i}_i=1^N and γ represent the Lagrange multipliers associated with the aforementioned constraints. The probability of the generic configuration 𝐀 reads


    P_GRM(𝐀)=∏_i=1^N∏_j(>i)(p_ij^→)^a_ij^→(p_ij^←)^a_ij^←(p_ij^↔)^a_ij^↔(p_ij^×)^a_ij^×

where


    p_ij^→=e^-α_i-β_j/1+e^-α_i-β_j+e^-α_j-β_i+e^-α_i-β_j-β_i-α_j-γ

is the probability that a non-reciprocated link points from node i towards j,


    p_ij^←=e^-α_j-β_i/1+e^-α_i-β_j+e^-α_j-β_i+e^-α_i-β_j-β_i-α_j-γ

is the probability that a non-reciprocated link points from node j towards node i,


    p_ij^↔=e^-α_i-β_j-β_i-α_j-γ/1+e^-α_i-β_j+e^-α_j-β_i+e^-α_i-β_j-β_i-α_j-γ

is the probability that nodes i and j are connected by a reciprocated link and p_ij^×=1-p_ij^→-p_ij^←-p_ij^↔ is the probability that i and j are disconnected.

To tune the unknown parameters defining the Global Reciprocity Model to ensure that ⟨ k_i⟩_GRM=k_i(𝐀^*), ∀ i, ⟨ h_i⟩_GRM=h_i(𝐀^*), ∀ i and ⟨ L^↔⟩_GRM=L^↔(𝐀^*), ∀ i, we maximise the likelihood function ℒ_GRM=ln P_GRM(𝐀^*) with respect to them. Such a recipe leads us to solve



    k_i(𝐀^*)   =∑_j(≠ i)e^-α_i-β_j+e^-α_i-β_j-β_i-α_j-γ/1+e^-α_i-β_j+e^-α_j-β_i+e^-α_i-β_j-β_i-α_j-γ,   ∀ i
    
    h_i(𝐀^*)   =∑_j(≠ i)e^-α_j-β_i+e^-α_i-β_j-β_i-α_j-γ/1+e^-α_i-β_j+e^-α_j-β_i+e^-α_i-β_j-β_i-α_j-γ,   ∀ i
    
    L^↔(𝐀^*)   =∑_i=1^N∑_j(≠ i)e^-α_i-β_j-β_i-α_j-γ/1+e^-α_i-β_j+e^-α_j-β_i+e^-α_i-β_j-β_i-α_j-γ


with obvious meaning of the symbols.

In the case of the Global Reciprocity Model, π_1 is the spectral radius of the N× N matrix of i.n.i.d. random variables 𝐏≡𝐏_GRM={p_ij}_i,j=1^N, p_ij=p_ij^→+p_ij^↔=(e^-α_i-β_j+e^-α_i-β_j-β_i-α_j-γ)/(1+e^-α_i-β_j+e^-α_j-β_i+e^-α_i-β_j-β_i-α_j-γ), ∀ i≠ j.



 §.§ Density-Corrected Gravity Model


The density-corrected Gravity Model <cit.> is a two-step model inducing a probability for the generic configuration 𝐀 reading


    P_dcGM(𝐀)=∏_i=1^N∏_j(≠ i)p_ij^a_ij(1-p_ij)^1-a_ij

where




    p_ij=za_il_j/1+za_il_j

is the probability that a link points from node i towards node j and a_i=∑_j(≠ i)w_ij is the out-strength of node i, l_i=∑_j(≠ i)w_ji is the in-strength of node i and z is a free parameter, determined by fixing the value of the total number of links[Analogously, one could have fixed the connectance, or link density, defined as c=L/N(N-1).], i.e. by solving the equation


    L(𝐀^*)=∑_i=1^N∑_j(≠ i)za_il_j/1+za_il_j.


The second step of the density-corrected Gravity Model, instead, is a conditional one, prescribing loading the link a_ij=1 with the value


    w_ij=a_il_j/Wp_ij,

where W=∑_i=1^N∑_j(≠ i)w_ij=∑_i=1^Na_i=∑_i=1^Nl_i is the total network volume. As a consequence of such a prescription, one recovers the result


    ⟨ w_ij⟩=a_il_j/W;

in other words, the dcGM ensures that the (financial equivalent of the) Gravity Model prescription is recovered on average.



  §.§.§ Expected value of the spectral radius


According to eq. (<ref>), ϕ_1 is the spectral radius of the N× N matrix of i.n.i.d. random variables 𝐐≡𝐐_dcGM={⟨ w_ij⟩}_i,j=1^N, ⟨ w_ij⟩=a_il_j/W, ∀ i≠ j.


As for the Binary Configuration Model, a more explicit expression can also be derived for the density-corrected Gravity Model. To this aim, let us consider that a way to identify ϕ_1 in case ⟨ w_ij⟩=a_il_j/W, ∀ i,j rests upon the relationship


    𝐐=𝐚⊗l/W=|a⟩⟨l|/W,

indicating that the matrix 𝐐 characterising the dcGM can be obtained as the direct product of the vector of out-strengths, 𝐚, and the vector of in-strengths, 𝐥. Employing the bra-ket formalism allows the calculations to be carried out quite easily, as


    𝐐|a⟩=|a⟩⟨l|/W|a⟩=⟨a||l⟩/W|a⟩

where ⟨a||l⟩=∑_i=1^Na_il_i. Since 𝐐 obeys the GPF Theorem, the equation above allows us to identify the value of its spectral radius[Notice that ⟨l|𝐐=⟨l||a⟩⟨l|/W=⟨l|⟨a||l⟩/W as well.] quite straightforwardly as ϕ_1=⟨a||l⟩/W=∑_i=1^Na_il_i/W. The density-corrected Gravity Model is, however, defined by the position ⟨ w_ij⟩=a_il_j/W, ∀ i≠ j, a piece of evidence leading us to write


    ϕ_1≃⟨a||l⟩/W=∑_i=1^Na_il_i/W≡ω_1^CL.


As the considered matrix is deterministic, the variance of its spectral radius is, by definition, zero.



§ DATA DESCRIPTION




 §.§ Dutch Interbank Network


The Dutch Interbank Network (DIN) is represented as a binary, directed network whose nodes are anonymised, Dutch banks and links represent exposures (from contractual obligations to swaps) up to one year and larger than 1.5 millions of euros. Data are reported quarterly from 1998Q1 to 2008Q4, hence consisting of 44 snapshots. Notice that the last four ends of quarters correspond to 2008, i.e. the first year of the global financial crisis <cit.>. Given the nature of the available data, a link pointing from bank i to bank j at time t indicates the existence of a total exposure of more than 1.5 million euros, directed from i to j, registered at the end of the particular quarter t.





 §.§ Electronic Market for Interbank Deposit


The Electronic Market for Interbank Deposit (e-MID) is represented as a weighted, directed network whose nodes are anonymised, Italian banks and weights represent exposures in million euros[e-MID is a centralised interbank market for trading unsecured deposits, working as follows: a bank quotes an offer to lend or borrow money (minimum quote: 1.5 million euros) at a certain maturity and interest rate; a second bank chooses (at least a part of) the quoted order (minimum quote: 50.000 euros), and the trade is registered if and only if both counterparties have agreed on it. The following information is available for each active bank during the period: an anonymous ID identifying the bank and the country where it is legally settled. In <cit.>, Fricke and Lux have highlighted i) how the number of active, foreign banks largely varies over the considered period, experiencing a dramatic drop in correspondence of the Lehman-Brothers bailout; ii) how the number of active Italian banks is quite stable over the period, although it decreases after the global financial crisis.]. Reported data cover the period January 1999-December 2014, on a daily frequency: a link with weight w_ij, pointing from bank i to bank j at time t indicates the existence of the total exposure w_ij≥ 50.000 euros, directed from i to j, registered at the end of the particular period t. Considering that ≃98% of banks are Italian and that the volume of their transactions covers ≃85% of the total volume (as of 2011), our analysis solely focuses on the subgraph induced by such a subset of nodes. We also examine all aggregation periods ranging from daily to yearly - although the figures will depict e-MID on a quarterly basis.



 §.§ International Trade Network


The International Trade Network (ITN) is represented as a weighted, directed network whose nodes are countries and weights represent imports/exports in million euros. Data on yearly trade flows during the period 2000-2020 have been downloaded from the UN-COMTRADE website[https://comtradeplus.un.org/https://comtradeplus.un.org/]. To consistently compare data, a panel of 112 countries for which trade information was available for the entire period has been selected <cit.>. Given the nature of the available data, a link whose weight is w_ij, pointing from country i to country j during the year y indicates the existence of an exported amount of commodities whose value matches w_ij, directed from i to j, during that year.



§ RESULTS




 §.§ Inspecting the accuracy of our approximations


The derivation of our results rests upon several approximations whose accuracy must be explicitly checked case by case.


The first one concerns the expected value of the trace of the exponential of 𝐀 - which has been proven to satisfy the relationship ⟨Tr[e^𝐀]⟩≥Tr[e^𝐏], hence being strictly larger than the trace of the exponential of 𝐏 for any network with positive reciprocity, i.e. having r=L^↔/L>0. In order to check how close the two terms above are, we have explicitly computed the ratio Tr[e^𝐏]/⟨Tr[e^𝐀]⟩ for all the snapshots of our systems. The results are reported in the seventh column of tables <ref> and <ref> in Appendix AppEE. As evident, Tr[e^𝐏]/⟨Tr[e^𝐀]⟩≲1 irrespectively from the structural details of our configurations - in particular, even for configurations with a non-negligible level of reciprocity such as those constituting the DIN, for which r≃ 0.3. In other words, the trace of the matrix 𝐏 describing a random network model provides a quite accurate approximation of the expected value of the trace of the adjacency matrix 𝐀 under the same model. As the 2008Q1, 2008Q2, 2008Q3 and 2008Q4 snapshots of the DIN confirm, the accuracy of the approximation above increases as r decreases.

Analogously, Tr[e^𝐐]/⟨Tr[e^𝐖]⟩≲1, as the seventh column of table <ref> in Appendix AppEE shows.


The second one concerns the hypothesis that the trace of the exponential of 𝐀 and the trace of the exponential of 𝐏 are both dominated by their largest addendum, i.e. Tr[e^𝐀]≳ e^λ_1 and Tr[e^𝐏]≳ e^π_1. In order to check how close the two pairs above of terms are, we have explicitly computed the ratios e^λ_1/Tr[e^𝐀] and e^π_1/Tr[e^𝐏] for all the snapshots of our systems. The results are reported in the fifth and sixth columns of tables <ref> and <ref> in Appendix AppEE. As evident, e^λ_1/Tr[e^𝐀]≲1 and e^π_1/Tr[e^𝐏]≲1 irrespectively from the structural details of our configurations. In words, the trace of the matrix 𝐀 is exponentially dominated by the addendum e^λ_1 and the trace of the matrix 𝐏 is exponentially dominated by the addendum e^π_1. The accuracy of the approximation remains steadily high.

Analogously, e^ω_1/Tr[e^𝐖]≲1 and e^ϕ_1/Tr[e^𝐐]≲1, as the fifth and sixth column of table <ref> in Appendix AppEE show.



 §.§ Expected value and variance of the
spectral radius


After having checked the goodness of our approximations, let us investigate the accuracy of the estimations of the expected value and variance of the spectral radius of our random network models.




  
*Erdös-Rényi Model. As the last column of tables <ref> and <ref> shows, the expected value of the spectral radius of 𝐀, evaluated numerically as the average over |𝔸|=10^3 configurations reading


    ⟨λ_1⟩=∑_𝐀∈𝔸λ_1(𝐀)/|𝔸|,

is always very well approximated by the spectral radius of 𝐏, i.e. π_1. The accuracy of such an estimation is pictorially confirmed by the left panels of fig. <ref>, showing the related scatter plot for each of the 44 snapshots constituting the DIN and for each of the 64 snapshots constituting the quarterly e-MID.



The central panels of fig. <ref>, instead, provide information about the explicit functional form of π_1, that matches the estimation reading λ_1^ER=(N-1)p=L/N.

The right panels of fig. <ref> provide information about the explicit functional form of the variance of the spectral radius by comparing


    Var[λ_1]=∑_𝐀∈𝔸[λ_1(𝐀)-⟨λ_1⟩]^2/|𝔸|

with Var[λ_1^ER]=p(1-p): as it can be appreciated, such an expression slightly underestimates the ensemble variance of the spectral radius.




  
*Binary Configuration Model. As the last column of tables <ref> and <ref> shows, the expected value of the spectral radius of 𝐀, evaluated numerically as the average over |𝔸|=10^3 configurations reading ⟨λ_1⟩=∑_𝐀∈𝔸λ_1(𝐀)/|𝔸|, is always very well approximated by the spectral radius of 𝐏, i.e. π_1. The accuracy of such an estimation is pictorially confirmed by the left panels of fig. <ref>, showing the related scatter plot for each of the 44 snapshots constituting the DIN and for each of the 64 snapshots constituting the quarterly e-MID.

The central panels of fig. <ref>, instead, provide information about the explicit functional form of π_1 which is (overall) well approximated by the Chung-Lu estimation reading λ_1^CL=∑_i=1^Nk_ih_i/L for what concerns the e-MID and overestimated by the same expression for what concerns the DIN.

The right panels of fig. <ref> provide information about the explicit functional form of the variance of the spectral radius, by comparing Var[λ_1]=∑_𝐀∈𝔸[λ_1(𝐀)-⟨λ_1⟩]^2/|𝔸| with


    Var[λ_1^CL]=∑_𝐀∈𝔸[λ_1^CL(𝐀)-⟨λ_1^CL⟩]^2/|𝔸|;

as it can be appreciated, such an expression either underestimates (for what concerns the e-MID) or overestimates (for what concerns the DIN) the ensemble variance of the spectral radius. Notice also that such an expression calculates the variance of the spectral radius by evaluating λ_1^CL(𝐀), i.e. the numerical value of the Chung-Lu approximation, for each matrix in the sampled ensemble. As fig. <ref> shows, these discrepancies seem to be due to a systematic mismatch caused by the configuration-specific values of the spectral radius - the DIN, for instance, obeys the relationship λ_1^CL(𝐀)>λ_1(𝐀), ∀ 𝐀, a result potentially explaining the differences between λ_1^CL and π_1 and between Var[λ_1^CL] and Var[λ_1] - in words, the numbers λ_1^CLs are not only larger than their ensemble counterparts but are also more dispersed (see also fig. <ref> in Appendix AppFF).




  
*Reciprocal Configuration Model. The Reciprocal Configuration Model performs similarly to the Binary Configuration Model. While the last column of tables <ref> and <ref> shows that the expected value of the spectral radius of 𝐀, evaluated numerically as the average over |𝔸|=10^3 configurations reading ⟨λ_1⟩=∑_𝐀∈𝔸λ_1(𝐀)/|𝔸|, is always very well approximated by the spectral radius of 𝐏, i.e. π_1, the left panels of fig. <ref>, show the scatter plot concerning the two sets of quantities ⟨λ_1^↔⟩ and π_1^↔ for each of the 44 snapshots constituting the DIN and for each of the 64 snapshots constituting the quarterly e-MID.

The central panels of fig. <ref>, instead, provide information about the explicit functional form of π_1^↔ which is (overall) well approximated by the Chung-Lu estimation reading λ_1^CL=∑_i=1^Nk_i^↔ k_i^↔/2L for what concerns the e-MID and overestimated by the same expression for what concerns the DIN.



The right panels of fig. <ref> provide information about the explicit functional form of the variance of the spectral radius by comparing Var[λ_1^↔]=∑_𝐀∈𝔸[λ_1^↔(𝐀)-⟨λ_1^↔⟩]^2/|𝔸| with


    Var[λ_1^CL^↔]=∑_𝐀∈𝔸[λ_1^CL^↔(𝐀)-⟨λ_1^CL^↔⟩]^2/|𝔸|;

as it can be appreciated, such an expression overestimates the ensemble variance of the spectral radius. As for the Binary Configuration Model, such an expression calculates the variance of the spectral radius by evaluating λ_1^CL^↔(𝐀), i.e. the numerical value of the Chung-Lu approximation, for each matrix in the sampled ensemble. These discrepancies may, thus, be imputable to a systematic mismatch caused by the configuration-specific values of the spectral radius.




  
*Density-Corrected Gravity Model. The last column of tables <ref> and <ref> shows that the expected value of the spectral radius of 𝐖, evaluated numerically as the average over |𝕎|=10^3 configurations reading ⟨ω_1⟩=∑_𝐖∈𝕎ω_1(𝐖)/|𝕎|, is always very well approximated by the spectral radius of 𝐐, i.e. ϕ_1, as the left panels of fig. <ref> pictorially confirm. Besides, the right panels of the same figure provide information about the explicit functional form of ϕ_1 which is (overall) well approximated by the Chung-Lu estimation reading ω_1^CL=∑_i=1^Na_il_i/W for each of the 16 snapshots constituting the yearly e-MID and for each of the 21 snapshots constituting the yearly ITN.





 §.§ Spectral signature of structural changes in financial networks


Now, let us inspect the presence of structural changes affecting our networked configurations. To this aim, we will plot the evolution of z[λ_1] across the periods covered by our datasets; we will proceed numerically by explicitly sampling the network ensemble induced by each of the benchmarks considered here per snapshot.



  §.§.§ Dutch Interbank Network


As fig. <ref> clearly shows, the structural change undergone by the system in 2008 is signalled by several quantities: the total number of active Dutch banks sharply decreases as well as the total number of links, whose number diminishes in corresponding of the last year covered by our dataset; this, in turn, causes the connectance to rise. As already discussed in <cit.>, one of the most evident signals of the global financial crisis is provided by reciprocity: for most of the period, it is characterised by an essentially constant trend, with small fluctuations around an average value of ≃0.26; the last, four snapshots are, then, characterised by a drop of ≃40%, causing the empirical values to lie almost three sigmas away from the sample average - a trend indicating that the reciprocity of the DIN is anomalously low during the critical period and imputable to a decrease of the level of trust characterising the Dutch system.

An additional signal of the global financial crisis is provided by the empirical value of the spectral radius itself, which decreases in correspondence with 2008Q1 and remains constant across the last four snapshots of our dataset. As it is related to the number of closed walks in a network, its decrease may be related to the decline of reciprocity. However, the latter's trend appears as (much) less affected by the statistical fluctuations characterising the evolution of the DIN throughout its entire history.

Let us now comment on the signal provided by z[λ_1]. Even if the Erdös-Rényi Model is, from a merely financial perspective, an unlikely benchmark (its homogeneous nature forces the banks to be similar in size), employing it still allows us to conclude that the DIN is characterised by two structural changes - the first one taking place across 2005 and the second one taking place across 2008. More specifically, after a (more or less) stationary trend characterising the evolution of the DIN from 1998 to 2005 - in correspondence of which the number of closed walks is significantly large - a smooth trend characterising the pre-crisis phase is recovered; afterwards, an abrupt drop connecting the last quarter of 2007 with the first quarter of 2008 emerges. Such a result complements the ones presented in <cit.> where such behaviour could have been revealed only by employing a heterogeneous benchmark (specifically, the Binary Configuration Model).

Employing the heterogeneous benchmarks - preserving the heterogeneity of banks by constraining the observed (reciprocal) degrees - leads to the same qualitative result. More quantitatively, instead, all such null models reveal that the number of closed walks is perfectly compatible with their predictions during the stationary phase of the system. Such a consistency confirms that, in the absence of distress, the topology of the DIN can be reconstructed quite accurately, solely employing the information provided by the number of (inward, outward and reciprocated) partners of each bank. It is noticed that the explanatory power of the Reciprocal Configuration Model is larger than that of the Global Reciprocity Model, which, in turn, is (only slightly) larger than that of the Binary Configuration Model.

As the build-up phase of the crisis began, a decreasing trend led to 2008, indicating that the local connectivity of banks became less and less informative about the network as a whole - emerges. Under the same benchmarks, the second regime shift is preceded by a short, rising trend. As already noticed in <cit.>, maximum-entropy techniques yield a realistic guess of the real network only in tranquil times: when the network is under stress, instead, these models provide a sort of distorted picture of it, whose differences from the empirical situation constitute the structural changes we are looking for.

Apart from model-specific differences, however, the degree of informativeness about the changes affecting the DIN carried by the spectral radius seems quite independent of the model employed to spot the differences above.



  §.§.§ Electronic Market for Interbank Deposit


For what concerns the e-MID, instead, the evolution of the total number of active Italian banks steadily decreases, hence not providing any clear indication about the presence of structural changes. On the contrary, the evolution of the total number of links provides a quite clear indication of the presence of two regime shifts as L drops in correspondence of 2008 and 2012. Overall, the connectance and the reciprocity provide a very similar indication - the global financial crisis being characterised by a stronger signal than the one characterising the long-term refinancing operation (LTRO) promoted by the European Central Bank at the end of 2011[The two LTRO measures date December the 22nd, 2011 and February the 29th, 2012.].

The evolution of the empirical value of the spectral radius is characterised by a drop in correspondence of the first crisis, originating a slightly fluctuating trend that lasts until 2012, the year in correspondence of which a second, decreasing trend can also be observed.



Let us now comment on the signal provided by z[λ_1]. Employing a homogeneous benchmark such as the Erdös-Rényi Model allows us to conclude that the e-MID is characterised by three structural changes, the first one taking place across 2000, the second one taking place between 2007 and 2008 and the third one taking place across 2012.

More specifically, the evolution of the e-MID starts with a drop of the z-score of the spectral radius, indicating that the number of closed walks has become significantly smaller than expected during 2001. Afterwards, an increasing trend leading to a phase characterised by several closed walks compatible with the output of the prediction by the Erdös-Rényi Model becomes visible. Such a period is interrupted by the so-called pre-crisis phase, during which the trend of λ_1 reverts and becomes again significantly smaller than expected. From 2009 on, a second, increasing trend lasting until 2012 becomes visible: afterwards, the system stabilises.

Employing the heterogeneous benchmarks leads to quite different results: more quantitatively, the first regime shift disappears, replaced by a stationary trend lasting until 2003; afterwards, a rising trend leading the system to its (pre-)critical phase appears. Since 2009 on, a decreasing trend lasting a couple of years emerges to be followed, once more, by an increasing one. From this perspective, the DIN and the e-MID behave, somehow, oppositely: while the global financial crisis induces a statistically significant signal in the case of the DIN, it does not in the case of the e-MID. In a sense, maximum-entropy techniques can be used to reconstruct the e-MID when the system is under stress, while this should be avoided in tranquil times - e.g. the first years of the dataset - when the picture of it inferred from local constraints departs the most from the empirical one.

Differently from the DIN, the explanatory power of the Reciprocal Configuration Model (still larger than the one of the Global Reciprocity Model, which, however, performs similarly to the Binary Configuration Model) is so large that the measurements carried out on the e-MID (practically) always compatible with the predictions. Although such a piece of evidence speaks against the use of the Reciprocal Configuration Model to detect deviations from the average behaviour, statistical tendencies can still be revealed, confirming once more that a dichotomous yes/no answer to the question is this pattern statistically significant? may be quite unsatisfactory to gain a sufficiently deep insight into system behaviour.



§ DISCUSSION


The so-called stability analysis represents an application of particular interest in the study of financial networks, a topic whose popularity has steadily increased since the turmoil due to the mortgage crisis <cit.>. The objective of this kind of analysis is to understand the relationship(s) between the topological structure of financial networks and their resilience to events like shocks, cascading failures, etc., by employing real data <cit.>, reconstructed configurations <cit.> or (simple) toy models <cit.>. A direct way to explore this connection is by running stress tests on several different topological structures by measuring the effects of a simulated shock and the subsequent propagation of losses ex post <cit.>: later works have related these results to the magnitude of the spectral radius of the so-called leverage matrix <cit.> although no algorithm has been devised to estimate its magnitude from the (partial) information that is usually available in financial contexts.

With the present contribution, we have tackled a more general challenge, i.e. that of estimating the spectral radius of random network models calibrated on real-world evolving networks. To this aim, we have adopted several approximations that have led to the surprisingly simple recipe ⟨λ_1⟩≃π_1 for estimating the expected value of λ_1, with π_1 representing the spectral radius of the probabilistic matrix describing the chosen model. Despite our result is based on an approximation[We have explicitly verified that the properties of existence, reality, positivity, maximality and uniqueness of the spectral radius hold for each, considered configuration.], it turns out to be extremely accurate for any directed (binary or weighted) random network model considered.

Besides the theoretical relevance of such a result, its usefulness lies in spotting the structural changes separating a (financial) regime from another by exploiting the interplay between distress and topological changes. As the case studies of the DIN and the e-MID illustrate, deviations from the average behaviour can happen in both directions, either moving away from a less structured configuration (hence becoming a less typical member of an equilibrium ensemble of graphs) or moving towards a less structured configuration (hence becoming a more typical member of an equilibrium ensemble of graphs): from this perspective, each quantity characterising the original network can be straightforwardly assigned a level of significance - which is sensitive to the direction - by computing the related z-score, i.e. an index comparing the measured value with the one expected under a null model preserving some properties of the observed network but, otherwise, being maximally random.

Although our results become exact in case a perfectly non-reciprocal network is observed, future research calls for a more accurate evaluation of our approximations - hopefully, in terms of the reciprocity itself. Besides, extending the results of the present analysis to undirected, binary or weighted networks would enlarge their applicability beyond the economic and financial domains.



§ ACKNOWLEDGMENTS


VM acknowledges support from the project NetRes - `Network analysis of economic and financial resilience', Italian DM n. 289, 25-03-2021 (PRO3 2021-2023 University joint program `Le Scuole Superiori ad Ordinamento Speciale: istituzioni a servizio del Paese'), CUP D67G22000130001 (<https://netres.imtlucca.it>) funded by the Italian Ministry of University and Research (MUR). This work is also supported by the European Union - NextGenerationEU - National Recovery and Resilience Plan (PNRR Research Infrastructures), project `SoBigData.it - Strengthening the Italian RI for Social Mining and Big Data Analytics' - Grant IR0000013 (Avviso MUR D.D. n. 3264, 28/12/2021)  (<https://pnrr.sobigdata.it/>) and the project “Reconstruction, Resilience and Recovery of Socio-Economic Networks” RECON-NET EP_FAIR_005 - PE0000013 “FAIR” - PNRR M4C2 Investment 1.3, financed by the European Union – NextGenerationEU.

We thank Anna Gallo for useful discussions.







AppA


§ APPENDIX A.
DYADIC EARLY-WARNING SIGNALS


Upon defining 


    X   =∑_i=1^N∑_j=1^Na_ija_ji=∑_i=1^N[𝐀^2]_ii=Tr[𝐀^2],

we are left with the task of calculating its expected value and variance. The evidence that the expected value is a linear operator (i.e. ⟨ aX+bY⟩=a⟨ X⟩+b⟨ Y⟩) and that the entries of a binary, directed network are treated as independent random variables under any of the random network models considered here, makes such a calculation straightforward. In fact,


    ⟨ X⟩=⟨∑_i=1^N∑_j=1^Na_ija_ji⟩=∑_i=1^N∑_j=1^N⟨ a_ija_ji⟩=∑_i=1^N∑_j=1^N⟨ a_ij⟩⟨ a_ji⟩=∑_i=1^N∑_j=1^Np_ijp_ji.


In order to calculate the variance of X, let us consider that X can be re-written as


    X=∑_i=1^N∑_j=1^Na_ija_ji=2∑_i=1^N∑_j(>i)a_ija_ji

i.e. as a sum over dyads, treated as independent random variables under any random network models considered here. Since the variance of a sum of independent random variables coincides with the sum of their variances, one can write


    Var[X]   =Var[2∑_i=1^N∑_j(>i)a_ija_ji]=4·∑_i=1^N∑_j(>i)Var[a_ija_ji];

then, since a_ija_ji∼Ber[p_ijp_ji], one finds that


    Var[X]   =4·∑_i=1^N∑_j(>i)p_ijp_ji(1-p_ijp_ji).


It is nevertheless instructive to follow an alternative road and consider that


    Var[X]=Var[∑_i=1^N∑_j=1^Na_ija_ji]   =∑_i=1^N∑_j=1^NVar[a_ija_ji]+2·∑_i=1^N∑_j(>i)Cov[a_ija_ji,a_ija_ji]
       =∑_i=1^N∑_j=1^Np_ijp_ji(1-p_ijp_ji)+2·∑_i=1^N∑_j(>i)p_ijp_ji(1-p_ijp_ji)
       =2·∑_i=1^N∑_j(>i)p_ijp_ji(1-p_ijp_ji)+2·∑_i=1^N∑_j(>i)p_ijp_ji(1-p_ijp_ji)
       =4·∑_i=1^N∑_j(>i)p_ijp_ji(1-p_ijp_ji).


The comparison between the analytical estimations of the expected value and the variance of the number of dyads and the numerical counterparts, obtained by explicitly sampling the ensembles induced by the Erdös-Rényi Model and the Binary Configuration Model is illustrated in figs. <ref> and <ref>.











AppB


§ APPENDIX B.
TRIADIC EARLY-WARNING SIGNALS


Upon defining 


    X   =∑_i=1^N∑_j=1^N∑_k=1^Na_ija_jka_ki=∑_i=1^N[𝐀^3]_ii=Tr[𝐀^3],

we are left with the task of calculating its expected value and variance. Analogously to the dyadic case, calculating the expected value is straightforward. In fact,


    ⟨ X⟩=⟨∑_i=1^N∑_j=1^N∑_k=1^Na_ija_jka_ki⟩=∑_i=1^N∑_j=1^N∑_k=1^N⟨ a_ija_jka_ki⟩=∑_i=1^N∑_j=1^N∑_k=1^N⟨ a_ij⟩⟨ a_jk⟩⟨ a_ki⟩=∑_i=1^N∑_j=1^N∑_k=1^Np_ijp_jkp_ki.


In order to calculate the variance of X, let us, first, consider that X can be re-written as


    X   =∑_i=1^N∑_j=1^N∑_k=1^Na_ija_jka_ki=3·∑_i=1^N∑_j(>i)∑_k(>j)(a_ija_jka_ki+a_ika_kja_ji)≡3·∑_i<j<k(a_ija_jka_ki+a_ika_kja_ji)

i.e. as a sum over triads. Then, let us notice that


    Var[X]=3^2·[∑_𝐈Var[a_𝐈]+2·∑_𝐈<𝐉Cov[a_𝐈,a_𝐉]]

where we have employed the multi-index notation, i.e. 𝐈≡(i,j,k) and 𝐉≡(l,m,n). More explicitly,


    Var[a_𝐈]   =Var[a_ija_jka_ki]+Var[a_ika_kja_ji]+Cov[a_ija_jka_ki,a_ika_kja_ji]
       =p_ijp_jkp_ki(1-p_ijp_jkp_ki)+p_ikp_kjp_ji(1-p_ikp_kjp_ji)+Cov[a_ija_jka_ki,a_ika_kja_ji]

with Cov[a_ija_jka_ki,a_ika_kja_ji] depending on the adopted benchmark: under both the Erdös-Rényi Model and the Binary Configuration Model, it amounts at zero. Overall, thus,


    ∑_𝐈Var[a_𝐈]=∑_i<j<k[p_ijp_jkp_ki(1-p_ijp_jkp_ki)+p_ikp_kjp_ji(1-p_ikp_kjp_ji)].


Moreover,


    Cov[a_𝐈,a_𝐉]   =⟨(a_ija_jka_ki+a_ika_kja_ji)·(a_lma_mna_nl+a_lna_nma_ml)⟩-⟨ a_ija_jka_ki+a_ika_kja_ji⟩·⟨ a_lma_mna_nl+a_lna_nma_ml⟩
       =⟨(a_ija_jka_ki+a_ika_kja_ji)·(a_lma_mna_nl+a_lna_nma_ml)⟩-(p_ijp_jkp_ki+p_ikp_kjp_ji)·(p_lmp_mnp_nl+p_lnp_nmp_ml)

is different from zero, i.e. any two triads co-variate as long as they share an edge. In this case, they form a diamond whose vertices can be labelled either as i≡ l, j≡ m, k, n or as i≡ m, j≡ l, k, n and induce the expression


    Cov[a_𝐈,a_𝐉]   =p_ijp_jkp_kip_jnp_ni-(p_ij)^2p_jkp_kip_jnp_ni+p_jip_ikp_kjp_inp_nj-(p_ji)^2p_ikp_kjp_inp_nj
       =p_ij(1-p_ij)p_jkp_kip_jnp_ni+p_ji(1-p_ji)p_ikp_kjp_inp_nj.


Let us now, calculate the number of times such an expression appears, i.e. the number of triples sharing an edge: since we need to first, choose the pair of nodes individuating the common edge and, then the pair of nodes individuating the `free' vertices of the two triads, such a number amounts at N2N-22=N(N-1)(N-2)(N-3)/4; in case N=4, it amounts at 3!=6 - indeed, let us concretely focus on the triads (1,2,3), (1,2,4), (1,3,4), (2,3,4): (1,2,3) co-variates with (1,2,4), (1,3,4), (2,3,4); (1,2,4) co-variates with (1,3,4), (2,3,4); (1,3,4) co-variates with (2,3,4). Overall, then,


    ∑_𝐈<𝐉Cov[a_𝐈,a_𝐉]=3!·∑_i<j<k<n[p_ij(1-p_ij)p_jkp_kip_jnp_ni+p_ji(1-p_ji)p_ikp_kjp_inp_nj].


The comparison between the analytical estimations of the expected value and the variance of the number of triads and the numerical counterparts, obtained by explicitly sampling the ensembles induced by the Erdös-Rényi Model and the Binary Configuration Model is illustrated in figs. <ref> and <ref>.











AppC


§ APPENDIX C.
DIAGONALISATION AND TRACE OF THE MATRIX EXPONENTIAL


In this Appendix, we will provide a sketch of the proof that


    f(𝐀)=𝐅f(Λ)𝐅^-1

and that


    Tr[f(𝐀)]=Tr[𝐅f(Λ)𝐅^-1]=Tr[f(Λ)𝐅^-1𝐅]=Tr[f(Λ)],

i.e. that the trace is invariant under a cyclic permutation of matrices, in the special case f(·)≡ e^(·) and where 𝐅 is the matrix that diagonalises 𝐀, i.e. the one ensuring that 𝐅^-1𝐀𝐅=Λ.


Since the function of a matrix is formally identical to its series expansion, one can write that


    e^𝐀≡𝐈+𝐀+𝐀^2/2!+𝐀^3/3!+…+𝐀^n/n!+…;

let us now diagonalise it:


    𝐅^-1e^𝐀𝐅   ≡   𝐅^-1𝐈𝐅+𝐅^-1𝐀𝐅+𝐅^-1𝐀^2𝐅/2!+𝐅^-1𝐀^3𝐅/3!+…+𝐅^-1𝐀^n𝐅/n!+…
       =   𝐈+Λ+(𝐅^-1𝐀𝐅)(𝐅^-1𝐀𝐅)/2!+(𝐅^-1𝐀𝐅)(𝐅^-1𝐀𝐅)(𝐅^-1𝐀𝐅)/3!+…
       =   𝐈+Λ+(𝐅^-1𝐀𝐅)^2/2!+(𝐅^-1𝐀𝐅)^3/3!+…+(𝐅^-1𝐀𝐅)^n/n!+…
       =   𝐈+Λ+Λ^2/2!+Λ^3/3!+…+Λ^n/n!+…≡ e^Λ.


Since all matrices appearing in the last row are diagonal, e^Λ also has diagonal entries. As a consequence,


    Tr[e^Λ]   =   ∑_i=1^N(e^Λ)_ii=∑_i=1^Ne^λ_i=∑_i=1^N1+∑_i=1^Nλ_i+∑_i=1^Nλ_i^2/2!+∑_i=1^Nλ_i^3/3!+…+∑_i=1^Nλ_i^n/n!+…
       =   Tr[𝐈]+Tr[Λ]+Tr[Λ^2]/2!+Tr[Λ^3]/3!+…+Tr[Λ^n]/n!+…
       =   Tr[𝐈]+Tr[𝐅^-1𝐀𝐅]+Tr[(𝐅^-1𝐀𝐅)(𝐅^-1𝐀𝐅)]/2!+Tr[(𝐅^-1𝐀𝐅)(𝐅^-1𝐀𝐅)(𝐅^-1𝐀𝐅)]/3!+…
       =   Tr[𝐈]+Tr[𝐅^-1𝐀𝐅]+Tr[𝐅^-1𝐀^2𝐅]/2!+Tr[𝐅^-1𝐀^3𝐅]/3!+…+Tr[𝐅^-1𝐀^n𝐅]/n!+…
       =   Tr[𝐈]+Tr[𝐀𝐅𝐅^-1]+Tr[𝐀^2𝐅𝐅^-1]/2!+Tr[𝐀^3𝐅𝐅^-1]/3!+…+Tr[𝐀^n𝐅𝐅^-1]/n!+…
       =   Tr[𝐈]+Tr[𝐀]+Tr[𝐀^2]/2!+Tr[𝐀^3]/3!+…+Tr[𝐀^n]/n!+…=Tr[e^𝐀],

where we have exploited the property of the trace of being invariant under circular shifts.



AppD


§ APPENDIX D.
ENSEMBLE DISTRIBUTION OF THE SPECTRAL RADIUS






AppE


§ APPENDIX E.
DUTCH INTERBANK NETWORK








§ ELECTRONIC MARKET FOR INTERBANK DEPOSIT








§ INTERNATIONAL TRADE NETWORK










AppF


§ APPENDIX F.
INSPECTING THE ACCURACY OF THE CHUNG-LU APPROXIMATION




