Salem Said

Simon Heuveline

Cyrus Mostajeran

Riemannian statistics meets random matrix theory : towards learning from high-dimensional covariance matrices

Keywords: Riemannian Gaussian distribution, random matrix theory, covariance matrix, block-Toeplitz covariance matrix, high-dimensional data I

Riemannian Gaussian distributions were initially introduced as basic building blocks for learning models which aim to capture the intrinsic structure of statistical populations of positive-definite matrices (here called covariance matrices). While the potential applications of such models have attracted significant attention, a major obstacle still stands in the way of these applications : there seems to exist no practical method of computing the normalising factors associated with Riemannian Gaussian distributions on spaces of high-dimensional covariance matrices. The present paper shows that this missing method comes from an unexpected new connection, with random matrix theory.

Its main contribution is to prove that Riemannian Gaussian distributions of real, complex, or quaternion covariance matrices are equivalent to orthogonal, unitary, or symplectic log-normal matrix ensembles. This equivalence yields a highly efficient approximation of the normalising factors, in terms of a rather simple analytic expression. The error due to this approximation decreases like the inverse square of dimension. Numerical experiments are conducted which demonstrate how this new approximation can unlock the difficulties which have impeded applications to realworld datasets of high-dimensional covariance matrices. The paper then turns to Riemannian Gaussian distributions of block-Toeplitz covariance matrices. These are equivalent to yet another kind of random matrix ensembles, here called "acosh-normal" ensembles. Orthogonal and unitary "acosh-normal" ensembles correspond to the cases of block-Toeplitz with Toeplitz blocks, and block-Toeplitz (with general blocks) covariance matrices, respectively.

Over the past few years, Riemannian Gaussian distributions, and mixtures of these distributions, were introduced as a means of modeling the intrinsic structure of statistical populations of covariance matrices [START_REF] Cheng | A novel dynamic system in the space of SPD matrices with applications to appearance tracking[END_REF]- [START_REF] Said | Gaussian distributions on Riemannian symmetric spaces : statistical learning with structured covariance matrices[END_REF]. They have found successful applications in fields such as brain-computer interface analysis and artificial intelligence [START_REF] Zanini | Transfer learning: A Riemannian geometry framework with applications to brain-computer interfaces[END_REF] [START_REF] Mathieu | Continuous hierarchical representations with Poincaré variational auto-encoders[END_REF].

However, such applications could not be pursued for covariance matrices of relatively larger dimension (for example, 50 × 50). Indeed, there seemed to exist no practical method of computing the normalising factors associated with Riemannian Gaussian distributions on spaces of high-dimensional covariance matrices.

The theoretical contribution of this paper is to show that this missing method arises quite naturally, as soon as one realises that a strong connection exists between Riemannian Gaussian distributions and random matrix theory.

Roughly speaking, Riemannian Gaussian distributions on spaces of real, complex, or quaternion covariance matrices correspond to log-normal orthogonal, unitary or symplectic matrix ensembles. These are similar to the classical, widely-known Gaussian orthogonal, unitary and symplectic ensembles, but with the normal weight function replaced with a log-normal weight function. Thanks to this new connection, the powerful tools of random matrix theory can be employed to uncover several original properties of Riemannian Gaussian distributions, especially for higherdimensional covariance matrices.

The present paper also extends the equivalence between Riemannian Gaussian distributions and random matrix ensembles to the case of block-Toeplitz matrices. Instead of a log-normal weight function exp(-log 2 (x)), the corresponding weigh function is "acosh-normal" exp(-acosh 2 (x)), where acosh is the inverse hyperbolic cosine.

No attempt is made to develop these additional "acosh-normal" ensembles, at least not for now. In fact, there are already plenty of results to derive and apply with the log-normal ensembles, which cover the important cases of real and complex covariance matrices. This is in part because log-normal ensembles (also called Sitltjes-Wigert ensembles) are amenable to a direct analytic treatment, as will be seen below.

It is also because of a rather surprising connection also noted in [START_REF] Heuveline | Gaussian distributions on Riemannian symmetric spaces in the large N limit[END_REF]. The log-normal unitary ensemble had appeared in the theoretical physics literature about twenty years ago, as a random matrix model for the Chern-Simons quantum field theory [START_REF] Mariño | Chern-Simons theory, matrix models, and topological strings[END_REF]. A seen in the present paper, this ensemble corresponds to Riemannian Gaussian distributions on the space of complex covariance matrices. Of course, such Riemannian Gaussian distributions have nothing to do with quantum field theory, but some of their valuable properties can be obtained by carefully re-adapting already existing results, found within this theory.

The present paper was developed independently from a theoretical physics paper, published only very recently [START_REF] Santilli | Riemannian Gaussian distributions, random matrix ensembles and diffusion kernels[END_REF].

Both papers focus on the connection between Riemannian Gaussian distributions and log-normal matrix ensembles, but significantly differ in terms of their contribution and focus, as will be discussed below. The reader is also referred to [START_REF] Forrester | Vicious random walkers in the limit of a large number of walkers[END_REF] for previous results on asymptotic computations of log-partition functions from the random matrix theory literature, as well as the recent paper [START_REF]Global and local scaling limits for the β= 2 stieltjes-wigert random matrix ensemble[END_REF], which considers global and local scaling limits for the β = 2

Stieltjes-Wigert random matrix ensemble and associated physical interpretations.

The paper is organised as follows. Section II recalls basic background material on Riemannian Gaussian distributions. Section III presents the main original results, obtained by applying random matrix theory to the study of Riemannian Gaussian distributions. Section IV illustrates the importance of these original results to learning from datesets of high-dimensional covariance matrices. Finally, Section V considers Riemannian Gaussian distributions on the space of block-Toeplitz covariance matrices, and introduces the corresponding "acosh-normal" ensembles.

Proofs of the propositions stated in Section III and V are provided in Appendix A.

II. NOTATION AND BACKGROUND

Let P β N denote the space of N × N covariance matrices which are either real (β = 1), complex (β = 2) or quaternion (β = 4). Precisely, the elements of P 1 N are real, symmetric positive-definite matrices, while those of P 2

N

and P 4 N are respectively complex and quaternion, hermitian positive definite matrices. Of course, P β

N is an open convex cone, sitting inside the real vector space S β N of real, complex or quaternion self-adjoint matrices (the adjoint of a matrix being its conjugate transpose). Therefore, P β N is a real differentiable manifold, whose tangent space

T Y P β N at any point Y ∈ P β N is isomorphic to S β N .
In particular, the dimension of

P β N is dim β N = N N β where N β = β 2 (N -1) + 1. The elements Y of P β
N are in one-to-one correspondence with the centred (zero-mean) N -variate normal distributions (real, complex circular, or quaternion circular, according to the value of β). Thus, P β N can be equipped with the Rao-Fisher information metric of the centred N -variate normal model [START_REF] Amari | Information geometry and its applications[END_REF]. This is identical to the so-called affine-invariant metric introduced in [START_REF] Pennec | Intrinsic statistics on Riemannian manifolds: Basic tools for geometric measurements[END_REF],

⟨u, v⟩ Y = ℜ tr ((Y -1 u) (Y -1 v)) u, v ∈ T Y P β N ≃ S β N (1)
where ℜ denotes the real part and tr denotes the trace. The Riemannian geometry of the metric ( 1) is quite wellknown in the geometric information science community (see the recent book [START_REF] Pennec | Riemannian geometric statistics in medical image analysis[END_REF]). Recall here the associated geodesic distance

d 2 (X , Y ) = tr log 2 X -1 2 Y X -1 2 X , Y ∈ P β N (2)
where matrix logarithms and powers are understood as self-adjoint matrix functions, obtained by taking logarithms and powers of eigenvalues. The main advantage of the metric (1) is its invariance under the action of the group G β N of N × N invertible matrices with real, complex or quaternion entries (according to the value of β). For example,

if A ∈ G β N and X is replaced by A • X = AXA † while Y is replaced by A • Y = AYA † ,
then the distance d(X , Y ) remains unchanged (note that † denotes the adjoint, or conjugate transpose).

In terms of the distance (2), we define Riemannian Gaussian distributions on the space P β N to be given by parameterised probability density function

p(Y | Ȳ , σ), as in [2] [3], p(Y | Ȳ , σ) = (Z(σ)) -1 exp - d 2 (Y , Ȳ ) 2σ 2 (3)
where Ȳ ∈ P β N is the centre-of-mass parameter, and σ > 0 the dispersion parameter. The normalising factor Z(σ) is given by the integral

Z(σ) = P β N exp - d 2 (Y , Ȳ ) 2σ 2 dv(Y ) (4)
with respect to the Riemannian volume dv(Y ) = det(Y ) -N β {dY }. Here,

{dY } = i≤j dY ij (β = 1) {dY } = i≤j dY (a) ij i<j dY (b) ij (β = 2) {dY } = i≤j dY (a) ij i<j dY (b) ij i<j dY (c) ij i<j dY (d) ij (β = 4)
where

Y ij = Y (a) ij + Y (b) ij i if Y ij is complex and Y ij = Y (a) ij + Y (b) ij i + Y (c) ij j + Y (d) ij k if Y ij
is quaternion (here, i, j, k denote complex or quaternion imaginary units).

Both formulae (3) and ( 4) are greatly simplified by introducing "polar coordinates". Each Y ∈ P β N can be diagonalised as Y = U e r U † where e r is a diagonal matrix, with diagonal elements e ri for (r 1 , . . . , r N ) ∈ R N , and

U ∈ K β N , which means U U † = I N (the N × N identity matrix). Note that K 1 N is the orthogonal group O(N ), K 2 N
is the unitary group U (N ), and K 4 N is the symplectic (quaternion unitary) group Sp(N ) (for deeper insight, see the recent review [START_REF] Edelman | Fifty three matrix factorizations: A systematic approach[END_REF]). Now, as in [START_REF] Said | Gaussian distributions on Riemannian symmetric spaces : statistical learning with structured covariance matrices[END_REF], let Y follow the Riemannian Gaussian density (3), and

Y ′ = Ȳ -1 2 Y Ȳ -1 2 . If Y ′ = U e r U † ,
then U is uniformly distributed on K β N (for a rigorous definition of this "uniform distribution", see [START_REF] Meckes | The random matrix theory of the classical compact groups[END_REF]), while (r 1 , . . . , r N ) have joint probability density function

p(r|σ) ∝ i exp - r 2 i 2σ 2 i<j sinh β r i -r j 2 (5)
where ∝ indicates proportionality and sinh the hyperbolic sine. In addition, the normalising factor Z(σ) in ( 4) reduces to a certain multiple integral z β (σ). Specifically,

Z(σ) = Ω N z β (σ) (6)
where Ω N is a numerical constant, which appears after the uniformly distributed matrix U is integrated out of (4), and where

z β (σ) = 1 N ! R N i<j sinh β r i -r j 2 i exp - r 2 i 2σ 2 dr i (7)
One of the main issues addressed in the present paper is the efficient approximation of the multiple integral z β (σ).

Until now (in [START_REF] Said | Riemannian Gaussian distributions on the space of symmetric positive definite matrices[END_REF] [3]), this was done using a Monte Carlo technique which involved a smoothing method (containing certain arbitrarily fixed parameters) and which failed to produce coherent results when the dimension N increased beyond N = 20.

Before proceeding to present our main results, we should briefly note the existence of several alternative proposals for the extension of Gaussian distributions to Riemannian manifolds, including constructions based on heat flows and diffusion processes [START_REF] Sommer | Anisotropic distributions on manifolds: Template estimation and most probable paths[END_REF]. Further discussion of these alternative formulations is beyond the scope of this paper and the interested reader is referred to the literature on the topic for further information. See [START_REF] Pennec | Riemannian geometric statistics in medical image analysis[END_REF] (Section 3.4.3 and Chapter 10) and the references therein for a recent and comprehensive account.

III. MAIN RESULTS

The main results of the present paper stem from the equivalence between Riemannian Gaussian distributions on the spaces of real, complex, or quaternion covariance matrices and log-normal orthogonal, unitary, and symplectic matrix ensembles.

Let Y be a random matrix in P β N which follows the Riemannian Gaussian density p(Y | Ȳ , σ) of (3). It is always possible to assume that Ȳ = I N , since this holds after replacing Y by

Y ′ = Ȳ -1 2 Y Ȳ -1 2 .
Then, the probability distribution of Y is described by the following proposition.

Proposition 1. Let Y follow the Riemannian Gaussian density (3) with Ȳ = I N . If X = e N β σ 2 Y , then the probability distribution of X is given by

P(X ∈ B) ∝ B etr - log 2 (X) 2σ 2 {dX} (8)
for any measurable subset B of P β N . Here, etr(•) = exp(tr(•)) and the notation {dX} was introduced after (4).

In plain words, this proposition states that X follows a log-normal matrix ensemble. If X is diagonalised as

X = U xU † , then U is uniformly distributed on K β N
, and the eigenvalues (x 1 , . . . , x N ), which are all positive, have joint probability density function

p(x|σ) ∝ |V (x)| β i ρ(x i , 2σ 2 ) ( 9
)
where

V (x) = i<j (x j -x i ) is the Vandermonde determinant, and ρ(x, k) = exp(-log 2 (x)/k) is the log-normal weight function.
In essence, Proposition 1 is already contained in [START_REF] Santilli | Riemannian Gaussian distributions, random matrix ensembles and diffusion kernels[END_REF]. As a consequence of this proposition, the multiple integral z β (σ) of ( 7) may be expressed as follows

z β (σ) = 2πσ 2 N/2 × exp -N N 2 β (σ 2 /2) × 1 N ! R N + |V (x)| β ω(dx) (10)
where

ω(dx) = ρ(dx i ) (product over i = 1, . . . , N ) with ρ(dx i ) = 2πσ 2 -1/2 ρ(x i , 2σ 2 )dx i .
In [START_REF] Santilli | Riemannian Gaussian distributions, random matrix ensembles and diffusion kernels[END_REF], integrals as the one in [START_REF]Global and local scaling limits for the β= 2 stieltjes-wigert random matrix ensemble[END_REF] are expressed using the Andreev or De Bruijn identities, often employed in random matrix theory.

These yield somewhat cumbersome formulae, involving determinants of size N × N . On the other hand, the present paper focuses on an alternative approach, which turns out to be more suitable from a practical point of view. Instead of expressing z β (σ) exactly, by means of complicated formulae, the aim is to use a highly efficient approximation, which involves a single analytic expression. This is indicated by the following proposition. Proposition 2. In the limit where N → ∞ and σ → 0, while the product t = N σ 2 remains constant,

1 N 2 log z β (σ) -→ β 2 Φ β 2 t Φ(ξ) = ξ 6 - Li 3 (e -ξ ) -Li 3 (1) ξ 2 (11)
where Li 3 is the trilogarithm function, Li

3 (η) = ∞ k=1 η k /k 3 .
In [START_REF] Santilli | Riemannian Gaussian distributions, random matrix ensembles and diffusion kernels[END_REF], the limit in ( 11) is only mentioned in passing, and not stated under the same form. Here, this limit will be given centre stage. Proposition 2 states that ( 11) is valid in the "double-scaling regime" (N → ∞ and σ → 0), but numerical experiments have shown that log z β (σ) can be replaced by the expression afforded by [START_REF] Amari | Information geometry and its applications[END_REF] without notable loss of accuracy, whenever σ is small in comparison with N (this is further illustrated below).

At least informally, this can be justified by appealing to arguments originating in theoretical physics [START_REF] Mariño | Chern-Simons theory, matrix models, and topological strings[END_REF].

Considered as a function of t = N σ 2 , F (t) = log z β (σ) is called the Free energy (log of partition function).

This free energy can be expanded in an asymptotic series (see Section 1.3 of [START_REF] Mariño | Chern-Simons theory, matrix models, and topological strings[END_REF]),

F (t) ∼ ∞ g=0 F g (t) t N 2g-2 (12)
which is obtained by summing over Feynman diagrams, or so-called fatgraphs. Each coefficient

F g (t) is itself a series F g (t) ∼ h F g,h t h
, where F g,h counts fatgraphs which are said to have h holes and genus g (this means that one thinks of a fatgraph as a graph with h loops, drawn on a surface of genus g, such as a sphere or torus, etc). Now, accepting [START_REF] Pennec | Intrinsic statistics on Riemannian manifolds: Basic tools for geometric measurements[END_REF], it follows that for each fixed value of t,

1 N 2 F (t) = F 0 (t) + O 1 N 2 (13)
so that F 0 (t) approximates the left-hand side up to an error of the order of 1/N 2 . Finally, recalling that F (t) = log z β (σ), it is clear that F 0 (t) is the right-hand side of ( 11) -this follows by uniqueness of asymptotic expansions.

In addition to the asymptotic form of log z β (σ), another quantity of interest is the asymptotic empirical distribution of eigenvalues, of a random matrix Y which follows the Riemannian Gaussian density p(Y | Ȳ , σ) with Ȳ = I N .

Let (y 1 , . . . , y N ) denote the eigenvalues of Y , and consider their empirical distribution

νβ (I) = E |y i ∈ I| N ( 14
)
for any open interval I ⊂ R + , where E denotes expectation, and |y i ∈ I| the number of y i which belong to I.

In [START_REF] Santilli | Riemannian Gaussian distributions, random matrix ensembles and diffusion kernels[END_REF], the probability density function of ν2 was expressed as a weighted sum of Gaussian distributions, by a direct application of the Christoffel-Darboux formula, as in [START_REF] Deift | Orthogonal polynomials and random matrices : a Riemann-Hilber approach[END_REF]. It is possible to do so for any value of N , only because the β = 2 case can be studied using a well-known family of orthogonal polynomials, called Stieltjes-Wigert polynomials. No such analytic tool is available when β = 1 or 4.

To make up for this issue, the following proposition provides an asymptotic expression of the distribution νβ , valid for all values β = 1, 2, 4. Proposition 3. In the limit where N → ∞ and σ → 0, while the product t = N σ 2 remains constant, the empirical distribution νβ converges weakly to a distribution with probability density function n(y|βt/2), where

n(y|ξ) = 1 πξy arctan 4e ξ y -(y + 1) 2 y + 1 ( 15
)
on the interval a(ξ) ≤ y ≤ b(ξ), where a(ξ

) = c(1 + √ 1 -c) -2 and b(ξ) = c(1 - √ 1 -c) -2 , with c = e -ξ .
One hopes that, similar to the situation discussed after Proposition 2, the asymptotic density (15) approximates the finite-N empirical distribution νβ to such a good accuracy that one can replace νβ by this asymptotic density, for many practical purposes. This possibility will not be further investigated in the present paper.

IV. TOWARDS LEARNING APPLICATIONS

Riemannian Gaussian distributions were initially proposed as basic building blocks for learning models which aim to capture the intrinsic structure of statistical populations of covariance matrices. These include the mixture models and hidden Markov models, introduced in [2] [START_REF] Said | Hidden Markov chains and fields with observations in Riemannian manifolds[END_REF] and further developed in [START_REF] Tupker | Online learning of Riemannian hidden Markov models in homogeneous Hadamard spaces[END_REF]. In order to make use of these models in real-world applications, it is indispensable to know how to effectively compute the logarithm of the multiple integral z β (σ) of ( 7), for a dimension N which may be in the tens or hundreds.

Knowledge of log z β (σ) is already crucial in the simplest situation, where one tries to fit a single Riemannian Gaussian density p(Y | Ȳ , σ) (rather than a whole mixture) to data Y 1 , . . . , Y M ∈ P β N . Indeed, this will require setting Ȳ = ŶM and σ = σM , where ŶM and σM are the maximum-likelihood estimates, given in [START_REF] Said | Riemannian Gaussian distributions on the space of symmetric positive definite matrices[END_REF] [START_REF] Said | Gaussian distributions on Riemannian symmetric spaces : statistical learning with structured covariance matrices[END_REF]. Specifically,

ŶM = argmin Y ∈P β N M m=1 d 2 (Y m , Y ) (16)
is the Fréchet mean of the data Y m , with respect to the distance (2), and σM is the solution of the nonlinear equation

ϕ(σ M ) = 1 M M m=1 d 2 (Y m , ŶM ) ϕ(σ) = σ 3 d dσ log z β (σ) (17)
Therefore, it is already impossible to solve a toy problem, with a single Riemannian Gaussian density, without having some kind of hold on log z β (σ).

Until now (in [START_REF] Said | Riemannian Gaussian distributions on the space of symmetric positive definite matrices[END_REF] [3]), z β (σ) was approximated using an ad hoc Monte Carlo technique, which failed to produce coherent results for a dimension N just above N = 20. In the present section, the aim will be to show that significantly improved results can be obtained by solving Equation ( 17) after approximating log z β (σ) with the expression afforded by [START_REF] Amari | Information geometry and its applications[END_REF], according to Proposition 2.

First, numerical experiments were conducted to verify the validity of this new approximation. According to [START_REF] Pennec | Riemannian geometric statistics in medical image analysis[END_REF], the right-hand side of (11) should approximate log z β (σ)/N 2 up to an error of the order of 1/N 2 . To see that this is correct, the right-hand side of ( 11) was compared to certain exact expressions of log z β (σ). Namely, for the β = 2 case, one has the following expression, obtained in [START_REF] Said | Statistical models and probabilistic methods on Riemannian manifolds[END_REF],

log z 2 (σ) = N 2 log(2πσ 2 ) + N (N 2 -1)(σ 2 /6) + N -1 n=1 (N -n) log 1 -e -nσ 2 (18)
and for the β = 1 case, when the dimension N is even, the following expression, based on [START_REF] Santilli | Riemannian Gaussian distributions, random matrix ensembles and diffusion kernels[END_REF],

log z 1 (σ) = N 2 log(2πσ 2 ) -N (N + 1) 2 (σ 2 /8) + log Pf [M (σ)] (for N even) (19)
where Pf denotes the Pfaffian, equal to the square root of the determinant, and the matrix M (σ) has entries

M ij (σ) = exp (i 2 + j 2 )(σ 2 /2) erf((j -i)(σ/2)) i, j = 1, . . . , N (20)
with erf the error function.

Numerical evaluation of [START_REF] Said | Hidden Markov chains and fields with observations in Riemannian manifolds[END_REF] for large values of N or σ (up to σ = 10) is quite straightforward. Moreover, it immediately shows that ( 11) and ( 18) agree very closely when σ is smaller than N , and then gradually diverge away from one another as σ increases. Figure 1 provides graphical illustration for N = 10 and 20. Still larger values of N yield an even stronger match between [START_REF] Amari | Information geometry and its applications[END_REF] and [START_REF] Said | Hidden Markov chains and fields with observations in Riemannian manifolds[END_REF].

It is not equally straightforward to numerically evaluate [START_REF] Tupker | Online learning of Riemannian hidden Markov models in homogeneous Hadamard spaces[END_REF]. Even at N = 10, we begin encountering overflow problems for moderate values of σ when performing the computations in MATLAB. Still, as long as these overflow problems do not appear, it is possible to observe a close agreement between [START_REF] Amari | Information geometry and its applications[END_REF] and [START_REF] Tupker | Online learning of Riemannian hidden Markov models in homogeneous Hadamard spaces[END_REF], as in the β = 2 case. This is shown in Figure 2 for N = 6 and 12.

Based on the numerical results summarised in Figures 1 and2, it seems possible to use the right-hand side of [START_REF] Amari | Information geometry and its applications[END_REF], multiplied by N 2 , as a substitute for log z β (σ). While this is only an approximation, it is a highly efficient one, and has the advantage of being given by a rather simple analytic expression. In the β = 1 case, direct numerical evaluation of log z β (σ) is unstable for larger values of σ, and (11) offers a practical way out of this problem.

This can be verified by using symbolic computation (Mathematica), which allows us to extend the evaluation of [START_REF] Tupker | Online learning of Riemannian hidden Markov models in homogeneous Hadamard spaces[END_REF] to larger values of σ for larger N . It is nonetheless associated with the drawback that the resulting curves tend to be artificially non-smooth as in Figure 3a, due to numerical artifacts. These non-smooth features also exist in the curve depicted in Figure 3b, but do not appear visible at the given resolution. This behaviour is particularly problematic in the context of Equation ( 17), due to the presence of the derivative of log z β (σ). Furthermore, symbolic computation of ( 19) becomes exceedingly slow for sufficiently large N , and this can only be overcome by relying on approximations such as [START_REF] Amari | Information geometry and its applications[END_REF].

The second set of experiments directly addressed Equation [START_REF] Deift | Orthogonal polynomials and random matrices : a Riemann-Hilber approach[END_REF], in the β = 1 case, for a dimension N ranging between 10 and 25. This equation was solved using the Newton method, after its right-hand side was approximated according to [START_REF] Amari | Information geometry and its applications[END_REF]. Here, the solution obtained in this way will be denoted σM . This σM is an approximation of the maximum-likelihood estimate σM (the exact solution of ( 17)). If this approximation is any good, σM should approach the true value of σ for sufficiently large M (the number of data points Y m ). This can already be observed at M = 10 3 when N = 10, as reported in the following Table I. Each entry in this table gives the average value and standard deviation of σM , calculated over 20 independent trials (average ± standard deviation). It is clear that increasing M from 10 3 to 10 4 only reduces the standard deviation of σM , without really affecting its average.

When N = 10, the right-hand side of (17) can still be approximated using the Monte Carlo technique mentioned When N is above 20, Monte Carlo approximation of the right-hand side of ( 17) is not feasible anymore, and one is left only with the approximation using [START_REF] Amari | Information geometry and its applications[END_REF]. The solution σM obtained from this approximation is shown in Table III, for N = 20 and M = 10 4 . Here (contrary to Table II), σM systematically underestimates the true value of σ. In fact, identical behavior was observed for N between 20 and 25, along with similar values of σM .

A practical means of overcoming this issue would be to include a penalty term into Equation [START_REF] Deift | Orthogonal polynomials and random matrices : a Riemann-Hilber approach[END_REF], in order to enforce greater values of its solution. Then, the approximation using (11) can be successfully implemented for larger dimension N , where Monte Carlo approximation becomes useless (the present discussion stopped at N = 25, because the sampling algorithms used to generate the data points Y m could not be taken any further).

In conclusion, the present section has demonstrated the new approximation of log z β (σ), based on [START_REF] Amari | Information geometry and its applications[END_REF], significantly improves the existing Monte Carlo approximation, and also extends it to higher dimension. In upcoming work, the experiments conducted in this section will be generalised to more realistic learning models which may be applied to real-world data (for example, the mixture models considered in [START_REF] Zanini | Transfer learning: A Riemannian geometry framework with applications to brain-computer interfaces[END_REF]).

V. BLOCK-TOEPLITZ COVARIANCE MATRICES

Similar to the spaces P β N introduced in Section II, consider the spaces D β N , defined as follows.

• D 2 N is the space of N × N matrices with complex entries, whose operator norm is < 1.

• D 1
N is the space of N × N symmetric matrices with complex entries, whose operator norm is < 1.

Here, operator norm means the largest singular value. The space D 2 N will be called the Hermitian Siegel domain, and D 1 N the symmetric Siegel domain. The focus will be restricted to these β = 1 or 2 cases, since they are closely related to time-series analysis and signal processing (the β = 4 case is similar, but involves quaternion matrices).

Specifically, assume a wide-sense stationary N -variate time series of length T is described by its autocovariance matrix Γ. Wide-sense stationarity implies Γ has a block-Toeplitz structure, with T × T blocks of size N × N .

When solving an optimal prediction problem, one may apply a multidimensional Szegö-Levinson algorithm to the autocovariance Γ [21] [START_REF] Cabanes | Multidimensional complex stationary centred Gaussian autoregressive model classification : Applications for audio and radar machine learning in hyperbolic and Siegel spaces[END_REF], and obtain a family of matrices (Γ 0 , Ω 1 , . . . , Ω T -1 ), where Γ 0 is the zero-lag autocovariance of the original time series (this is a complex covariance matrix, Γ 0 ∈ P 2 N ), and each Ω t is a socalled matrix reflection coefficient, which belongs to D 2 N . If the autocovariance Γ has Toeplitz blocks, then each Ω t moreover belongs to D 1 N (which is a subspace of D 2 N ). Riemannian geometry enters this picture in the following way [START_REF] Jeuris | The Kähler mean of block Toeplitz matrices with Toeplitz structured blocks[END_REF] [START_REF] Cabanes | Multidimensional complex stationary centred Gaussian autoregressive model classification : Applications for audio and radar machine learning in hyperbolic and Siegel spaces[END_REF]. Consider an information metric on the space of block-Toeplitz autocovariance matrices Γ, equal to the Hessian of the entropy function S(Γ) = log det(Γ).

In terms of the new coordinates (Γ 0 , Ω 1 , . . . , Ω T -1 ), this information metric is a direct product of the affine-invariant metric (1) on the first coordinate Γ 0 , and of a scaled copy of the Siegel domain metric on each of the remaining coordinates Ω t (see formula (3.11) in [START_REF] Jeuris | The Kähler mean of block Toeplitz matrices with Toeplitz structured blocks[END_REF]). Dropping the subscript t, this Siegel domain metric is given by,

⟨u, v⟩ Ω = ℜ tr I N -ΩΩ † -1 u I N -ΩΩ † -1 v † u, v ∈ T Ω D β N ( 21
)
where the tangent space T Ω D β N is isomorphic to the space of N × N complex matrices if β = 2, and to the space of N × N symmetric complex matrices if β = 1 (this isomorphism shows that the dimension of D β N equals 2N N β ). The geodesic distance associated with the Siegel metric [START_REF] Jeuris | The Kähler mean of block Toeplitz matrices with Toeplitz structured blocks[END_REF] has the following expression [21] [22],

d 2 (Ξ, Ω) = tr arctanh 2 R 1 2 (Ξ, Ω) Ξ, Ω ∈ D β N ( 22
)
where

R(Ξ, Ω) = (Ξ -Ω)(I N -Ω † Ξ) -1 (Ξ † -Ω † )(I N -ΩΞ † ) -1 is called the matrix cross-ratio.
Riemannian Gaussian distributions on the space D β N are given by their probability density function p(Ω| Ω, σ) which is of the same form as in ( 3), but with the distance d(Ω, Ω) determined by [START_REF] Cabanes | Multidimensional complex stationary centred Gaussian autoregressive model classification : Applications for audio and radar machine learning in hyperbolic and Siegel spaces[END_REF], and the normalising factor

Z(σ) = D β N exp - d 2 (Ω, Ω) 2σ 2 dv(Ω) (23)
with respect to the Riemannian volume

dv(Ω) = det(I N -ΩΩ † ) -2N β {dΩ}. Here, {dΩ} = ij ℜ dΩ ij ℑ dΩ ij
where ℜ and ℑ denote the real and imaginary parts (the product is over i ≤ j if β = 1 and over all i, j if β = 2).

As shown in [START_REF] Said | Gaussian distributions on Riemannian symmetric spaces : statistical learning with structured covariance matrices[END_REF], a Riemannian Gaussian distribution on the space of block-Toeplitz covariance matrices Γ is just a product of independent Riemannian Gaussian distributions, one for each coordinate (Γ 0 , Ω 1 , . . . , Ω T -1 ).

For Γ 0 , this is a Riemannian Gaussian distribution on P 2 N , already considered in Section II. Thus, to understand Riemannian Gaussian distributions of block-Toeplitz covariance matrices, it only remains to study Riemannian Gaussian distributions on the Siegel domain D β N . In Section III, it was seen that Riemannian Gaussian distributions on P β N are equivalent to log-normal matrix ensembles. On the other hand, Riemannian Gaussian distributions on D β N are equivalent to "acosh-normal" ensembles, which will be described in Proposition 4 below.

Note first that each matrix Ω ∈ D β N can be factorised in the following way

Ω = U λV † (if β = 2) Ω = U λU T (if β = 1) (24)
where U and V are unitary, T denotes the transpose, and λ is diagonal, with diagonal elements λ i = tanh(r i ) for (r 1 , . . . , r N ) ∈ R N . The β = 2 case follows from the singular value decomposition of Ω, and the β = 1 case follows from the Takagi decomposition of Ω [START_REF] Edelman | Fifty three matrix factorizations: A systematic approach[END_REF]. In either case, λ i is of the form tanh(r i ) because the singular values of Ω are all < 1.

Proposition 4.

Let Ω follow a Riemannian Gaussian distribution on the Siegel domain D β N with Ω = 0 N (N × N zero matrix). If Ω is factorised as in [START_REF] Helgason | Differential geometry and symmetric spaces[END_REF], then U and V are uniformly distributed on the unitary group U (N ).

Moreover, if x i = cosh(2r i ) then (x 1 , . . . , x N ) have joint probability density function

p(x|σ) ∝ |V (x)| β i ρ(x i , 8σ 2 ) (x 1 , . . . , x N ) ∈ (1, ∞) N (25)
where V (x) = i<j (x j -x i ), and where ρ(x, k) = exp(-acosh 2 (x)/k) is the "acosh-normal" weight function.

Note that the assumption that the centre-of-mass parameter Ω is equal to 0 N does not entail any loss of generality.

The transformation [START_REF] Cabanes | Multidimensional complex stationary centred Gaussian autoregressive model classification : Applications for audio and radar machine learning in hyperbolic and Siegel spaces[END_REF]]

Ω -→ Ψ(Ω) = (I N -Ω Ω † ) -1 2 (Ω -Ω)(I N -Ω † Ω) -1 (I N -Ω † Ω) 1 2 (26)
maps Ω to 0 N , while preserving the Siegel domain metric [START_REF] Jeuris | The Kähler mean of block Toeplitz matrices with Toeplitz structured blocks[END_REF] and the associated distance and Riemannian volume.

Thus, if Ω follows a Riemannian Gaussian density p(Ω| Ω, σ), it is enough to replace Ω by Ψ(Ω), which will have a Riemannian Gaussian density with Ω = 0 N and with the same σ.

Proposition 4 implies that the transformed singular values x i = cosh(2r i ) follow the classical eigenvalue distribution of an orthogonal (β = 1) or unitary (β = 2) matrix ensemble with "acosh-normal" weight function.

Therefore, in particular, the normalising factor Z(σ) of ( 23) reduces to a multiple integral (compare to ( 6) and ( 10))

z β (σ) = 1 N ! (1,∞) N |V (x)| β ω(dx) (27)
where ω(dx) = ρ(dx i ), with ρ(dx i ) = ρ(x i , 8σ 2 )dx i for the weight function ρ(x, k) = exp(-acosh 2 (x)/k).

Thus, Riemannian Gaussian distributions on the Siegel domain D β N can be studied by applying the tools of random matrix theory to new "acosh-normal" matrix ensembles. Hopefully, in the near future, this will lead to similar results to the ones obtained with log-normal ensembles in Section III, paving the way to learning from datesets of high-dimensional block-Toeplitz covariance matrices.

P(Y ∈ B) = (Z(σ)) -1 B exp - d 2 (Y , I N ) 2σ 2 dv(Y ) (28)
Let (y 1 , . . . , y N ) denote the eigenvalues of Y . Using (2) and the fact that dv(Y ) = det(Y ) -N β {dY }, (28) becomes

P(Y ∈ B) = (Z(σ)) -1 B N i=1 exp - log 2 (y i ) 2σ 2 y -N β i {dY } (29)
Recall that X = e N β σ 2 Y . Accordingly, if (x 1 , . . . , x N ) are the eigenvalues of X, an elementary calculation yields

exp - log 2 (y i ) 2σ 2 y -N β i = exp[N 2 β (σ 2 /2)] × exp - log 2 (x i ) 2σ 2
Therefore, (29) can be written

P(Y ∈ B) = exp[N N 2 β (σ 2 /2)] × (Z(σ)) -1 B etr - log 2 (X) 2σ 2 {dY } (30)
To conclude, it is enough to use once more the definition X = e N β σ 2 Y . This implies

{dY } = exp[-2N N 2 β (σ 2 /2)]{dX}
Thus, using the fact that

P(X ∈ B) = P e N β σ 2 Y ∈ B
and changing the variable of integration from Y to X in (30), it follows that

P(X ∈ B) = exp[-N N 2 β (σ 2 /2)] × (Z(σ)) -1 B etr - log 2 (X) 2σ 2 {dX}
as required in [START_REF] Santilli | Riemannian Gaussian distributions, random matrix ensembles and diffusion kernels[END_REF].

B. Proof of Proposition 2

The β = 2 case follows from [START_REF] Said | Hidden Markov chains and fields with observations in Riemannian manifolds[END_REF], by an elementary calculation, after noting that

1 N 2 N -1 n=1 (N -n) log 1 -e -nσ 2
is a Riemann sum for the improper integral

1 0 (1 -x) log 1 -e -tx dx = Li 3 (1) -Li 3 (e -t ) t 2
For other values of β, the result will be obtained from the reasoning presented in the proof of Proposition 3, based on the scaling equation (33).

C. Proof of Proposition 3

The β = 2 case has already been proved in [START_REF] Said | Statistical models and probabilistic methods on Riemannian manifolds[END_REF]. To deal with the general case, write [START_REF]Global and local scaling limits for the β= 2 stieltjes-wigert random matrix ensemble[END_REF] under the form

z β (σ) = c β (σ) I β (σ), where c β (σ) = 1 N ! exp -N N 2 β (σ 2 /2)
and I β (σ) is the multiple integral

I β (σ) = R N + exp   - 1 2σ 2 N i=1 log 2 (x i ) + β i<j log |x j -x i |   dx Now, if µ N = (1/N ) N i=1
δ xi is the empirical distribution of x i (where δ xi denotes the Dirac measure at x i ), then

I β (σ) = R N + exp -N 2 E β (µ N , t) dx (31)
where E β (µ, t) is the so-called energy functional

E β (µ, t) = 1 2N t R+ log 2 (x)µ(dx) -β R+ R+ log |y -x|µ(dx)µ(dy) (32)
defined for any probability distribution µ on R + . This energy functional satisfies the following scaling equation

E β (µ, t) = (β/2)E 2 (µ, (β/2)t) (33)
which is easily obtained after dividing (32) by β/2. The proof can now be completed by applying to (31) and (33) the arguments in [START_REF] Deift | Orthogonal polynomials and random matrices : a Riemann-Hilber approach[END_REF]. First [START_REF] Deift | Orthogonal polynomials and random matrices : a Riemann-Hilber approach[END_REF] (Corollary 6.90, Page 155),

1 N 2 log I β (σ) -→ -E β (µ * , t)
where µ * = µ * β (t) is the so-called equilibrium distribution, the unique minimiser of the energy (32) among probability distributions on R + . From the scaling equation (33), it is now clear

1 N 2 log I β (σ) -→ -(β/2)E 2 (µ * , (β/2)t) (34)
Therefore, by adding the limit of log c β (σ)/N 2 ,

1 N 2 log z β (σ) -→ -(β/2) ((βt/4) + E 2 (µ * , (β/2)t)) (35)
However, when β = 2, it is already known that this limit is Φ(t), where Φ is defined in [START_REF] Amari | Information geometry and its applications[END_REF]. This provides E 2 (µ * , t) = -t/2 -Φ(t), which can be replaced back into (35), yielding the general case of Proposition 2.

To prove Proposition 3, it is possible to use the scaling equation (33), once more. From [START_REF] Deift | Orthogonal polynomials and random matrices : a Riemann-Hilber approach[END_REF] (Section 6.4), if μN is defined as in ( 14), but with the x i instead of the y i , then μN converges weakly to the equilibrium distribution µ * = µ * β (t). On the other hand, the scaling equation (33) implies µ * = µ * 2 (βt/2), and µ 2 (t) is already known to be the image of the distribution with density n(y|t), defined as in [START_REF] Meckes | The random matrix theory of the classical compact groups[END_REF], under the change of variables x = e β 2 t y. This provides µ * for any value of β, and the Proposition then follows by changing the variables back from x to y.

D. Proof of Proposition 4

The proof relies on the general theory of Gaussian distributions on Riemannian symmetric spaces, as outlined in [3] [20]. Here, to keep the proof self-contained, it will be helpful to briefly recall certain aspects of this theory (for a more detailed, in-depth discussion, the reader is referred to [START_REF] Said | Statistical models and probabilistic methods on Riemannian manifolds[END_REF], Sections 1.9 and 3.3).

A Gaussian distribution on a Riemannian symmetric space M is defined by its probability density function

p(x|x, σ) = (Z(σ)) -1 exp - d 2 (x, x) 2σ 2 x ∈ M (36)
with respect to the Riemannian volume dv(x) on M , with x ∈ M and σ > 0.

It is always assumed M is a Riemannian symmetric space of non-positive curvature, associated to a symmetric pair (G, K). Precisely, G is a connected Lie group which acts transitively and isometrically on M (this action is denoted x → g • x, for g ∈ G and x ∈ M ), and K is a compact subgroup of G, made up of those elements k ∈ G

which fix a certain point o ∈ M (that is, k ∈ K if and only if k • o = o).
For a concrete understanding of Gaussian distributions on M , it is necessary to take a closer look at the Lie algebras of G and K, denoted g and k, respectively. These are related together by the so-called Cartan decomposition, g = k + p (direct sum), where the subspace p of g can be identified with the tangent space T o M (tangent space to M at at o). In terms of this decomposition, the main construction needed for the proof can be described as follows.

Without any loss of generality, G and K are taken to be matrix Lie groups, and g and k their matrix Lie algebras.

Let a be a maximal abelian subspace of p (that is, all the matrices a ∈ a commute with each other). Then, each matrix u ∈ p can be written under the form u = kak -1 , where k ∈ K and a ∈ a (this factorisation is the general template for the fifty three matrix factorisations, outlined in [START_REF] Edelman | Fifty three matrix factorizations: A systematic approach[END_REF]). Moreover, any x ∈ M admits a representation

x = exp kak -1 • o k ∈ K and a ∈ a (37)
where exp denotes the matrix exponential. Incidentally, this representation is not unique, but for almost all x ∈ M (for all x ∈ M , except a subset of zero volume) x has exactly |W | couples (k, a) which satisfy (37), where |W | is the number of elements of the Weyl group of the symmetric couple (G, K).

It is now possible to state the following Lemma 1 [START_REF] Said | Gaussian distributions on Riemannian symmetric spaces : statistical learning with structured covariance matrices[END_REF] [20], which will yield the entire proof, by direct application.

Lemma 1. Let x follow the Gaussian density (36) on M , with x = o. If x is represented as in (37), then k is uniformly distributed on the compact group K. Moreover, a has the following probability density function on a (recall that a is a real vector space, so the density is with respect to the usual Lebesgue measure on a)

p(a|σ) ∝ exp - ∥a∥ 2 o 2σ 2 ρ∈∆+ (sinh |ρ(a)|) mρ (38)
where ∥a∥ o is the Riemannian norm of a (since a ∈ p, it can be identified with a vector in T o M ), and ∆ + is a set of positive roots ρ on a, with respective multiplicities m ρ (each ρ ∈ ∆ + is a certain linear function ρ : a → R).

Recall that positive roots ρ : a → R are any set of linear functions on a such that, for any a ∈ a, the eigenvalues of the linear operator ad a : g → g, given by ad a (ξ) = [a, ξ] (this is aξ -ξa), are equal to ±ρ(a) with respective multiplicities m ρ .

The proof may now begin in earnest. It merely consists in identifying G, K, g, k, p, and a, for each Siegel domain D β N , and then writing down the corresponding version of Lemma 1, which directly yields Proposition 4. Fortunately, all the necessary information can be found in [START_REF] Piatetskii-Shapiro | Automorphic functions and the geometry of classical domains[END_REF] [START_REF] Helgason | Differential geometry and symmetric spaces[END_REF].

The symmetric pair : D β N is associated to the symmetric pair (G, K), where G and K are groups of 2N × 2N complex matrices g, defined in the following way (here, * denotes the complex conjugate). ξ = (ω ij , τ ij , ω ij ) =⇒ ad a (ξ) = (r i -r j )ξ ξ = i (τ ij , ω ij , τ ij ) =⇒ ad a (ξ) = (r i -r j )ξ ξ = i (-τ ij , τ ij , τ ij ) =⇒ ad a (ξ) = (r i + r j )ξ ξ = (ω ji , ω ij , ω ij ) =⇒ ad a (ξ) = (r i + r j )ξ

where τ ij = e ij + e ji and ω ij = e ij -e ji , with e ij a matrix all of whose entries are zero, except the one on line i and column j, which is equal to 1, and where r = diag(r 1 , . . . , r N ).

This shows that the positive roots are ρ(a) = r i -r j where i < j and ρ(a) = r i + r j where i ≤ j, which all have multiplicity m ρ = 2 (this is β), except for ρ(a) = 2r i (this is r i + r j when i = j), which has multiplicity m ρ = 1. Moreover, of the above eigenvectors, only the ones in the left column belong to g β=1 . Thus, in the β = 1 case, all the multiplicities m ρ are equal to 1. Finally, note that, from the power series of the matrix exponential,

From the above form of the action of G on D β N , it is then straightforward that the representation (37) is the same as [START_REF] Helgason | Differential geometry and symmetric spaces[END_REF]. This implies the first part of Proposition 4 (U, V are uniformly distributed on U (N )). For the second part, it follows from (38) that, in the present case, p(r|σ) ∝

2 Fig. 1 : 2 Fig. 2 :

 2122 Fig. 1: Proposition 2 (β = 2) : (11) in dashed red and (18) in solid blue

10 Fig. 3 :

 103 Fig. 3: Proposition 2 (β = 1) : (11) in dashed red and (19) in solid blue for N = 50

  APPENDIX A PROOFS OF PROPOSITIONS 1 TO 4 A. Proof of Propositon 1 If Y follows the Riemannian Gaussian density (3) with Ȳ = I N , then



  : gPg † = P and gSg T = S These groups G and K act on D β N by matrix fractional transformationsg • Ω = (AΩ + B)(CΩ + D) -1 g ∈ G and Ω ∈ D β NThe Lie algebras : these are given byg β=2 = {(γ, δ, ε) : γ + γ † = ε + ε † = 0 N } ; g β=1 = {(γ, δ, ε) : γ + γ † = ε + ε † = δ -δ T = γ + ε T = 0 N } k β=2 = {(γ, 0 N , ε) : γ + γ † = ε + ε † = 0 N } ; k β=1 = {(γ, 0 N , γ * ) : γ + γ † = 0 N } in terms of the notation (γ, δ, ε) =The subspaces p and a : these are given byp β=2 = {(0 N , δ, 0 N )} ; p β=1 = {(0 N , δ, 0 N ) : δ -δ T = 0 N } a β=2 = a β=1 = {a = (0 N , r, 0 N ) : r real diagonal}Lemma 1 : the only information still needed is the positive roots. Matrix multiplication shows that for any a ∈ a, ad a has the following eigenvalues and eigenvectors in g β=2 .

  0 N , r, 0 N )

i exp - r 2 i 2σ 2

 22 sinh(2r i ) i<j (sinh |r i -r j | sinh |r i + r j |) β (40)This yields (25), by using 2 sinh |r i -r j | sinh |r i + r j | = cosh(2r i ) -cosh(2r j ), and introducing the change of variables x i = cosh(2r i ).

  [START_REF] Amari | Information geometry and its applications[END_REF] in dashed red and[START_REF] Tupker | Online learning of Riemannian hidden Markov models in homogeneous Hadamard spaces[END_REF] in solid blue for N = 50

	true σ	1	2	3	4	5	6	7
	M = 10 3	1.09 ± 0.00 2.01 ± 0.01 3.11 ± 0.30 3.56 ± 0.60 4.58 ± 0.60 5.13 ± 0.50 5.60 ± 0.75
	M = 10 4	1.09 ± 0.00 2.01 ± 0.00 3.11 ± 0.15 3.80 ± 0.16 4.50 ± 0.23 5.37 ± 0.31 5.60 ± 0.37

TABLE I :

 I

	true σ	1	2	3	4	5	6	7
	M = 10

The solution σM of

[START_REF] Deift | Orthogonal polynomials and random matrices : a Riemann-Hilber approach[END_REF]

, for β = 1 and N = 10 (r.h.s. approximated using

[START_REF] Amari | Information geometry and its applications[END_REF]

) in

[START_REF] Said | Riemannian Gaussian distributions on the space of symmetric positive definite matrices[END_REF]

[START_REF] Said | Gaussian distributions on Riemannian symmetric spaces : statistical learning with structured covariance matrices[END_REF]

. The solution obtained with this approximation will be denoted [σ M ] MC . The following Table

II

shows that the estimation error from σM is quite improved, in comparison with the estimation error from [σ M ] MC . Moreover, [σ M ] MC seems to systematically overestimate the true value of σ.

3

1.16 ± 0.01 2.90 ± 0.01 5.63 ± 0.93 7.38 ± 1.35 8.35 ± 1.50 9.40 ± 1.11 9.95 ± 1.12 M = 10 4 1.16 ± 0.01 2.90 ± 0.00 5.40 ± 0.37 7.00 ± 0.60 9.00 ± 0.50 9.50 ± 0.31 9.77 ± 0.35

TABLE II :

 II

The solution [σ M ] MC of (

17

), for β = 1 and N = 10 (MC approximation of r.h.s.)

TABLE III :

 III The solution σM of[START_REF] Deift | Orthogonal polynomials and random matrices : a Riemann-Hilber approach[END_REF], for β = 1 and N = 20 (r.h.s. approximated using (11))

ACKNOWLEDGMENT S.H. is supported by the Science and Technology Facilities Council (STFC) and St. John's College, Cambridge. S.H. also benefited from partial support from the Cambridge Mathematics Placement (CMP) Programme and the European Research Council under the Advanced ERC Grant Agreement Switchlet n.670645. C.M. was supported by a Henslow Fellowship from the Cambridge Philosophical Society as well as a Presidential Postdoctoral Fellowship at NTU and is grateful for support from Fitzwilliam College, Cambridge.


Michaël Goujon
email: michael.goujon@udamail.fr.

Aristide Mabali

Cerdi

THE DIFFERENT IMPACTS OF DIFFERENT TYPES OF NATURAL RESOURCES ON POLITICAL INSTITUTIONS IN DEVELOPING COUNTRIES

Keywords: Natural Resources, Rents, Political Institutions, Panel Data Models JEL classification: P48, Q34, C23

Rents generated by natural resources are usually thought to weaken the quality of institutions, particularly in developing countries. Our hypothesis is that this effect may differ depending on the types of natural resources characterized by their different degree of appropriability. We test this hypothesis using panel data covering 90 developing countries for the period 1970-2010. We find that total rents weaken the quality of institutions. However, while oil rents have a significant negative effect, forest and mineral rents do not, after controlling for the other relevant determinants of institutional quality, institutional persistence, neighbor effect, and endogeneity of rents.

INTRODUCTION

There is now a large consensus that the quality of institutions is a pillar of economic development, alongside other important determinants such as geography, human capital, and openness [START_REF] Knack | Institutions and Economic Performance: Crosscountry Tests Using Alternative Institutional Measures[END_REF][START_REF] Mauro | Corruption and Growth[END_REF][START_REF] Porta | The Quality of Government[END_REF][START_REF] Acemoglu | Institutions as the Fundamental Cause of Long-run Growth[END_REF][START_REF] Rodrik | Growth Strategies[END_REF][START_REF] Kanyama | Quality of Institutions: Does Intelligence Matter?[END_REF]. 1 However, the determinants of the quality of institutions are still being debated. [START_REF] Acemoglu | Institutions as the Fundamental Cause of Long-run Growth[END_REF] argue that groups which win political power choose institutions that serve their interest, and that political power mainly depends on the distribution -or capture -of economic resources. 2 In this paper we focus on natural resource endowment as a strong structural characteristic of economic resource distribution and institutional change in developing countries.

The literature broadly supports the idea that rents from natural resources are harmful to institutional quality in developing countries (see [START_REF] Frankel | The Natural Resource Curse: A Survey[END_REF][START_REF] Van Der Ploeg | Natural resources: Curse or blessing?[END_REF] for recent reviews). Firstly, appropriation disputes explain why natural resources can encourage rent-seeking, corruption, or conflicts [START_REF] Boschini | Resource Curse or Not: A Question of Appropriability[END_REF]. Secondly, governments can also suffer from the so-called "rent effect" of natural resources through (i) taxation, (ii) spending, and (iii) group formation [START_REF] Ross | Does Oil Hinder Democracy?[END_REF][START_REF] Omgba | On the Duration of Political Power in Africa: The Role of Oil Rents[END_REF]. The taxation effect suggests that when governments derive sufficient revenues from natural resources, they tax their population less heavily. In turn, the population is less likely to demand accountability from, and representation in, the government, following the so-called principle of "no representation without taxation" [START_REF] Tilly | The Formation of National States in Western Europe[END_REF]. The spending effect suggests that rents lead to patronage, which dampens latent pressures for democratization [START_REF] Atkinson | Savings, Growth and the Resource Curse Hypothesis[END_REF][START_REF] Vicente | Does oil corrupt? Evidence from a natural experiment in West Africa[END_REF]. The group formation effect suggests that the government will use its largesse to prevent the formation of contestant groups [START_REF] Mahdavy | The Patterns and Problems of Economic Development in Rentier States: The Case of Iran[END_REF][START_REF] Anderson | The State in the Middle East and North Africa[END_REF][START_REF] Ross | Does Oil Hinder Democracy?[END_REF][START_REF] Andersen | Constitutions and the Resource Curse[END_REF].

However, the concept of "rentier states" suggests that different types of natural resources could have different impacts on the quality of institutions. The idea is that natural resources differ in their ease of appropriation (or capturability of the rent), inducing a more or less unequal distribution of the rent. Thus, [START_REF] Bulte | Resource Intensity, Institutions, and Development[END_REF] distinguish between "point resources" and "diffuse resources". Point resources generate rents that can be more easily captured, and so are more harmful to the quality of institutions than diffuse resources. Similarly, [START_REF] Boschini | Resource Curse or Not: A Question of Appropriability[END_REF] focus on the technical appropriability of a resource. They argue that resources such as diamonds or precious metals are easily transportable, storable and saleable, and are therefore more detrimental than agricultural products. [START_REF] Petermann | Mining and Corruption[END_REF] consider that rents from oil are more harmful than those derived from forests, because oil involves contracts of higher value and is more technology-intensive. [START_REF] Engerman | Factor Endowments, Institutions, and Differential Paths of Growth Among New World Economies: A View From Economic Historians of the United States[END_REF] and [START_REF] Ross | Does Oil Hinder Democracy?[END_REF] pointed out that oil rents are more likely to degrade institutional quality because they are easily captured by the central government (see also [START_REF] Robinson | Political Foundations of the Resource Curse: A Simplification and a Comment[END_REF]. Conversely, mineral or agricultural rents may contribute to improving the quality of institutions because they can be captured by a larger share of the population.

We investigate whether rents generated by natural resources have a detrimental effect on institutions, by considering different types of natural resources characterized by their different degree of appropriability. Taking advantage of a large dataset covering developing countries, we use appropriate techniques for panel data (Fixed Effects model, GMM-System and 2SLS estimators) in order to take into account unobserved heterogeneities, endogeneity of rents, and institutional persistence. Panel data offer a solution to the endogeneity problem through the use of lagged values as instruments for institutions" where historical events determined institutions that have persisted over time.

endogenous variables [START_REF] Lederman | Natural Resources Neither Curse nor Destiny[END_REF].

This paper is organized as follows. Section 2 presents the empirical model and data. Section 3 discusses the results of model estimations. Section 4 presents robustness checks, and Section 5 concludes by offering policy discussions.

EMPIRICAL STRATEGY AND DATA

The few previous studies which have focused on the direct relationship between natural resources and institutions suffer from a number of limitations 3 : (i) the distinction between point resources and diffuse resources is not done or done imperfectly; (ii) the measures for rents, usually based on the ratio of primary commodity exports to GDP, are imperfect; (iii) they use cross-sectional data that do not allow to properly take into account unobserved heterogeneities [START_REF] Lederman | Natural Resources Neither Curse nor Destiny[END_REF][START_REF] Borge | Local Natural Resource Curse?[END_REF]. In this paper we attempt to resolve these shortcomings by using appropriate econometric techniques and measures for rents.

We use annual data covering the period from 1970 to 2010, for 90 developing countries (see country list in Table A1 in Appendix). Our baseline specification is as follows:

= + + + + ∑ + , (1)
where is a measure of institutional quality in country at time . Inertia in institutions is modeled by lagged among the regressors. is a measure of natural resource rents (total, oil, mineral, and forest), is a vector of control variables, is a country dummy variable for country fixed effects, is a time dummy variable which controls for common shocks (e.g. international context), and is an idiosyncratic errors term. The parameter of interest is and is expected to be negative. is lagged by one period to model a delayed effect of a change in rents on institutions [START_REF] Bhattacharyya | Natural Resources, Democracy and Corruption[END_REF][START_REF] Brückner | Rain and the Democratic Window of Opportunity[END_REF].

However, the estimation of Equation 1 faces two main issues. Firstly the presence of the lagged dependent variable among the regressors. We take this into account by using the Generalized Method of Moments (GMM) procedure. Secondly, the probable endogeneity of Rent. Endogeneity might come from measurement errors in Rent series 3 See for instance [START_REF] Bulte | Resource Intensity, Institutions, and Development[END_REF][START_REF] Brunnschweiler | Cursing the Blessings? Natural Resource Abundance, Institutions, and Economic Growth[END_REF]Brunnsschweiler and Bulte, 2008. Other empirical works explore institutions as a transmission channel or a state variable conditioning the link between natural resources and economic performances [START_REF] Mehlum | Institutions and the Resource Curse[END_REF][START_REF] Brunnschweiler | Cursing the Blessings? Natural Resource Abundance, Institutions, and Economic Growth[END_REF][START_REF] Bulte | Resource Intensity, Institutions, and Development[END_REF][START_REF] Boschini | Resource Curse or Not: A Question of Appropriability[END_REF], or rent-seeking activities [START_REF] Torvik | Natural Resources, Rent Seeking and Welfare[END_REF][START_REF] Bhattacharyya | Natural Resources, Democracy and Corruption[END_REF]. In some ways the "political survival models" in which resource rents affect the duration in office of political leaders are closer to our theme [START_REF] Caselli | Power Struggles and the Natural Resource Curse[END_REF][START_REF] Omgba | On the Duration of Political Power in Africa: The Role of Oil Rents[END_REF][START_REF] Cuaresma | Oil and the duration of dictatorships[END_REF][START_REF] Andersen | Constitutions and the Resource Curse[END_REF].

(attenuation bias). These errors are also likely to be correlated with (attrition bias). Endogeneity might also come from the omission of significant factors in the vector . Following [START_REF] Carmignani | The Social Development Effects of Primary Commodity Export Dependence[END_REF], we address endogeneity by using both the instrumental variables estimator, with the hypothesis = 0, and dynamic panel estimates, relaxing the restriction, = 0, which allows us to take into account the persistence in institutions. In the dynamic panel model we use five year average data to eliminate short-term variability. According to [START_REF] Roodman | How to Do Xtabond2: an Introduction to "Difference" and "System" GMM in Stata[END_REF], the GMM-System estimator is only relevant if the cross-country dimension is greater than the time dimension, which is the case when using five-year average data. According to the World Bank's definition, "natural resources rents are calculated as the difference between the price of a commodity and the average cost of producing it. This is done by estimating the world price of units of specific commodities and subtracting estimates of average unit costs of extraction or harvesting costs (including a normal return on capital). These unit rents are then multiplied by the physical quantities countries extract or harvest to determine the rents for each commodity as a share of gross domestic product (GDP)". 4 We use this measure of rents for the following reasons: (i) it is a good proxy for revenues from natural resources that can be captured by the State, politicians, or social groups, (ii) the time and country coverage is broad, thus reducing selection bias, (iii) the use of exogenous world prices limits the endogeneity problem, (iv) this measure is widely used in the literature (e.g. Ross, 2006;[START_REF] Collier | Testing the Neocon Agenda: Democracy in Resource-rich Societies[END_REF][START_REF] Bhattacharyya | Natural Resources, Democracy and Corruption[END_REF]. To test the different impact of different types of natural resources on institutional quality, we use data series for total rents, and oil, mineral and forest rents. Unfortunately, data on other raw materials rents are very limited. Figure 1 shows a scatter plots with rents-to-GDP ratio (and different types of rents-to-GDP ratio) on the X-axis and polity2 score on the Y-axis, and Figure 2 shows the evolution of total rents over the period 1970-2010 for our sample of countries.

Figure 2. Average Rents and Polity2 Over Time, 90 Developing Countries

For the quality of institutions variable, we use the democracy score polity2 from the polityIV database (Marshall and Jaggers, 2002). The choice of this indicator was motivated by the following reasons: (i) it captures both de facto and de jure political power as discussed by [START_REF] Acemoglu | Institutions as the Fundamental Cause of Long-run Growth[END_REF] 5 , (ii) it covers a large number of countries 4 For more details, refer to [START_REF] Hamilton | Genuine Savings Rates in Developing Countries[END_REF], [START_REF] Collier | Testing the Neocon Agenda: Democracy in Resource-rich Societies[END_REF], or [START_REF] Bhattacharyya | Natural Resources, Democracy and Corruption[END_REF]. Sources and methods are described in [START_REF] Jarvis | The Changing Wealth of Nations: Measuring Sustainable Development in the New Millennium[END_REF].

5 [START_REF] Acemoglu | Institutions as the Fundamental Cause of Long-run Growth[END_REF] distinguish two types of political power. The first one is de jure political power whose legitimacy comes from political institutions (e.g. constitution). The second is de facto political power (Glaeser et al., 2004, Persson and[START_REF] Persson | Democracy and Development. The Devil in Detail[END_REF]Ross, 2006;[START_REF] Collier | Testing the Neocon Agenda: Democracy in Resource-rich Societies[END_REF][START_REF] Bhattacharyya | Natural Resources, Democracy and Corruption[END_REF]. The polity2 score varies between -10 and +10, and increases with the quality of institutions. Figure 1 shows a high heterogeneity of institutional quality between countries. The average polity2 score, reported in Figure 2, increases over the period 1970-2010 with an acceleration at the beginning of the 1990s. This is consistent with global political events at that time, marked by the wave of democratization which followed the fall of the Berlin Wall in 1989, and in Africa, the La Baule Conference in 1990. Below we rescale the polity2 variable so that it ranges from 0 to 1, with higher values indicating better democratic institutions. We also use the Political Rights index (PR) from Freedom House as an alternative measure of institutional quality for robustness tests. We use a set of control variables commonly used in the literature concerning the factors affecting institutions (see Table A2 in appendix for data description): inflation that mainly relies on the distribution of economic resources between social groups, and which determines their ability to use and modify political institutions, even though social conflict. The group that exercises de facto political power is tempted to establish new political institutions in order to legitimize its power and maintain the capture of resources, thus introducing changes in institutions. [START_REF] Haggard | The Political Economy of Democratic Transitions[END_REF]; official development assistance [START_REF] Svensson | Aid, Growth and Democracy[END_REF][START_REF] Collier | Aid, Policy and Growth in Post-conflict Societies[END_REF]; remittances [START_REF] Abdih | Remittances and Institutions: Are Remittances a Curse?[END_REF]; GDP per capita (Barro, 1999;[START_REF] Przeworski | Democracy and Development[END_REF][START_REF] Epstein | Democratic Transitions[END_REF][START_REF] Brückner | Rain and the Democratic Window of Opportunity[END_REF]; educationmeasured by the gross enrolment rate [START_REF] Leite | Does Mother Nature Corrupt? Natural Resources, Corruption and Economic Growth[END_REF][START_REF] Kanyama | Quality of Institutions: Does Intelligence Matter?[END_REF]; openness [START_REF] Leite | Does Mother Nature Corrupt? Natural Resources, Corruption and Economic Growth[END_REF][START_REF] Brunnschweiler | The Resource Curse Revisited and Revised: A Tale of Paradoxes and Red Herrings[END_REF]. Geography and ethnic fragmentation, which are constant over time and country specific, are captured using country fixed effects.

We add a proxy for the institutional quality in neighbor countries which aims to capture a diffusion effect of institutions at the regional level [START_REF] Barro | Economic Growth in a Cross-Section of Countries[END_REF][START_REF] Ades | Thy Neighbor's Curse: Regional Instability and Economic Growth[END_REF][START_REF] Kelejian | Spatial Spillovers in the Development of Institutions[END_REF]. This variable is computed as the mean of the polity2 scores of contiguous countries. The variance-covariance matrix in Table 1 shows that there is no obvious collinearity between the control variables.

EMPIRICAL RESULTS

The Effect of Total Rents on Institutional Quality

Table 2 reports the results of the estimates of Equation 1. 6 In Column 1, the OLS estimate shows a statistically significant negative effect of total rents on the polity2 score, suggesting that natural resources are associated with a poor quality of institutions. But this association may be driven by unobserved country characteristics. To tackle this issue, we use the fixed effects estimator in Columns 2 and 3, which controls for country characteristics that do not vary over time, and we add time varying common shocks in Column 3. The negative relationship remains but the magnitude of the impact coefficient falls by half, suggesting that there is a significant impact of country fixed effects on the quality of institutions. In Column 4 we use the Between estimator that uses only inter-country variability. The estimated coefficient for Rent is then statistically significant, and is higher than the one obtained with the fixed effects estimator. This seems to corroborate the results of [START_REF] Bhattacharyya | Natural Resources, Democracy and Corruption[END_REF] who showed that inter-country variability of institutions outweighs intra-country variability, due to the inertia of institutions.

Columns 5 to 7 address the endogeneity of Rent and the persistence of institutions. The 2SLS estimator in Column 5 is performed using two instrumental variables, which are expected to be exogenous and are expected to be significant factors in the exploitation of natural resources without affecting institutions directly: the ratio of international oil price to population size, and population density. The international oil price is an instrument for rents in [START_REF] Collier | Testing the Neocon Agenda: Democracy in Resource-rich Societies[END_REF] and [START_REF] Bhattacharyya | Natural Resources, Democracy and Corruption[END_REF] considering the oil price to be exogenous and correlated to the price of other natural resources. We divide this variable by the population size to generate country variability. Results from the first stage equation are reported in Table 3. The validity of the instruments is verified using three statistics: (i) The Kleibergen-Paap rk F-statistic which is well above the critical values compiled by [START_REF] Stock | Testing for Weak Instruments in Linear IV Regression[END_REF], suggesting that the instruments are not weak, (ii) our model is identified because the p-value associated with the over-identification test is zero, (iii) the J-test of over-identifying restrictions has a p-value above 0.1 supporting the choice of instruments because they are not correlated with the error terms (Van der Ploeg and Poelhekke, 2010). In column (6), even if [START_REF] Stock | Testing for Weak Instruments in Linear IV Regression[END_REF]'s test previously rejected the instruments' weakness, we also perform the LIML estimator.

The Generalized Method of Moments (GMM) procedure addresses two concerns: (i) the presence of the lagged dependent variable among the regressors, and (ii) the correlation between fixed effects and explanatory variables. We use the two-step [START_REF] Arellano | Some Tests of Specification for Panel Data: Monte Carlo Evidence and an Application to Employment Equations[END_REF] estimator, where the first-differenced lagged dependent variable is also instrumented with its past levels [START_REF] Roodman | How to Do Xtabond2: an Introduction to "Difference" and "System" GMM in Stata[END_REF]. Here, we take care to correct heteroscedasticity (Windjmeier, 2005). The validity of the GMM estimates is verified since there is no first-order or second-order serial correlation in the residuals in the first-differences (AR (1) and AR (2) tests) and the instruments are exogenous (Hansen J-statistic).

In Columns 5 to 7, the instrumental variable estimators generate a higher score for the Rent coefficient than the other estimators, confirming the attenuation bias. A doubling (an increase of 100%) of total rents is associated with a significant decrease of 0.28 points in the polity2 score, as shown in Columns 5 and 6. The GMM-System results reveal the inertia in institutions. The size of the Rent coefficient is lower with the GMM-System than with the 2SLS estimator but remains significant, showing that results are robust to changes in the estimators and data frequency.

With regard to the other regressors, trade openness, gross primary enrollment, and inflation contribute to improvements in the quality of institutions in most estimates, which is consistent with the baseline hypothesis. Trade openness is likely to facilitate the transfer of institutional standards, and an educated population is more aware of its rights and is more capable of asking for better institutions [START_REF] Leite | Does Mother Nature Corrupt? Natural Resources, Corruption and Economic Growth[END_REF][START_REF] Bulte | Resource Intensity, Institutions, and Development[END_REF][START_REF] Kanyama | Quality of Institutions: Does Intelligence Matter?[END_REF]. Higher inflation can be regarded as the failure by a leader, resulting in demands for institutional reforms [START_REF] Aisen | Does Political Instability Lead to Higher Inflation? A Panel Data Analysis[END_REF]. Moreover, the diffusion effect of neighboring countries' institutions is highly significant (except in the GMM). However, this phenomenon deserves a deeper analysis through, for instance, the use of spatial econometric models. Notes: The dependent variable is polity2 score transformed so that it ranges from 0 to 1. F(2, 1,476) 39.11 Prob > F 0.0000

Effects of Different Resource Rents on Institutional Quality

Next we disaggregate total rents into 3 types of rents: forest, mineral and oil. The estimates are reported in Table 4. 7 Regressions 1 to 4 were carried out with the fixed effects estimator. In Column 1, the estimates show that oil rents exert a highly negative effect on the quality of institutions: an increase of 100% in oil rents is associated with a lower polity2 score of 0.15 points. However, forest and mineral rents have no significant impact, as shown in Columns 2 and 3. The same results emerge in Column 4 where all rents are simultaneously introduced into the same regression. This confirms the baseline hypothesis about the specificity of oil rents. Finally, in Columns 5 to 7, we use the GMM-System on 5-year data which allows us to take into account the inertia in institutions. As expected, the size of the oil rents parameter is higher than the one obtained with the fixed effects estimator.

Taken together, the results in Table 4 suggest that different rents have different effects on the quality of institutions. The negative association between total rents and institutional quality seems to be driven by oil rents. These results are consistent with our baseline hypothesis, that in contrast to forest and mineral rents, oil rents are more easily appropriated by elites and are more likely to be a source of social conflicts thereby reducing the quality of institutions.

ROBUSTNESS TESTS

We performed several robustness tests to check whether our results are sensitive to sample selection and model specification. The results are summarized in Table 5. In Columns 1 to 8 we use the variable of Political Rights (PR) from the Freedom House database as an alternative proxy for institutions. This variable is based on an evaluation of three sub-categories of institutional features: (i) the electoral process, (ii) political pluralism and participation, (iii) functioning of the government; its range is between 0 and 7. We transform the PR index so that it ranges from 0 to 1, with higher values indicating better democratic institutions. The results remain broadly unchanged. In particular, total rents and oil rents are negatively associated with the PR index. However, the positive coefficient of forest rents becomes statistically significant with the fixed effects estimator, in Column 2, even though it loses its significance with the System-GMM estimator.

African countries and civil law countries have a lower quality of institutions than the other developing countries. It therefore seems reasonable to test whether the detrimental impact of rents on institutions are specific for these two groups of countries. The African continent is usually characterized by resource richness but poor economic and institutional performances [START_REF] Mehlum | Institutions and the Resource Curse[END_REF]. For this purpose, we use two dummy variables: "Africa" which equals 1 if the country is in Africa, 0 otherwise; and "Legal" which equals 1 if the country has a civil law legal tradition, 0 otherwise (Data from La [START_REF] Porta | The Quality of Government[END_REF], see Table A2). Columns 5 and 6 show that the coefficient of total rents remains negative and statistically significant, while the coefficients of the interaction terms (Rent ´ Africa and Rent ´ Legal) are not statistically significant. These results show that the effect of rents on institutions does not differ from the general case in African or "civil law" countries (see also [START_REF] Mehlum | Institutions and the Resource Curse[END_REF].

In Column 7, we test for the existence of a simple nonlinear relationship between rents and the quality of institutions by including rents squared as an additional regressor. The coefficient of rents remains negative and statistically significant, while the coefficient of rents squared is not statistically different from zero, invalidating this form of nonlinear relationship. Notes: Dependent variable used for regressions of columns (1) to ( 8) is Political rights from Freedom House.

Those of columns ( 9) to ( 11) use the polity2 score as dependent variable. Regressions of columns (1) to ( 4) and ( 9) to ( 11) are carried out with the fixed effects estimator. Those of columns ( 5) to ( 8) are performed with GMM-System estimator. p-values in brackets; * p<0.10, ** p<0.05, *** p<0.01.

CONCLUSION

In this paper, we have investigated the impact of different rents on the quality of institutions using data covering the period 1970-2010 for 90 developing countries. To measure institutional quality we use the polity2 democracy score, and we distinguish total, oil, forest, and mineral rents. We find that total rents are negatively associated with the quality of institutions, after controlling for the other relevant determinants of institutional quality, neighbor effect, persistence of institutions, and endogeneity of rents. However, while oil rents have significant negative effects on institutional quality, forest and mineral rents do not seem to affect institutions. This result holds when we control for time varying common shocks, legal origins, regional dummies, diffusion effects, an alternative measure of institutional quality, and a quadratic specification. Thus, our findings support the hypothesis that oil, unlike mining and forests, generates rents that are easily captured by elites, and leads to the weakening of institutional quality. This could be more thoroughly explored by the use of more detailed data on rents, particularly for mining (Mining has different types of products with different degrees of appropriability), when they are available.

Resource curse literature usually explains that resource richness generates bad performances when institutional quality is poor, suggesting that improving institutions must be -unconditionally -at the top of the agenda. Here we find that the task of improving institutions may face different challenges depending on the kind of resource: typically, oil-rich developing countries seem to be caught in a bad institution trap. For the other resource-rich countries, official development assistance, and a regional strategy (given the neighbor effects), seem to be able to play a key role, whatever the region or the legal origins.

Table A2. Data description and sources

Variable

Description Source polity2 Polity2 scores is defined as the difference between democracy and autocracy scores and it ranges from -10 to 10. It proxies the degree of democracy and measures institutional quality based on (i) the freedom of suffrage, (ii) operational constraints, (iii) balances on executives, and (iv) respect for other basic political rights and civil liberties. We have transformed it so that it ranges from 0 to 1 with higher values indicating better democratic institutions.

legal

A dummy variable that equals 1 if the country has a legal tradition of "Civil law", 0 otherwise. "Civil law" characterizes countries of continental Europe, their former colonies, and Latin America... where the sovereign strongly controls economic life. In contrast the English "Common law" is more based on the Parliament and judges. Legal origins are widely used in the literature on the determinants of institutional quality.

La [START_REF] Porta | The Quality of Government[END_REF]

Figure 1 .

 1 Figure 1. Polity2 and Rents in Percent of GDP

  1800 to 2011), (iii) it is widely used in the literature

::

  Instruments are the ratio of international oil price to population size, and population density.d Instruments are period dummies. The estimate in the last column is performed on five-year average data, reducing the sample size. p-values in brackets; * p<0.10, ** p<0.05, *** p<0.01.

Table 1 .

 1 Variance-covariance Matrix

		polity2 rent	oil	forest Mineral school trade	rem	oda	inflation gdpcapita neighbor
	polity2 1.00
	rent	-0.28* (0.00)	1.00
	oil	-0.38* (0.00)	0.97* (0.00)	1.00
	forest	-0.11* (0.00)	0.16* (0.00)	0.02 (0.49)	1.00
	mineral	0.06* (0.00)	0.22* (0.00)	-0.13* (0.00)	0.09* (0.00)	1.00
	school	0.22* (0.00)	0.13* (0.00)	0.13* (0.00)	-0.33* (0.00)	-0.07* (0.00)	1.00
	trade	-0.01 (0.64)	0.24* (0.00)	0.38* (0.00)	-0.03 (0.17)	0.28* (0.00)	0.16* (0.00)	1.00
	rem	0.03 (0.23)	-0.17* (0.00)	-0.23* (0.00)	-0.01 (0.54)	-0.01* (0.000)	0.07* (0.00)	0.33* (0.00)	1.00
	oda	-0.02 (0.28)	-0.08* (0.00)	0.03 (0.25)	0.55* (0.00)	0.11* (0.00)	-0.33* (0.00)	0.08* (0.00)	0.18* (0.00)	1.00
	inflation	0.01 (0.68)	0.01 (0.64)	-0.02 (0.41)	0.06* (0.01)	0.02 (0.30)	0.01 (0.99)	-0.03 (0.19)	-0.03 (0.12)	-0.01 (0.53)	1.00
	gdpcapita	-0.14* (0.00)	0.41* (0.00)	0.32* (0.00)	-0.36* (0.00)	-0.05* (0.03)	0.18* (0.00)	0.152* (0.00)	-0.113* (0.00)	-0.32* (0.00)	-0.021 (0.31)	1.00
	neighbor	0.53* (0.00)	-0.27* (0.00)	-0.36* (0.00)	-0.15* (0.00)	0.12* (0.00)	0.28* (0.00)	0.08* (0.00)	0.1322* (0.00)	-0.09* (0.00)	0.0329 (0.10)	-0.18* (0.00)	1.00

Table 2 .

 2 Impact of Total Rents on the Quality of Institutions

			Dependent variable: Polity2
		(1)	(2)	(3)	(4)	(5)	(6)	(7)
		OLS	FE	FE	Between	2SLS c	LIML	GMM d
								0.799***
	Polity2 (-1)
								(0.000)
		-0.053***	-0.025***	-0.022***	-0.060***	-0.276***	-0.282***	-0.035**
	log Rent
		(0.000)	(0.001)	(0.006)	(0.009)	(0.000)	(0.000)	(0.037)
		0.003***	0.001	0.002*	0.005*	0.001	0.001	0.004**
	Oda
		(0.005)	(0.364)	(0.080)	(0.090)	(0.656)	(0.661)	(0.024)
		-0.003***	-0.007***	-0.007***	-0.004	-0.006***	-0.006***	-0.006**
	Rem
		(0.000)	(0.000)	(0.000)	(0.278)	(0.000)	(0.000)	(0.046)
		0.073***	0.064**	-0.021	0.128	0.150***	0.152***	0.049
	log school
		(0.001)	(0.019)	(0.462)	(0.186)	(0.000)	(0.000)	(0.513)
		-0.045***	0.145***	0.085***	-0.003	0.231***	0.232***	-0.004
	log trade
		(0.001)	(0.000)	(0.000)	(0.960)	(0.000)	(0.000)	(0.962)
		0.100***	0.029	0.046**	0.183	0.067***	0.068***	-0.069
	log inflation
		(0.000)	(0.182)	(0.021)	(0.376)	(0.009)	(0.009)	(0.595)
		0.021**	0.050***	-0.038**	0.019	-0.028	-0.029	-0.015
	log gdpcapita
		(0.020)	(0.000)	(0.021)	(0.619)	(0.192)	(0.176)	(0.556)
		0.418***	0.407***	0.118***	0.494***	0.353***	0.351***	0.279*
	Neighbor (-1)
		(0.000)	(0.000)	(0.002)	(0.000)	(0.000)	(0.000)	(0.066)
		0.048	-0.852***	0.225	-0.383			-0.062
	constant
		(0.776)	(0.000)	(0.270)	(0.351)			(0.879)
	Adjusted R2	0.280	0.167	0.125	0.246
	Kleibergen-Paap					39.11	39.114
	rk/Wald F statistic a
	Stock-Yogo critical value					19.93	8.68
	AR(1), p-value							0.001
	AR(2), p-value							0.660
	Hansen J stat (p-value)					0.1864	0.1897	0.542
	Overidentification test,
	(p-value) b					0.0000	0.0000
	Year dummies	yes	yes	no	no	no	no	yes
	Observations	1898	1898	1898	1898	1892	1892	445
	Countries	84	84	84	84	84	84	84

Table 3 .

 3 2SLS First-stage Regression

	Dependent variable	log Rent	Number of obs	1,892
			F (9, 1,476)	22.33
			Prob > F	0.000
	Total (centered) SS	594.868	Centered R2	0.114
	Total (uncentered) SS	594.868	Uncentered R2	0.114
	Residual SS	527.008	Root MSE	0.541
		Coef.	t	P>|t|
	log trade	0.360	6.15	0.000
	oda	-0.002	-0.66	0.508
	rem	0.006	2.15	0.032
	log inflation	0.014	2.29	0.029
	log gdpcapita	-0.303	-7.22	0.000
	log school	0.483	5.83	0.000
	lag neighbor	-0.094	0.86	0.392
	Oil price	0.375	6.92	0.000
	density	-0.003	-3.36	0.001
	Included instruments	Log trade, oda, rem, log inflation, log gdpcapita, log school, lag
		neighbor, oil price, density
	F test of instruments

Table 4 .

 4 Impacts of different rents on the quality of institutions

			Dependent variable: polity2
			Fixed Effects			GMM-System
		(1)	(2)	(3)	(4)	(5)	(6)	(7)
	Polity2 (-1)					0.704*** (0.000)	0.701*** (0.000)	0.744*** (0.000)
	log oil rent	-0.012**			-0.012**	-0.033***
		(0.017)			(0.016)	(0.006)
	log forest rent		0.004 (0.772)		0.007 (0.614)		0.061 (0.245)
				-0.001	0.001			0.008
	log mineral rent			(0.661)	(0.681)			(0.547)
	constant	0.122 (0.548)	0.132 (0.540)	0.160 (0.431)	0.093 (0.667)	-0.812 (0.131)	-1.159* (0.053)	-0.931* (0.092)
	Adjusted R2	0.331	0.329	0.329	0.331
	Control	all	all	all	all	all	all	all
	Year dummies	yes	yes	yes	yes	yes	yes	yes
	AR(1), p-value					0.0059	0.0159	0.0134
	AR(2), p-value					0.626	0.761	0.778
	Hansen p-value					0.445	0.178	0.373
	Observations	1,898	1,898	1,898	1,898	445	445	445
	Countries	84	84	84	84	84	84	84
	Notes: The dependent variable is polity2 score transformed such that it ranges from 0 to 1. The GMM-System
	results are two-step estimates with heteroskedasticity-consistent standard errors (Windmeijer, 2005
	correction). AR (1) and AR (2) are tests for null of no first-order and no second-order serial correlation in the
	first-differenced residuals. Hansen J-statistic tests the joint validity of instruments. p-values in parentheses *
	p<0.10, ** p<0.05, *** p<0.01.

Table 5 .

 5 Robustness tests

					Political Rights					Polity2
	Dependent
		(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)
	variable
			Fixed Effects			GMM-System		Fixed Effects
						0.500***	0.517***	0.577***	0.567***
	PolRights (-1)					(0.000)	(0.000)	(0.000)	(0.000)
		-0.0156*				-0.083**				-0.030**	-0.036***	-0.023***
	log Rent	(0.060)				(0.014)				(0.010)	(0.002)	(0.005)
			-0.0082**				-0.045**
	log oil rent		(0.027)				(0.020)
				0.0271*				0.0469
	log forest rent			(0.067)				(0.350)
					0.0012				0.0156
	log mineral rent				(0.625)				(0.218)
										0.0135
	Rent*Africa									(0.353)
											0.0264*
	Rent*Legal										(0.072)
												0.00133
	Rent squared											(0.573)
		0.173	0.526***	-0.0170	0.110	-1.445*	-1.316*	-1.687**	-1.483**	0.238	0.217	0.225
	contant	(0.413)	(0.000)	(0.939)	(0.602)	(0.071)	(0.074)	(0.030)	(0.014)	(0.244)	(0.291)	(0.270)
	Adjusted R2	0.0731 0.100 0.0531 0.0571					0.332	0.334 0.332
	control	all	all	all	all	all	all	all	all	all	all	all
	Year dummies	yes	yes	yes	yes	yes	yes	yes	yes	yes	yes	yes
	AR(1), p-value					0.0070 0.0040 0.0027 0.0010
	AR(2), p-value					0.673 0.467 0.566 0.607
	Hansen,					0.565 0.632 0.555 0.536
	p-value
	Observations	1,902 1,902 1,902	1,902	447	447	447	447	1,898	1,898 ,1898
	Countries	85	85	85	85	85	85	85	85	84	84	84

  PRPolitical Rights index (pr) is a score that aggregates three institutional features: (i) Electoral process, (ii) Pluralism, (iii) Participation and functioning of government. Initially pr index varies between 1 and 7. Transformed such that it ranges from 0 to 1 with higher values indicating better democratic institutions. Rents are defined as the world market price minus the average extraction costs. Unit rents are multiplied by quantities and expressed as a percentage of GDP. "Rent" include all natural resources available in a country. In log. Gross enrolment rate is the share of children of any age that are enrolled in primary school. In log. oda Net official development assistance received as a percentage of GDP. density Density is defined as the number of people per square kilometer.AfricaA dummy variable that equals 1 if the country is in Africa, 0 otherwise rem Remittances received as a percentage of GDP. The world oil prices in current US dollars. Divided by population size to get inter-individual variability.

							PolityIV
							http://www.systemicpeac
							e.org/polity/polity4.htm
	neighbor	Neighbor proxies the institutional quality of neighboring countries
		(measured by polity2 scores).
						n
		it neighbor	=	(1/ ) n	å	inst	jt
					j	1, = ¹ j i
		where	is the number of neighbor countries.
							Freedom House Index
							http://www.freedomhous
							e.org
	rent						Adjusted Net Savings
							Dataset, World Bank
	oil	Log oil rents.
	mineral	Log mineral rents
	forest	Log forest rents
	Trade	The sum of exports and imports as a percentage of GDP. In log.	WDI
	gdpcapita	GDP per capita PPP in current international $. In log.	http://databank.worldban
	inflation	Inflation as measured by the consumer price index annual percentage	k.org/data/views/variable
		change. We have transformed it as follows:	Selection/selectvariables.
			=	(		)	.In log.	aspx?source=world-devel
							opment-indicators
							CNUCED
	oil_price						http://unctadstat.unctad.or
							g/ReportFolders/reportFo
							lders.aspx?sRF_ActivePa
							th=P,5,27&sRF_Expande
							d=,P,5,27

school

We also performed estimates where all control variables are lagged by one period. The results remain unchanged and are available from the authors upon request.

All control variables are maintained but the estimated parameters are not reported in Table4to save space. Complete results are available from the authors upon request.

* We would like to thank the anonymous referee and the editor for their valuable comments. This research was supported by the Agence Nationale de la Recherche of the French government through the program 'Investissements d'avenir' (ANR-10-LABX-14-01), through the IDGM+ initiative led by Ferdi (Fondation pour les Etudes et Recherches sur le Développement International). The usual disclaimers apply. 1

APPENDIX