T Espinasse

F Gamboa

J-M Loubes

Parametric estimation for Gaussian fields indexed by graphs

Keywords: . . . . . . . . . . . 1.2 The . . . . . 1.3 Time series, spectral representation, and MA ∞ . . . .

In this paper, using spectral theory of Hilbertian operators, we study ARM A Gaussian processes indexed by graphs. We extend Whittle maximum likelihood estimation of the parameters for the corresponding spectral density and show their asymptotic optimality.

Contents

1 Definitions and useful properties for spectral analysis and Toeplitz operators 1.1 Graphs, adjacency operator, and spectral representation

Introduction

In the past few years, much interest has been paid to the study of random fields over graphs. It has been driven by the growing needs for both theoretical and practical results for data indexed by graphs. On the one hand, the definition of graphical models by J.N. Darroch, S.L. Lauritzen and T.P. Speed in 1980 [START_REF] Darroch | Markov fields and log-linear interaction models for contingency tables[END_REF] fostered new interest in Markov fields, and many tools have been developed in this direction (see, for instance [START_REF] Verzelen | Tests for Gaussian graphical models[END_REF] and [START_REF] Verzelen | Adaptive estimation of stationary Gaussian fields[END_REF]). On the another hand, the industrial demand linked to graphical problems has risen with the apparition of new technologies. In very particular, the Internet and social networks provide a huge field of applications, but biology, economy, geography or image analysis also benefit from models taking into account a graph structure.

The analysis of road traffic is at the root of this work. Actually, prediction of road traffic deals with the forecast of speed of vehicles which may be seen as a spatial random field over the traffic network. Some work has been done without taking into account the particular graph structure of the speed process (see for example [START_REF] Gamboa | Semi-parametric estimation of shifts[END_REF] and [START_REF] Loubes | Road trafficking description and short term travel time forecasting, with a classification method[END_REF] for related statistical issues). In this paper, we build a new model for Gaussian random fields over graphs and study statistical properties of such stochastic processes.

A random field over a graph is a spatial process indexed by the vertices of a graph, namely (X i ) i∈G , where G is a given graph. Many models already exist in the probabilistic literature, ranging from Markov fields to autoregressive processes, which are based on two general kinds of construction. On the one hand, graphical models are defined as Markov fields (see for instance [START_REF] Guyon | Champs aléatoires sur un réseau[END_REF]), with a particular dependency structure. Actually, they are built by specifying a dependency structure for X i and X j , conditionally to the other variables, as soon as the locations i ∈ G and j ∈ G are connected. For graphical models, we refer for instance to [START_REF] Darroch | Markov fields and log-linear interaction models for contingency tables[END_REF] and references therein. On the other hand, the graph itself, through the adjacency operator, can provide the dependency. This is the case, for example, of autoregressive models on Z d (see [START_REF] Guyon | Champs aléatoires sur un réseau[END_REF]). Here, the local form of the graph is strongly used for statistical inference.

More precisely, the usual purpose of graphical models is to design an underlying graph which reflects the dependency of the data. This method has to be applied when this graph is not easily known (for instance social networks) or when it plays the role of a model which helps understanding the correlations between high complex data (for instance for biological purpose). Our approach differs since, in our case, the graph is known, and we aim at using a model with stationary properties. Indeed, in the case of road traffic, we can consider that the correlations of the process depend mainly on the local structure of the network. This assumption is commonly accepted among professionals of road trafficking speaking of capacity of the road.

In this paper, we extend some classical results from time series to spatial fields over general graphs and provide a new definition for regular ARMA processes on graphs. For this, we will make use of spectral analysis and extend to our framework some classical results of time series. In particular, the notion of spectral density may be extended to graphs. This will enable us to construct a maximum likelihood estimate for parametric models of spectral densities. This also leads to an extension of the Whittle's approximation (see [START_REF] Grenander | Toeplitz forms and their applications[END_REF], [START_REF] Azencott | Series of irregular observations[END_REF]). Actually, many extensions of this approximation have been performed, even in nonstationary cases (see [START_REF] Dahlhaus | Nonparametric quasi-maximum likelihood estimation for Gaussian locally stationary processes[END_REF], [START_REF] Robinson | Multiple local Whittle estimation in stationary systems[END_REF], [START_REF] Giraitis | Whittle estimation of arch models[END_REF]). The extension studied here concerns general ARMA processes over graphs. We point out that we will compare throughout all the paper our new framework with the case G = Z d , d ≥ 1.

Section 1 is devoted to some definitions of graphs and spectral theory for time series. Then we state the definition of general ARMA processes over a graph in Section 2. The convergence of the Whittle maximum likelihood estimate and its asymptotic efficiency are given in Theorems 3.1 and 3.2 in Section 3. Section 4 is devoted to a short discussion on potential applications and perspectives. Some simulations are provided in Section 5. The last section provides all necessary tools to prove the main theorems, in particular Szegö's Lemmas for graphs are given in Section 6.1, while the proofs of the technical Lemmas are postponed in Section 6.3.

1 Definitions and useful properties for spectral analysis and Toeplitz operators 1.1 Graphs, adjacency operator, and spectral representation

In the whole paper, we will consider a Gaussian spatial process (X i ) i∈G indexed by the vertices of an infinite undirected weighted graph. We will call G = (G, W ) this graph, where

• G is the set of vertices. G is said to be infinite as soon as G is infinite (but countable).

• W ∈ [-1, 1] G×G is the symmetric weighted adjacency operator. That is, |W ij | = 0 when i ∈ G and j ∈ G are connected.

We assume that W is symmetric (W ij = W ji , i, j ∈ G) since we deal only with undirected graphs.

For any vertex i ∈ G, a vertex j ∈ G is said to be a neighbor of i if, and only if, W ij = 0. The degree deg(i) of i is the number of neighbors of the vertex i, and the degree of the graph G is defined as the maximum degree of the vertices of the graph G :

deg(G) := max i∈G deg(i).
From now on, we assume that the degree of the graph G is bounded :

deg(G) < +∞. Assume now that W is renormalized : its entries belong to [-1 deg(G) , 1 deg(G)
]. This is not restrictive since re-normalizing the adjacency operator does not change the objects introduced later. In particular, the spectral representation of Hilbertian operator is not sensitive to a renormalization.

Notice that in the classical case G = Z, the renormalized adjacency operator is

W (Z) ij = 1 2 1 1 {|i-j|=1} , (i, j ∈ Z). (1)
Here, deg(Z) = 2. This case will be used in all the paper as an illustration example.

To introduce the spectral decomposition, consider the action of the adjacency operator on l 2 (G) as

∀u ∈ l 2 (G), (W u) i := j∈G W ij u j , (i ∈ G).
We denote by B G the set of all bounded Hilbertian operators on l 2 (G) (the set of square sommable real sequences indexed by G). The operator space B G will be endowed with the classical operator norm

∀A ∈ B G , A 2,op := sup u∈l 2 (G), u 2 ≤1 Au 2 ,
where . 2 stands for the usual norm on l 2 (G). Notice that, as the degree of G and the entries of W are both bounded, W lies in B G , and we have W 2,op ≤ 1.

Recall that for any bounded Hilbertian operator A ∈ B G , the spectrum Sp(A) is defined as the set of all complex numbers λ such that λ Id -A is not invertible (here Id stands for the identity on l 2 (G)). Since W is bounded and symmetric, Sp(W ) is a non-empty compact subset of R [START_REF] Rudin | Functional analysis[END_REF].

We aim now at providing a spectral representation of any bounded normal Hilbertian operator. For this, first recall the definition of a resolution of identity (see for example [START_REF] Rudin | Functional analysis[END_REF]): Definition 1.1. Let M be a σ-algebra over a set Ω. We call identity resolution (on M) a map

E : M → B G such that, 1. E(∅) = 0, E(Ω) = I.
2. For any ω ∈ M, the operator E(ω) is a projection operator.

3. For any ω, ω ′ ∈ M, we have

E(ω ∩ ω ′ ) = E(ω)E(ω ′ ) = E(ω ′ )E(ω). 4. For any ω, ω ′ ∈ M such that ω ∩ ω ′ = ∅, we have E(ω ∪ ω ′ ) = E(ω) + E(ω ′ ).
We can now recall the fundamental decomposition theorem (see for example [START_REF] Rudin | Functional analysis[END_REF])

Theorem 1.1 (Spectral decomposition). If A ∈ B G is symmetric, then there exists a unique identity resolution E over all Borelian subsets of Sp(A), such that

A = Sp(A) λdE(λ).
From the last theorem, we obtain the spectral representation of the adjacency operator W thanks to an identity resolution E over the Borelians of Sp(W )

W = Sp(W ) λdE(λ).
Obviously, we have

W k = Sp(W ) λ k dE(λ), k ∈ N.
Define now, for any i ∈ G, the sequences δ i in l 2 (G) by

δ i := (1 1 k=i ) k∈G .
For any i, j ∈ G, the sequences δ i and δ j define the real measure µ ij by

∀ω ⊂ Sp(W ), µ ij (ω) := E(ω)δ i , δ j l 2 (G) .
Hence, we can write :

∀k ∈ N, ∀i, j ∈ G, W k ij = Sp(W ) λ k dµ ij .
This family of measures µ ij , i, j ∈ G will be used in the whole paper. They convey both spectral information of the adjacency operator, and combinatorial information on the number of path and loops in G. Indeed, the quantity W k ij is the number of path (counted with their weights) going from i to j with length k.

Note also that all diagonals measures µ ii , i ∈ G are probability measures.

The adjacency operator of Z and its spectral decomposition

In the usual case of Z, an explicit expression for µ ij can be given. Denote T k (X) the k th -Chebychev polynomial (k ∈ N). We can provide the spectral decomposition of W (Z) (W (Z) has been defined in Equation 1).

∀i, j ∈ Z, W (Z) k ij = 1 π [-1,1] λ k T |j-i| (λ) √ 1 -λ 2 dλ.
This shows that, in this case, and for any i, j ∈ G, the measure dµ ij is absolutely continuous with respect to the Lebesgue measure, and its density is given by

dµ ij dλ = 1 π T |j-i| (λ) √ 1 -λ 2 .
Notice that we recover the usual spectral decomposition pushing forward µ ij by the function cos : ∀i, j ∈ G, dμ ij (t) := 1 2π cos ((ji)t) dt.

We get

∀i, j ∈ Z, W (Z) k ij = [0,2π] cos(t) k dμ ij (t).
1.3 Time series, spectral representation, and MA ∞

Our aim is to study some kind of stationary processes indexed by the vertices G of the graph G. To begin with, let us recall the usual case of Z. In particular, let us introduce Toeplitz operators associated to stationary time series. Let X = (X i ) i∈Z be a stationary Gaussian process indexed by Z. Since X is Gaussian, stationarity is equivalent to second order stationarity, that is, ∀i, k ∈ Z, Cov(X i , X i+k ) does not depend on i. Thus, we can define

r k := Cov(X i , X i+k ).
Aassume further that (r k ) k∈Z ∈ l 1 (Z). This leads to a particular form of the covariance operator Γ defined on l 2 (Z) by ∀i, j ∈ Z, Γ ij := r i-j .

Recall that B Z denotes here the set of bounded Hilbertian operators on l 2 (Z). Notice that, since (r k ) k∈Z ∈ l 1 (Z), we have Γ ∈ B Z (see for instance [START_REF] Brockwell | Introduction to time series and forecasting[END_REF] for more details). This bounded operator is constant over each diagonals, and is therefore called a Toeplitz operator (see also [START_REF] Böttcher | Analysis of Toeplitz operators[END_REF] for a general introduction to Toeplitz operators). As (r k ) k∈Z ∈ l 1 (Z), we have

∀i, j ∈ Z, T (g) ij := Γ ij = 1 2π [0,2π] g(t) cos ((i -j)t) dt,
where g is the spectral density of the process X, defined by

g(t) := 2 k∈N * r k cos(kt) + r 0 .
This expression can be written, using the Chebychev polynomials (T k ) k∈N ,

g(t) := 2 k∈N * r k T k (cos(t)) + r 0 T 0 (cos(t)) . Let, for λ ∈ [-1, 1], f (λ) := 2 k∈N * r k T k (λ) + r 0 T 0 (λ). (2)
We get, using the family (μ ij ) i,j∈Z defined above,

∀i, j ∈ Z, Γ ij = [0,2π] f (cos(t)) dμ ij (t).
Notice that the last expression may also be written as Γ = f (W (Z) ), and the convergence of the operator valued series defined by Equation 2 is ensured by the boundedness of W (Z) and of the Chebychev polynomials (T k ([-1, 1]) ⊂ [-1, 1], ∀k ∈ Z), together with the summability of the sequence (r k ) k∈Z .

We will extend usual MA processes to any graph, using this previous remark. This will be the purpose of Section 2.

Let us recall some properties about the moving average representation MA ∞ of a process on Z. This representation exists as soon as the log of the spectral density is integrable (see for instance [START_REF] Brockwell | Introduction to time series and forecasting[END_REF]). In this case, there exists a sequence (a k ) k∈N , with a 0 = 1, and a Gaussian white noise ǫ = (ǫ k ) k∈Z ., such that the process X may be written as

∀i ∈ Z, X i = k∈N a k ǫ i-k .
Defining the function h over the unit circle C by

∀x ∈ C, h(x) = k∈N a k x k ,
we recover, with a few computations, the spectral decomposition of the covariance operator Γ of X :

∀i, j ∈ Z, Γ ij = [0,2π] h(e it ) 2 dμ ij (t).
This implies the equality

f (cos(t)) = h(e it ) 2 .
Recall that when h is a polynomial of degree p (with non null first coefficient), the process is said to be MA p . In this case, f is also a polynomial of degree p. Reciprocically, if f is a real polynomial of degree p, and as soon as f (cos(t)) is even, and non-negative for any t ∈ [0, 2π], the Fejér-Riesz theorem provides a factorization of f (cos(t)) such that f (cos(t)) = |h(e it )| 2 (see for instance [START_REF] Kreȋn | The Markov moment problem and extremal problems[END_REF]). This proves that X is MA p if, and only if, its covariance operator may be written f (W (Z) ), where f is a polynomial of degree p.

This remark is fundamental for the construction we provide in the following section (see Definition 2.1).

Whittle maximum likelihood estimation for time series

Here, we recall briefly the Whittle's approximation for time series. Let Θ be a compact interval of R d , d ≥ 1, and (f θ ) θ∈Θ be a parametric family of spectral densities. Let θ 0 ∈ Θ, and assume that (X i ) i∈Z is a Gaussian time series whith spectral density f θ 0 .

If we observe X n := (X i ) i=1,•••n , n > 0, we can define the maximum lokelihood estimate θn of θ 0 as:

θn := arg max L n (θ, X n ), where L n (θ, X n ) := - 1 2 n log(2π) + log det (T n (f θ )) + X T n T n (f θ ) -1 X n .
This estimator is consistent as soon as the spectral densities are regular enough, and under assumptions on the function θ → f θ (see for instance [START_REF] Azencott | Series of irregular observations[END_REF]). However, in practical situations, it is hard to compute. The Whittle's estimate is built by maximizing an approximation of the likelihood instead of the likelihood itself:

θn := arg max Ln (θ, X n ), where Ln (θ, X n ) := - 1 2 n log(2π) + n [0,2π] log (f θ (λ)) dλ + X T n T n ( 1 f θ )X n .
The Whittle estimate is also consistent and asymptotically normal and efficient, as soon as the spectral densities are regular enough.

The consistency of the Whittle estimate relies on the Szegö's Lemma, which provide a bound on the error between 1 n log det (T n (f θ )) and [0,2π] log (f θ (λ)). There exists many versions of this Lemma (see for instance [START_REF] Azencott | Series of irregular observations[END_REF], [START_REF] Grenander | Toeplitz forms and their applications[END_REF]).

In this work, we are interested in a weak version given by Azencott and Dacunha-Castelle in [START_REF] Azencott | Series of irregular observations[END_REF]. The lemma relies on the following fondamental inequality: Let f (x) = k∈N f k x k and g(x) = k∈N g k x k be two analytics function on the complex unitar disk. Then we have

i,j=1,••• ,N T N (f )T N (g) -T N (f g) ij ≤ 1 2 k∈N (k + 1)f k k∈N (k + 1)g k . (3)
In the following, we aim at developing the same kind of tools for processes indexed by a graph.

Spectral definition of ARM A processes

In this section, we will define moving average and autoregressive processes over the graph G.

As explained in the last section, since W is bounded and self-adjoint, Sp(W ) is a nonempty compact subspace of R, and W admits a spectral decomposition thanks to an identity resolution E, given by

W = Sp(W ) λdE(λ).
We define here MA and AR Gaussian processes, with respect to the operator W , by defining the corresponding classes of covariance operators, since the covariance operator fully characterizes any Gaussian process. Definition 2.1. Let (X i ) i∈G be a Gaussian process, indexed by the vertices G of the graph G, and Γ its covariance operator.

If there exists an analytic function f defined on the convex hull of Sp(W ), such that

Γ = Sp(W ) f (λ)dE(λ),
we will say that X is

• MA q if f is a polynomial of degree q.
• AR p if 1 f is a polynomial of degree p which has no root in the convex hull of Sp(W ).

• ARMA p,q if f = P Q with P a polynomial of degree p and Q a polynomial of degree q with no roots in the convex hull of Sp(W ).

Otherwise, we will talk about the MA ∞ representation of the process X. We call f the spectral density of the process X, and denote its corresponding covariance operator by

Γ = K(f ).
Remark Actually, this last construction may also be understood as

Γ = K(f ) = f (W ),
in the sense of normal convergence of the associated power series. However, the spectral representation will be useful in the following. Even if we consider only regular processes in this works, the definition using the spectral representation allows weaker regularity than the definition using the normal convergence of the associated power series.

This kind of modeling is interesting when the interactions are locally propagated (that may be for instance a good modeling for traffic problems.).

The notation K(.) has to be understood by analogy with the notation T (.) used for Toeplitz operators.

Notice that, in the usual case of Z, and for finite order ARMA, we recover the usual definition as shown in Subsection 1.3. So, the last definition may be seen as an extension of isotropic ARMA for any graph G. Besides, note that this extension is given by the equivalence, for any

g ∈ L 2 ([0, 2π]), such that [0,2π] log(g) < +∞, ∀f ∈ L 2 ([-1, 1]), (g = f (cos(t)) ⇔ T (g) = K(f )) .
This means that, in the usual case G = Z, the definition of spectral density in our framework is the usual one, up to an change of variable λ = cos(t) (see Subection 1.3). Now, we get a representation of moving average processes over any graph G. The following section gives the main result of this paper. It deals with the maximum likelihood identification.

Convergence of maximum approximated likelihood estimators

In this section as before, G = (G, W ) is a graph with bounded degree. Let also (X i ) i∈G be a Gaussian spatial process indexed by the vertices of G with spectral density f θ 0 (defined in Section 2) depending on an unknown parameter θ 0 ∈ Θ. We aim at estimating θ 0 . For this, we will generalize classical maximum likelihood estimation of time series. We will also develop a Whittle's approximation for ARMA processes indexed by the vertices of a graph. We follow here the guidelines of the proof given in [START_REF] Azencott | Series of irregular observations[END_REF] for the usual case of time series.

Framework and Assumptions

Let us now specify the framework of our study. Let (G n ) n∈N be a growing sequence of finite nested subgraphs. This means that if

G n = (G n , W n ), we have G n ⊂ G n+1 ⊂ G and that for any i, j ∈ G n , it holds that W n (i, j) = W (i, j).
Let m n = Card(G n ). We set also

δ n = Card {i ∈ G n , ∃j ∈ G\G n , W ij = 0} .
The sequence (m n ) n∈Z may actually be seen as the "volume" of the graph G n , and δ n as the size of the boundary of G n . For the special case

G = Z d and G n = [-n, n] d , we get m n = (2n + 1) d and δ n = 2d(2n + 1) d-1 .
The ratio δn mn is a natural quantity associated to the expansion of the graph that also appears in isoperimetrical [START_REF] Pittet | On the isoperimetry of graphs with many ends[END_REF] and graph expander issues. We will assume here that this ratio goes to 0 when the size of the graph goes to infinity. In short, we set

Assumption 3.1. δ n = o(m n )
This assumption is a non-expansion criterion. The graph has to be amenable, which is satisfied for the last examples G = Z d and G n = [-n, n] d , but not for a homogeneous tree, whatever the choice of the sequence of subgraphs (G n ) n∈N is.

We will now choose a parametric family of covariance operators of MA processes as defined in the last section. First, let Θ be a compact interval of R.

We point out that for sake of simplicity, we choose a one-dimensional parameter space Θ. Nevertheless, all the results could be easily extended to the case Θ ⊂ R k , k ≥ 1.

Define F as the set of positive analytic functions over the convex hull of Sp(W ). Let also (f θ ) θ∈Θ be a parametric family of functions of F . They define a parametric set of covariances on G (see Section 2) by

K(f θ ) = f θ (W ).
As in [START_REF] Azencott | Series of irregular observations[END_REF], we will need a strong regularity for this family of spectral densities. Let us introduce a regularity factor for any analytic function

f ∈ F , f (x) = k f k x k (x ∈ Sp(W )) , by setting α(f ) := k∈N |f k | (k + 1). (4)
Now, let ρ > 0 and define,

F ρ := {f ∈ F , α(log(f )) ≤ ρ} . (5)
Notice that for any f ∈ F ρ , we have α(f ) ≤ e ρ , α( 1 f ) ≤ e ρ . We need the following assumption Assumption 3.2.

• The map θ → f θ is injective.

• For any λ ∈ Sp(W ), the map θ → f θ (λ) is continuous.

• ∀θ ∈ Θ, f θ ∈ F ρ .
From now on, consider θ 0 ∈ Θ. Let X be a centered Gaussian MA ∞ process over G with covariance operator K(f θ 0 ) (see Section 2).

We observe the restriction of this process on the subgraph G n defined before. Our aim is to compute the maximum likelihood estimator of θ 0 . Let X n = (X i ) i∈Gn be the observed process and K n (f θ ) be its covariance :

X n ∼ N (0, K n (f θ 0 )) .
The corresponding log-likelihood at θ is

L n (θ) := - 1 2 m n log(2π) + log det (K n (f θ )) + X T n K n (f θ ) -1 X n .
As discussed before, in the case G = Z, it is usual to maximize an approximation of the likelihood. The classical approximation is the Whittle's one ( [START_REF] Grenander | Toeplitz forms and their applications[END_REF]), where

1 n log det (T n (g)) is replaced by 1 2π [0,2π] log (g (t)) dt.
Back to the general case, we aim at performing the same kind of approximation. For this, we will need the following assumption to ensure the convergence of log det (K n (f θ )) (see Section 1 for the definition of µ ii ) : Assumption 3.3. There exists a positive measure µ, such that

1 m n i∈Gn µ ii D -→ n→∞ µ.
Here, D stands for the convergence in distribution

The limit measure µ is classically called the spectral measure of G with respect to the sequence of subgraphs (G n ) n∈Z (see [START_REF] Mohar | A survey on spectra of infinite graphs[END_REF] for example).

Actually, under Assumption 3.1, Assumption 3.3 is equivalent to the convergence of the empirical distribution of eigenvalues of W Gn (here, W Gn denotes the restriction of

W over the subgraph G n ) That is, if λ (n) 1 , • • • , λ (n)
mn denote the eigenvalues (written with their multiplicity orders) of W gn , Define µ [1] n :=

1 m n mn i=1 δ λ (n) i , and
µ [2] n = 1 m n i∈Gn µ ii ,
Then, under Assumption 3.1, the convergence of µ [START_REF] Aldous | Processes on unimodular random networks[END_REF] n to µ (i.e. Assumption 3.3) is equivalent to the convergence of µ [START_REF] Azencott | Series of irregular observations[END_REF] n to µ.

To prove this equivalence we just have to notice that :

Sp(W )
λ k dµ (1) n (λ) -

Sp(W ) λ k dµ (2) n (λ) = 1 m n mn i=1 λ (n) k i - 1 m n i∈Gn (W k ) ii = 1 m n Tr (W Gn ) k - 1 m n Tr (W k ) Gn .
So that, we get the result by Lemma 6.1 (see Section 6.1). As in the case of time series (for G = Z), we can approximate the log-likelihood. It avoids an inversion of a matrix and a computation of a determinant. Indeed, we will consider the two following approximations.

Ln (θ) := - 1 2 m n log(2π) + m n log(f θ (x))dµ(x) + X T n (K n (f θ )) -1 X n . Ln (θ) := - 1 2 m n log(2π) + m n log(f θ (x))dµ(x) + X T n K n 1 f θ X n .
Notice that approximated maximum likelihood estimators are not asymptotically normal in general (see for instance [START_REF] Guyon | Parameter estimation for a stationary process on a d-dimensional lattice[END_REF] for Z d ). Indeed, the score associated to the approximated log-likelihood has to be asymptotically unbiased [START_REF] Azencott | Series of irregular observations[END_REF].

To overcome this problem in Z d , the tapered periodogram can be used (see [START_REF] Guyon | Champs aléatoires sur un réseau[END_REF], [START_REF] Guyon | Parameter estimation for a stationary process on a d-dimensional lattice[END_REF], [START_REF] Dahlhaus | Edge effects and efficient parameter estimation for stationary random fields[END_REF]).

Let us consider graph extensions of standard time series models :

• The MA P case : There exists P > 0 such that the true spectral density f θ 0 is a polynomial of degree bounded by P .

• The AR P case : There exists P > 0 such that all the spectral densities (for any θ ∈ Θ) of the parametric set are such that 1 f θ is a polynomial of degree bounded by P . So, to define the good approximated log-likelihood, we first introduce the unbiased periodogram in each of the last cases. Now, let P > 0.

Define a subset V P of signed measures on R as

V P := {µ ij , i, j ∈ G, d G (i, j) ≤ P } ,
where d G (i, j), i, j ∈ G stands for the usual distance on the graph G, i.e. the length of the shortest path going from i to j. We will need the following assumption Assumption 3.4. The set V P of possible local measures over G is finite, and n is large enough to ensure that ∀v ∈

V P , ∃(i, j) ∈ G 2 n , µ ij = v.
Remark This assumption is quite strong, and holds for instance for quasi-transitive graphs (i.e. such that the quotient of the graph with its automorphism group is finite). This assumption may be relaxed, but it is a hard and technical work that will be the issue of a forthcoming paper.

Define now the matrix B (n) (the dependency on P is omitted, for clarity) by

B (n) ij := Card {(k, l) ∈ G n × G, µ kl = µ ij } Card {(k, l) ∈ G n × G n , µ kl = µ ij } , if , d G (k, l) ≤ P := 1 if d G (k, l) > P.
The matrix B (n) gives a boundary correction, comparing, for any v ∈ V P the frequency of the interior couples of vertices with local measure v with the boundary couples of vertices with local measure v. Actually, this way to deal with the edge effect is very similar to the one used for G = Z d (see [START_REF] Dahlhaus | Edge effects and efficient parameter estimation for stationary random fields[END_REF], [START_REF] Guyon | Parameter estimation for a stationary process on a d-dimensional lattice[END_REF]).

As example, let us now describe the case G = Z 2 , for P = 2. In this case

W (Z 2 ) is ∀i, j, k, l ∈ Z, W (Z 2 ) ((i, j), (k, l)) := 1 4 1 1 |i-j|+|k-l|=1 .
In this example, we set G n = [1, n] 2 , and we can compute the matrix B (n) . Indeed, it only is needed to notice that

µ (i 1 ,j 1 ),(i 1 +k,j 1 +l) = µ (i 2 ,j 2 ),(i 2 +ǫ 1 k,j 2 +ǫ 2 l) , i 1 , i 2 , j 1 , j 2 , k, l ∈ Z, ǫ 1 , ǫ 2 ∈ {-1, 1} .
This means that the local measure of a couple of vertices depends only of their relative positions (stationarity and isotropy of this set of measure). So, we need to count the configurations given by Figure 1 since we consider only couples of vertices u, v ∈ Z 2 such that

d Z 2 (u, v) ≤ 2.
We get, for any i, j ∈ Z,

• B (n) (i,j),(i,j) = n 2 n 2 = 1. • B (n) (i,j),(i,j±1) = B (n) (i,j),(i±1,j) = 4n(n-1) 4n 2 . • B (n) (i,j),(i±1,j±1) = 4(n-1) 2 n 2 . • B (n) (i,j),(i,j±2) = B (n) (i,j),(i±2,j) = 4n(n-2) 4n 2
One can notice that sup

ij B (n) ij -1 → n→∞ 0.
Assumption 3.5 ensure that this property holds for the graph we consider. Back to the general case, let f ∈ F ρ . We define the unbiased periodogram as

X T n Q n ( 1 f )X n .
where

Q n (f ) := B (n) ⊙ K n (f ).
Here, the operation ⊙ denotes the Hadamard product for matrices, that is

∀i, j ∈ G n , B (n) ⊙ K n (f ) ij = B (n) ij K n (f ) ij .
Notice that this is actually a way to extend the so called tapered periodogram (see for instance [START_REF] Guyon | Parameter estimation for a stationary process on a d-dimensional lattice[END_REF]).

We now define the unbiased empirical log-likelihood, for any θ ∈ Θ

L (u) n (θ) := - 1 2 m n log(2π) + m n log(f θ (x))dµ(x) + X T n Q n ( 1 f θ ) X n .
We denote by θn , θn , θn , θ (u) the maximum likelihood estimators associated to L n , Ln , Ln , L (u) n , respectively. We will need the following assumption, Assumption 3.5. There exists a positive sequence (u n ) n∈N such that,

u n → n→∞ 0,
and

sup ij B (n) ij -1 ≤ u n .
Notice that the last assumption holds for example in the case G = Z d , d > 1.

To prove asymptotic normality and efficiency of the estimator θ

(u)
n , we will also need the following assumption.

Assumption 3.6. Assume that

• There exists a positive sequence

(v n ) n∈N such that v n = o( 1 √ mn ) and ∀f ∈ F ρ , 1 m n Tr(K Gn (f )) -f dµ ≤ α(f )v n .
• For any θ ∈ Θ, f θ is twice differentiable on Θ and

d dθ (f θ ) ∈ F ρ , d 2 dθ 2 (f θ ) ∈ F ρ .
The first assumption means that the convergence of the empirical distribution of eigenvalues of K(f ) to the spectral measure µ is faster than 1 √ mn . It holds for instance for quasi-transitives graphs, with a suitable sequence of subgraphs. The second assumption is more classical. For example it is required in the case G = Z (see [START_REF] Azencott | Series of irregular observations[END_REF]).

Convergence and asymptotic optimality

Let ρ > 0. We can now state one of our main result: Theorem 3.1. Under Assumptions 3.1, 3.2 and 3.3, the sequences ( θn ) n∈N , ( θn ) n∈N , ( θn ) n∈N converge, as n goes to infinity, P f θ 0 -a.s. to the true value θ 0 . If moreover Assumption 3.5 holds, this is also true for (θ

(u) n ) n∈N .
Proof. The proof follows the guidelines of [START_REF] Azencott | Series of irregular observations[END_REF]. We highlight the main changes performed here. First, we define the Kullback information on G n of f θ 0 with respect to f ∈ F ρ , by

IK n (f θ 0 , f ) := E P f θ 0 -log( dP f dP f θ 0 ) .
and the asymptotic Kullback information (on G) by

IK(f θ 0 , f ) = lim n 1 m n IK n (f θ 0 , f )
whenever it is finite. The convergence of the estimators of the maximum approximated likelihood is a direct consequence of the following lemmas : Lemma 3.1. For any f ∈ F ρ , and under Assumptions 3.1, 3.2 and 3.3, the asymptotic Kullback information exists and may be written as

IK(f θ 0 , f ) = 1 2 -log( f θ 0 f ) -1 + f θ 0 f dµ.
Furthermore, if we set l n (θ, X n ) = 1 mn L n (θ, X n ), we have that P f θ 0 -a.s.,

l n (θ 0 , X n ) -l n (θ, X n ) → n→∞ IK(f θ 0 , f θ )
uniformly in θ ∈ Θ. This property also holds for ln := 1 mn Ln and ln := 1 mn Ln Furthermore, for P > 0, and for both the AR P or the MA P case (see above), this also holds for l (u)

n := 1 mn L (u)
n . Lemma 3.2. Let f θ 0 be the true spectral density, and (ℓ n ) n∈N be a deterministic sequence of continuous functions such that

∀θ ∈ Θ, ℓ n (θ 0 ) -ℓ n (θ) → n→∞ IK(f θ 0 , f θ )
uniformly as n tends to infinity. Then, if θ n = arg max θ ℓ n (θ), we have

θ n → n→∞ θ 0 .
The proofs of these lemmas are postponed in Appendix (Subsection 6.2). Theorem 3.2. In both the AR P or MA P cases, and and under all previous assumptions 3.1, 3.2, 3.3, 3.4, 3.5, 3.6, the estimator θ (u) n of θ 0 is asymptotically normal:

√ m n (θ (u) n -θ 0 ) D -→ n→∞ N 0, 1 2
f ′ θ 0 f θ 0 2 dµ -1
.

Furthermore, the Fisher information of the model is

J(θ 0 ) := 1 2 f ′ θ 0 f θ 0 2 dµ.
Hence, the previous estimator is asymptoticly efficient.

Proof. Here again, we mimic the usual proof by extending the result of [START_REF] Azencott | Series of irregular observations[END_REF] to the graph case. Using a Taylor expansion, we get

(l (u) n ) ′ (θ 0 ) = (l (u) n ) ′ (θ (u) n ) + (θ 0 -θ (u) n )(l (u) n ) ′′ ( θn ),
where θn ∈ θ

(u) n , θ 0 . As θ (u) n = arg max l (u) n , we have (l (u) n ) ′ (θ (u) n ) = 0. So that, √ m n (θ 0 -θ (u) n ) = (l (u) n ) ′′ ( θn ) -1 √ m n (l (u) n ) ′ (θ 0
). The end of the proof relies on three lemmas : Lemma 3.3 provides the asymptotic normality for √ m n (l (u) n ) ′ (θ 0 ). Combined with Lemma 3.4, we get the asymptotic normality for √ m n (θ 0 -θ (u) n ). Finally, Lemma 3.5 gives the Fisher information.

Lemma 3.3. √ m n (l (u) n ) ′ (θ 0 ) D -→ n→∞ N 0, 1 2
f ′ θ 0 f θ 0 2 dµ . Lemma 3.4. (l (u) n ) ′′ ( θn ) -1 → n→∞ 2 f ′ θ 0 f θ 0 2 dµ -1
, P f θ 0a.s.

Lemma 3.5. The asymptotic Fisher information is :

J(θ 0 ) = 1 2 f ′ θ 0 f θ 0 2 dµ.
The proofs of these lemmas are postponed in Appendix (Subsection 6.3)

Discussion

Note first that Theorem 3.1 provides consistency of the estimators under weak conditions on the graph. Indeed, amenability ensures Assumption 3.1, for a suitable sequence of subgraphs. Assumption 3.3 holds as soon as there is a kind of homogeneity in the graph. The simplest application is quasi-transitives graph. Note that if G is "close" to be quasi-transitive, Assumption 3.3 is still true. We also could adapt notions of unimodularity [START_REF] Aldous | Processes on unimodular random networks[END_REF] or stationarity [START_REF] Benjamini | Ergodic theory on stationary random graphs[END_REF] to our framework and prove the existence of a spectral measure. Furthermore, Assumption 3.3 holds for the real traffic network (this will be explained in a forthcomming paper).

To build the estimator θ

(u)
n , stronger assumptions on the graph G are needed. Let us discuss two very special cases. First, Theorem 3.2 may be applied in the Z d case with holes, that is in the presence of missing data, up to the condition that they remain few enough. Actually, Assumption 3.1 is required, so the boundary of the subgraphs (counting the holes) has to be small in front of the volume of this subgraphs.

We need furthermore a kind of homogeneity for these holes. For instance, we can assume that the data are missing completely at random. This particular case is interesting for prediction issues.

Another strong potential application is quasi-transitive graphs, as mentioned above. Indeed, take for instance a finite graph (the pattern) and reproduce it at each vertex of an infinite (amenable) vertex-transitive graph. The final graph is then quasi-transitive, and all the previous assumptions hold.

This seems to be a natural extension of what happens for Z d . Furthermore, in this situation as in Z d , our work may also be applied to a process with missing values.

Note also that conditions of both amenability of the graphs and regularity of spectral densities seem natural, looking at the Szegö's Lemmas (see Section 6.1). Indeed, the difference computed in Lemma 6.1 is only due to edge effects. Thus, there are two ways for relaxing this conditions. On the one hand, it could be interesting to deal with lower regularity (for instance to study long memory processes) for the spectral densities. On the other hand, it could be also interesting to relax conditions on the graph, for instance for more regular densities. In particular, we could investigate the case of random graphs, and try to pick up homogeneity conditions into the random structure. As mentioned above, another natural extension of this work could be done to graphs "close" to be quasi-transitive.

These two limits of our present work are actually two of our main perspectives in this framework.

Simulations

In this section, we give some simulations over a very simple case, where the graph G is built taking some rhombus connected by a simple edge both on the left and right (see Figure 2). The sequence of nested subgraphs chosen here is the growing neighborhood sequence (we chose a point x and we take G n = {y ∈ G, d G (x, y) ≤ n}). We study an AR 2 model, where,

Θ = ]-1, 1[ , f θ (x) = 1 1 -θx 2 (θ ∈ Θ).
Here, we take for W the adjacency operator of G normalized in order to get sup i,j∈G W ij ≤ 1 deg(G) . We choose θ 0 = 1 2 , m n = 724. We approximate the spectral measure of G by the spectral measure of a very large graph (around 10000 vertices) built in the same way. Figure 3 shows the empirical spectrum of the graph G with respect to the sequence of subgraphs (G n ) n∈N .

To compute (K n (f θ )) -1 , we use the power series representation of f θ , and truncate this expression after the 15 first coefficient. This choice ensures that the simulation errors are neglectible with respect to the theoretical ones.

Figure 4 gives the empirical distribution of

√ m n Sp(A) f ′ θ f θ 2 θn -θ 0 .

Szegö's Lemmas

Szegö's Lemmas [START_REF] Grenander | Toeplitz forms and their applications[END_REF] are useful in time series analysis. Indeed, they provide good approximations for the likelihood. As explained in Section 3, these approximations of the likelihood are easier to compute.

In this section, we generalize a weak version of the Szegö Lemmas, for a general graph, under Assumption 3.1 (non expansion criterion for G n ), and Assumption 3.3 (existence of the spectral measure µ).

For any matrix (B ij ) i,j∈Gn , we define the block norm

b N (B) = 1 δ N i,j∈G N |B ij | .
We can state the equivalent version of the first Szegö lemma for time-series Lemma 6.1. Asymptotic homomorphism Let k, n be positive integers, and let

g 1 , • • • , g k be analytic functions over [-1, 1] having finite regularity factors (i.e. α(g i ) < +∞, i = 1, • • • , k). Then, b n (K n (g 1 ) • • • K n (g k ) -K n (g 1 • • • g k )) ≤ k -1 2 α(g 1 ) • • • α(g k ).
Corollary 6.1. For any g ∈ F ρ (see the first page of Subsection 3.1 for the definition), and under Assumptions 3.1 and 3.3,

1 m n log det(K n (g)) → n→∞ log(g)dµ.
Proof. of Lemma 6.1 This proof follows again the one of [START_REF] Azencott | Series of irregular observations[END_REF]. We will prove the result by induction on k.

First we deal with the case k = 2. Let f and g analytic functions over [-1, 1] such that α(f ) < +∞ and α(g) < +∞. We write

b n (K n (f )K n (g) -K n (f g)) = 1 δ n i,j∈Gn k∈Gn (K n (f )) ik (K n (g)) kj - k∈G (K n (f )) ik (K n (g)) kj = 1 δ n i,j∈Gn k∈G\Gn |K(f ) ik | |K(g) kj | .
Using K(g) = ∞ h=0 g h W h , Fubini's theorem gives, since all the previous sequences are

in l 1 (G), b n (K n (f )K n (g) -K n (f g)) ≤ 1 δ n i,j∈Gn k∈G\Gn (K n (f )) ik (K n (g)) kj ≤ sup k∈G\Gn i∈Gn |K(f ) ik | × 1 δ n k∈G\Gn j∈Gn ∞ h=0 |g h | (W h ) kj ≤ sup k∈G i∈G |K(f ) ik | × ∞ h=0 |g h | 1 δ n k∈G\Gn j∈Gn (W h ) kj . Introducing ∆ h = sup N ∈N 1 δ N k∈G\G N j∈G N W h kj , we get b n (K n (f )K n (g) -K n (f g)) ≤ sup k∈G i∈G |K(f ) ik | ∞ h=0 |g h | ∆ h .
The coefficient ∆ h is a porosity factor. It measures the weight of the paths of length h going from the interior of G n to outside.

Note that ∆ h ≤ h + 1, so we get

∞ h=0 |g h | ∆ h ≤ α(g).
Now, we define another norm on B G :

B ∞,in := sup k∈G i∈G |B ik | , (B ∈ B G ) .
We thus obtain

K(f ) ∞,in = sup k∈G i∈G |K(f ) ik | ≤ ∞ h=0 |f h | W h ∞,in ≤ ∞ h=0 |f h | W h ∞,in ≤ ∞ h=0 |f h | := f 1,pol .
Finally, we get

b n (K Gn (f )K Gn (g) -K Gn (f g)) ≤ f 1,pol α(g).
To conclude the proof of the lemma, by symmetrization of the last inequality, and since 1 ≤ (h + 1), we have,

b n (K n (f )K n (g) -K n (f g)) ≤ 1 2 α(f )α(g). (6)
To perform the inductive step, we need the following inequalities [START_REF] Seber | A matrix handbook for statisticians[END_REF]:

α(f g) ≤ α(f )α(g), b n (BC) ≤ B ∞,in b n (C), b n (B + C) ≤ b n (B) + b n (C), K n (f ) ∞,in = f 1,pol ≤ α(f ).
Let k > 1, and assume that for all j ≤ k -1, Lemma 6.1 holds. Under the previous assumptions, and the inductive hypothesis for k -1 we get,

b n (K n (g 1 )× • • • ×K n (g k ) -K n (g 1 • • • g k )) ≤ K n (g 1 ) ∞,in b n (K n (g 2 ) • • • K n (g k ) -K n (g 2 • • • g k )) +b n (K n (g 1 )K n (g 2 • • • g k ) -K n (g 1 • • • g k )) ≤ α(g 1 ) k -2 2 α(g 2 ) • • • α(g k ) + 1 2 α(g 1 )α(g 2 • • • g k ) ≤ k -1 2 α(g 1 ) • • • α(g k ),
which completes the induction step and proves the result.

Proof. of Corollary 6.1 Let g ∈ F ρ , and k be a positive integer. Using Lemma 6.1, we have

Tr K n (g) k -K n (g k ) ≤ δ n m n b n K n (g) k -K n (g k ) . (7)
Thus, we have, thanks to Assumption 3.1

1 m n Tr K n (g) k -K n (g k ) → n→+∞ 0. Denote µ [1]
g the real measure whose k th -moment is given by

x k dµ [1] g = lim n 1 m n Tr K n (g) k ,
and µ

[2]

g the real measure whose k th -moment is given by

x k dµ [2] g = lim n 1 m n Tr K n (g k ) .
Notice that both of these measures have support between inf g ≥ e -ρ > 0 and sup g ≤ e ρ < +∞, since α(log(g)) < ρ (see Section 3). Therefore, the equality of the moments given by Equation 7gives the equality of the measures µ [START_REF] Aldous | Processes on unimodular random networks[END_REF] g and µ [START_REF] Azencott | Series of irregular observations[END_REF] g .

So that, we get

1 m n log (det (K n (g))) - 1 m n Tr (K n (log(g))) → n→+∞ 0. (8)
Assumption 3.3 completes the proof of the Corollary since it implies that

1 m n Tr (K n (log(g))) → n→+∞ log(g)dµ.
The following lemma enables to replace K n (g) by the unbiased version Q n (g) (see Section 3 for the definition). Lemma 6.2. Under Assumptions 3.1,3.3, 3.4 and 3.5, and if f or g is a polynomial having degree less than or equal to P , we have

1 m n Tr ((K n (f )K n (g)) p -(K n (f )Q n (g)) p ) ≤ 2 p u n α(f ) p α(g) p .
Proof. We define, for any f ,

f abs (x) = k |f k | x k .
Actually, the proof is based of the following idea: as soon as f or g is a polynomial having degree less than or equal to P , we have to control only the number of paths of length less than or equal to P (counted with their weights).

Let p be a positive number. Recall that Q n ( 1 g ) = B (n) ⊙ K n ( 1 g ) (see Section 3), we have,

1 mn Tr K n (f )K n ( 1 g ) p -K n (f )Q n ( 1 g ) p ≤ 1 m n i∈Gn i 0 =i,i 1 ,••• ,i 2p =i l=0•••p B (n) i 2l i 2l+1 K n ( 1 g ) i 2l i 2l+1 K n (f ) i 2l+1 i 2l+2 - 1 m n i∈Gn i 0 =i,i 1 ,••• ,i 2p =i l=0•••p K n ( 1 g ) i 2l i 2l+1 K n (f ) i 2l+1 i 2l+2 ≤ 1 m n sup i 1 ,i 2 ,••• ,i 2p+1 l=0•••p-1 B (n) i 2l+1 i 2l+2 -1 × i∈Gn i 0 =i,i 1 ,••• ,i 2p =i l=0•••p K n ( 1 g ) i 2l i 2l+1 K n (f ) i 2l+1 i 2l+2 ≤ 1 m n sup i 1 ,i 2 ,••• ,i 2p+1 l=0•••p-1 B (n) i 2l+1 i 2l+2 -1 × i∈Gn i 0 =i,i 1 ,••• ,i 2p =i l=0•••p K n (( 1 g ) abs ) i 2l i 2l+1 K n (f abs ) i 2l+1 i 2l+2 ≤ sup i 1 ,i 2 ,••• ,i 2p+1 l=0•••p-1 B (n) i 2l+1 i 2l+2 -1 K Gn (f abs )K Gn (( 1 g ) abs ) p 2,in ≤ sup i 1 ,i 2 ,••• ,i 2p+1 l=0•••p-1 B (n) i 2l+1 i 2l+2 -1 α(f ) p α( 1 g ) p .
Using Assumption 3.5, we get,

1 mn Tr K n (f )K n ( 1 g ) p -K n (f )Q n ( 1 g ) p ≤ |(1 + u n ) p -1| α(f ) p α( 1 g ) p ≤ (1 + u n -1) (1 + u n ) p-1 + (1 + u n ) p-2 + • • • + 1 α(f ) p α( 1 g ) p ≤ |u n (2 p -1)| α(f ) p α( 1 g ) p ≤ u n 2 p α(f ) p α( 1 g ) p .
This ends the proof of the Lemma.

Finally, the following lemma explains the choice of B (n) . The unbiased quadratic form Q n is no more than a correction of the error between K n (f )K n (g) and K n (f g). Lemma 6.3 (Exact correction). Let f, g ∈ F ρ , and assume that either f or g is a polynomial of degree less than or equal to P (see Section 3). Then, the unbiased quadratic form

Q n (f θ ) verify Tr (K n (f )Q n (g)) = Tr (K n (f g)) .
Proof. of Lemma 6.3 First, notice that

Tr (K n (f )Q n (g)) = i,j∈Gn K n (f ) ij K n (g) ij B (n) ij .
Since this expression is symmetric on f, g, we can now consider the case where f is a polynomial of degree less than or equal to P . Actually, since f is a polynomial, K n (f ) ij = 0 as soon as d(i, j) > P (i, j ∈ G). Then, if i, j, k, l ∈ G are such that µ ij = µ kl , we have

K n (f ) ij K n (g) ij = K n (f ) kl K n (g) kl .
So that, we may here denote, for convenience, K(f ) µ ij .

Using Assumption 3.4, this leads to

Tr (K n (f )Q n (g)) = i,j∈Gn K n (f ) ij K n (g) ij B (n) ij = v∈V P i,j∈Gn µ ij =v,d G (i,j)≤P K n (f ) v K n (g) v B (n) v = v∈V P K n (f ) v K n (g) v Card {(i, j) ∈ G n × G n , µ ij = v} × Card {(i, j) ∈ G n × G, µ ij = v} Card {(i, j) ∈ G n × G n , µ ij = v} , = v∈V P (i,j)∈Gn×G, µ ij =v,d G (i,j)≤P K n (f ) v K n (g) v B (n) v = (i,j)∈Gn×G K n (f ) ij K n (g) ij B (n) ij = Tr (K n (f g)) .
That ends the proof of Lemma 6.3.

Proofs of the lemmas of Theorem 3.1

Recall that the theorem relies on two lemmas. Lemma 3.2 states a condition on deterministic sequences to provide the convergence of the maximizer of these sequences.

Proof. of Lemma 3.2 Recall that f θ 0 denotes the true spectral density. Let (ℓ n ) n∈N be a deterministic sequence of continuous functions such that

∀θ ∈ Θ, ℓ n (θ 0 ) -ℓ n (θ) → n→∞ 1 2 -log( f θ 0 f θ ) -1 + f θ 0 f θ dµ. ( 9
)
uniformly as n tends to infinity. Denotes moreover θ n = arg max θ ℓ n (θ). We aim at proving that θ n → n→∞ θ 0 .

Using the compactness of Θ, let θ ∞ be an accumulation point of the sequence (θ n ) n∈N , and (θ n k ) k∈N be a subsequence converging to θ ∞ . As the function

θ → 1 2 -log( f θ 0 f θ ) -1 + f θ 0 f θ dµ
is continuous on Θ, and the convergence of (ℓ n (θ 0 )ℓ n (θ)) n∈N is uniform in θ, we have

ℓ n k (θ 0 ) -ℓ n k (θ n k ) k→∞ ---→ 1 2 -log( f θ 0 f θ∞ ) -1 + f θ 0 f θ∞ dµ. ( 10
)
But we can notice that, thanks to the definition of θ n , ℓ n k (θ 0 )ℓ n k (θ n k ) ≤ 0 So, since the function x →log(x) + x -1 is non negative and vanishes if, and only if, x = 1, we get that f θ 0 = f θ∞ . By injectivity of the function θ → f θ , we get θ ∞ = θ 0 , for any accumulation point θ ∞ of the sequence (θ n ) n∈N , which ends the proof of this first lemma.

Lemma 3.1 provides the uniform convergence of the contrasts of maximum likelihood and approximated maximum likelihood to the Kullback information. The proof may be cut into several lemmas.

Proof. of Lemma 3.1

First, notice that by construction, we have, for any θ ∈ Θ,

IK(f θ 0 , f θ ) = lim n E 1 m n (L n (f θ 0 , X n ) -L n (f θ , X n )) , (11)
when it exists. Then, we can compute

l n (f θ 0 , X n ) -l n (f θ , X n ) = - 1 2m n (log det(K n (f θ 0 )) -log det(K n (f θ ))) - 1 2m n X T n K n (f θ 0 ) -1 X n -X T n K n (f θ ) -1 X n
Corollary 6.1 of Lemma 6.1 provides the following convergence

1 m n (log det(K n (f θ 0 )) -log det(K n (f θ ))) → n→∞ log f θ 0 f θ dµ. ( 12
)
To prove the existence of IK(f θ 0 , f θ ), it only remains to prove the P f θ 0 -a.s. convergence

of 1 mn X T n K n (f θ ) -1 X n to f θ 0
f θ dµ as n goes to infinity. This is ensured by the following Lemma. Lemma 6.4 (Convergence lemma). For respectively

Λ = K n ( 1 f θ ), Λ = (K n (f θ )) -1 or Λ = Q n ( 1 f θ ), we have, 1 m n X T n ΛX n → n→∞ f θ 0 f θ dµ, P f θ 0 -a.s..
Lemma 6.4 combined with Corollary 6.1 ensures the P f θ 0a.s. convergence of ln (f θ 0 )ln (f θ ), ln (f θ 0 ) -ln (f θ ) to IK(f θ 0 , f θ ). It provides also the P f θ 0a.s. convergence of l

(u) n (f θ 0 ) - l (u)
n (f θ ) to IK(f θ 0 , f θ ) in the AR P or MA P cases (see Section 3). To complete the assertion of Lemma 3.1, it only remains to show the uniform convergences on Θ of the last quantities. This will be done using an equicontinuity argument given by the following Lemma. Lemma 6.5 (Equicontinuity lemma). For all n ≥ 0, the sequences of functions

(l n (f θ 0 , X n ) -l n (f θ , X n )) n∈N
is an P f θ 0 -a.s. equicontinuous sequence on ({f θ , θ ∈ Θ} , . ∞ ).This property also holds for ln , ln . Furthermore, the sequence l

(u) n (f θ 0 , X n -l (u) n (f θ , X n ) n∈N is also P f θ 0 -a.s. equicontin- uous, on {f θ , θ ∈ Θ} , . 1,pol .
We can now end the proof of Lemma 3.1: First, notice that the space {f θ , θ ∈ Θ} is compact for the topology of the uniform convergence. This also holds for {f θ , θ ∈ Θ} , . 1,pol . So, there exists a dense sequence (f θp ) p∈N .

Then, using Lemma 6.1 and Corollary 6.1, the sequence l

n (f θ 0 , X n ) -l n (f θp , X n ) n∈N con- verges P f θ 0 -a.s. to IK(f θ 0 , f θp ).
If a sequence of functions is equicontinuous and converges pointwise on a dense subset of its domain, and if its co-domain is a complete space, then the sequence converges pointwise on all the domain [START_REF] Rudin | Functional analysis[END_REF].

Using this well known property, we obtain, P f θ 0 -a.s., the pointwise convergence of

(l n (f θ 0 , X n ) -l n (f θ , X n )) n∈N
to IK(f θ 0 , f θ ), for any θ ∈ Θ. Furthermore, if a sequence of functions is equicontinuous and converges pointwise on its domain, then this convergence is uniform on any compact subspace of the domain [START_REF] Rudin | Functional analysis[END_REF].

Thus, we get, P f θ 0 -a.s., the uniform convergence on Θ of the sequence

(l n (f θ 0 , X n ) -l n (f θ , X n )) n∈N to IK(f θ 0 , f θ ).
Using the same kind of arguments, this uniform convergence also holds for ln , ln and l

n . This concludes the proof of Lemma 3.1.

Proof of the technical lemmas

Proof. of Lemma 6.4

Let θ ∈ Θ. First, consider the case Λ n = K n 1 f θ . We aim at proving that

1 m n X T n Λ n X n → n→∞ f θ 0 f θ dµ, P f θ 0 -a.s..
To do that, we make use of classical tools of large deviation (see [START_REF] Dembo | Large deviations techniques and applications[END_REF]). We compute the Laplace transform of X T n Λ n X n :

E P f θ 0 e λX T n Kn( 1 f θ )Xn = 1 ( √ 2π) mn det(K n (f θ 0 )) e 1 2 X T n (Kn(fθ 0 )) -1 -2λKn( 1 f θ ) Xn = 1 det(K n (f θ 0 )) det (K n (f θ 0 )) -1 -2λK n ( 1 f θ ) -1 = 1 det I Gn -2λK n (f θ 0 ) 1 2 K n ( 1 f θ )K n (f θ 0 ) 1 2
.

These last equalities hold as soon as

I Gn -2λK n (f θ 0 ) 1 2 K n ( 1 f θ )K n (f θ 0 ) 1 2
is positive. This is true whenever λ ≤ 0 or small enough. Now, for λ ≤ 0, define

φ n (λ) := 1 m n log ( E P f θ 0 e λX T n Kn( 1 f θ )Xn , This function verifies φ n (λ) = - 1 2m n log det I Gn -2λK n (f θ 0 ) 1 2 K n ( 1 f θ )K n (f θ 0 ) 1 2
.

Define also

φ(λ) = lim n φ n (λ),
We get, using Corollary 6.1,

φ(λ) = - 1 2 log 1 -2λ f θ 0 f θ .
We can also compute

φ ′′ (λ) = 2( f θ 0 f θ ) 2 (1 -2λ f θ 0 f θ ) 2 dµ > 0.
As very usual, we define the convex conjugate of φ by

φ * (t) := sup λ∈R - [λt -φ(λ)] , t ∈ R.
As soon as φ is strictly convex, φ * (t) > φ(0) = 0, for any t = φ ′ (0) = f g dµ.

We can now write, for λ ≤ 0,

1 m n log(P( 1 m n X T n Λ n X n ≥ t)) = 1 m n log(P(e λX T n ΛnXn ≥ e mnλt )) ≤ 1 m n log e -mnλt + 1 m n log E[e λX T n ΛnXn ] ≤ -λt + φ n (λ).
Then we get, ∀t > f g dµ, lim sup

n 1 m n log(P( 1 m n X T n Λ n X n ≥ t)) ≤ -λt + φ(λ)
So that, taking the infimum on λ, we get lim sup

n 1 m n log(P( 1 m n X T n Λ n X n ≥ t)) ≤ -φ * (t) < 0
We can obtain the same bound for t < f g dµ. By Borel-Cantelli theorem, we get the P f θ 0 -almost sure convergence of 1 mn X T n Λ n X n to f g dµ. To prove the same convergence with Λ n = (K n (f θ )) -1 , we have to show that the difference between the spectral empirical measure of K n (f θ 0 )

1 2 K n ( 1 f θ )K n (f θ 0 ) 1 2 and K n (f θ 0 ) 1 2 K n (f θ ) -1 K n (f θ 0 ) 1 2
converges weakly to zero. It is sufficient to control the convergence of every moment, because these two last measures both have compact support.

For this, we make use of the Schatten norms. For any A, B matrices of M mn (R), we define

A Sch,p = s k (A) p 1 p ,
where s k (A) are the singular values of A.

Note that

|T r(AB)| ≤ AB Sch,1 ≤ A Sch,1 B Sch,∞ .
Recall that since f θ ∈ F ρ , we have e -ρ ≤ f θ ≤ e ρ . Hence, for any p ≥ 1,

1 m n Tr K p n ( 1 f θ ) K p n (f θ 0 ) -K -p n (f θ )K p n (f θ 0 ) ≤ 1 m n K n (f θ ) -p K p n (f θ 0 ) Sch,∞ K p n ( 1 θ )K p n (f θ ) -I Gn Sch,1 ≤ δ n m n
e 2ρp e -2ρp α(f θ ) 2p α(

1 f θ ) 2p → n→∞ 0.
To obtain the same bound with Λ n = Q n ( 1 f θ ), we have to prove that the difference between the spectral empirical measures of K n (f θ 0 )

1 2 K n ( 1 f θ )K n (f θ 0 ) 1 2 and K n (f θ 0 ) 1 2 Q n ( 1 f θ )K n (f θ 0 ) 1 2
converge weakly to zero. This last assertion is a direct consequence of Lemma 6.2. So, we get 1 m n X T n Λ n X n → f θ 0 f θ , P f θ 0a.s. Recall that we aim at proving that, P f θ 0 -a.s., the sequence of functions (l n (f θ 0 , X n )l n (f θ , X n )) n∈N is equicontinuous on {f θ , θ ∈ Θ}, and that this property also holds for ln , ln and l

|λ i | = K n (f θ ′ ) -1 (K n (f θ ′ ) -K n (f θ )) 2,op ≤ e ρ f θ ′ -f θ ∞ .
So that, to prove the equicontinuity, we may assume that θ is close enough to θ ′ to ensure that sup i=1,••• ,n |λ i | ≤ Furthermore, the sequence ( log(f θ )dµ) n∈N is also equicontinuous since, using a Taylor formula, log(f θ ′ )dµlog(f θ )dµ ≤ e ρ f θ ′f θ ∞ . Now we tackle the equicontinuity of the sequences

X T n K n (f θ ) -1 X n n∈N , X T n K n ( 1 f θ )X n n∈N and X T n Q n ( 1 f θ )X n n∈N .
Notice first that, for any matrix B ∈ M n (R),

1 m n X T n BX n ≤ 1 m n B 2,op X T n X n .
It is thus sufficient to prove the equicontinuity of the sequences (K n (f θ ) -1 ) n∈N ,

(K n ( 1 f θ )) n∈N and (Q n (f θ ) -1 ) n∈N ,
for the norm . 2,op Note that

K n ( 1 f θ ′ ) -K n ( 1 f θ ) 2,op ≤ 1 f θ ′ - 1 f θ ∞ ≤ e 2ρ f θ ′ -f θ ∞ .
Then,

(K n (f θ ′ )) -1 -(K n (f θ )) -1 2,op ≤ (K n (f θ ′ )) -1 (K n (f θ )) -1 2,op (K n (f θ ′ )) -(K n (f θ )) 2,op ≤ e 2ρ f θ ′ -f θ ∞ .
Then, recall that, for any symmetric matrix B ∈ M n (R), we have

B 2,op ≤ B ∞,op . Recall also that Q n (f θ ) = B (n) ⊙ K n (f θ ). Denote Q n ( 1 f θ ′ ) -Q n ( 1 f θ ) 2,op ≤ Q n ( 1 f θ ′ ) -Q n ( 1 f θ ) ∞,op ≤ sup i,j=1,•••n B (n) ij K n ( 1 f θ ′ ) -K n ( 1 f θ ) ∞,op ≤ (1 + u n ) 1 f θ ′ - 1 f θ 1,pol
(see Assumption 3.5).

Since the map f θ → 1 f θ is continuous over F ρ , which is compact, we get the uniform equicontinuity of the map f θ → X T n Q n ( 1 f θ )X n (for the norm . 1,pol ). This concludes the proof of Lemma 6.5

Proof. of Lemma 3.3 We aim at proving the asymptotic normality of √ m n (l (u) n ) ′ (θ 0 ). Using the Fourier transform, it is sufficient to prove that The independence of Y k leads to log (E [exp (i

√ m n (Z n -E [Z n ]))]) = - mn k=1 i τ k √ m n + 1 2 log(1 -2i τ k √ m n ) .
The τ k are bounded, thanks to the following inequality:

M n 2,op = t 2 K n (f θ 0 ) 1 2 Q n ( f ′ θ 0 f 2 θ 0 )K n (f θ 0 ) 1 2 2,op ≤ t 2 K n (f θ 0 ) 1 2 2,op Q n ( f ′ θ 0 f 2 θ 0 ) 2,op K n (f θ 0 ) 1 2 2,op ≤ t 2 K n (f θ 0 ) 1 2 2,op Q n ( f ′ θ 0 f 2 θ 0 ) 1,op K n (f θ 0 ) 1 2 2,op ≤ e ρ α(f ′ θ 0 )α(f θ 0 ) 2 (1 + u n ).
The Taylor expansion of log This last convergence is a consequence of Lemmas 6.1 and 6.2. This provides the asymptotic normality of √ m n (l (u) n ) ′ (θ 0 ) and concludes the proof of Lemma 3.3:

√ m n (l (u) n ) ′ (θ 0 ) → n→∞ N (0, 1 2
f ′ θ 0 f θ 0 2 dµ).
Proof. of Lemma 3.4

We aim now at proving the P f θ 0 -a.s. following convergence:

(l (u) n ) ′′ ( θn )

-1 → n→∞ 1 2 (f ′ θ 0 ) 2 f 2 θ 0 dµ -1
We have (l (u) n ) ′′ (θ) = -

Figure 1 :

 1 Figure 1: Possible configurations for couple of vertices

Figure 2 :

 2 Figure 2: Graph G

Figure 3 Figure 4 :

 34 Figure 3: Empirical spectrum

.

  will prove the equicontinuity of the sequence1 m n log det(K n (f θ )) n∈N Let θ, θ ′ ∈ Θ. Denote λ i the eigenvalues of K n (f θ ′ ) -1 (K n (f θ ′ ) -K n (f θ )). Since f θ ∈ F ρ , we have e -ρ ≤ f θ ≤ e ρ .Notice that we have sup i=1,••• ,n

  det(K n (f θ ′ ))log det(K n (f θ )) = 1 m n log det I Gn -K n (f θ 0 ) -1 (K n (f θ ′ ) -K n (f θ )) ≤ 1 m n i∈Gn |log(1 + λ i )| ≤ 1 m n sup i∈Gn |log(1 + λ i )| ≤ 2 log(2) sup i∈Gn |λ i | ≤ 2 log(2)e ρ f θ ′f θ ∞ .

( 1 -|τ k | 3

 13 2 τ k √ mn ) gives log (E [exp (i √ m n (Z n -E [Z n ]))]) = -Since the τ k are bounded the assertion will be proved if we show that


François Beauducel
email: beauducel@ipgp.fr

Sara Bazin

Mendy Bengoubou-Valerius

Marie-Paule Bouin

Alexis Bosson

Christian Anténor-Habazac

Valérie Clouard

Jean-Bernard De Chabalier

Mendy Bengoubou-Valérius

Empirical model for rapid macroseismic intensities prediction in Guadeloupe and Martinique Modèle empirique pour la prédiction rapide des intensités macrosismiques en Guadeloupe et Martinique

Keywords: earthquake, ground motion, macroseismic intensity, Lesser Antilles, peak ground acceleration

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The Lesser Antilles arc is a zone of convergence between the American plate and the Caribbean plate at a rate of about 2 cm/yr [START_REF] Lopez | Is there a northern Lesser Antilles forearc block ?[END_REF]. This movement is absorbed by the subduction of the American plate below the Caribbean plate and deformation of the wedge of the upper plate on a 100-250 km-wide zone, producing an extended system of active crustal faults (Fig. 1 inset, [START_REF] Feuillet | Arc parallel extension and localization of volcanic complexes in Guadeloupe, Lesser Antilles[END_REF]). It results in a high seismicity level (about 1000 detected events per year) located on the subduction interface and within the slab with hypocentral depths ranging from 10 km up to 220 km, and within the deformed Caribbean plate with shallow crustal seismicity from 2 km up to 15-20 km in depth. Very shallow earthquakes occurring below/or very close to Guadeloupe archipelago islands can be felt sometimes with magnitude less than 2.0.

Since the french volcanological and seismological observatories (OVSG and OVSM) located in the Lesser Antilles are maintaining operational real-time seismic networks, they are responsible for detecting and informing local authorities and public of any felt earthquake occurrence and main event characteristics : location (epicenter and depth), type (tectonic or volcanic), magnitude, and maximum reported intensity in Guadeloupe and Martinique islands. Location and magnitude calculation are determined in a systematic way, using hand-picked phase arrivals and hypocenter inversion, and are available wi- and rock (triangles) site conditions (see [START_REF] Bengoubou-Valérius | CDSA : A New Seismological Data Center for the French Lesser Antilles[END_REF] for further details). Epicenters of some recent instrumental earthquakes are shown as black stars with dates (see text for details). Inset : Tectonic context of Lesser Antilles. CA = Caribbean plate, NA = North American plate, SA = South American plate. Black thick line with triangles : accretionary prism frontal thrust. Black lines : crustal faults from [START_REF] Feuillet | Arc parallel extension and localization of volcanic complexes in Guadeloupe, Lesser Antilles[END_REF]. NA Euler vector of [START_REF] Lopez | Is there a northern Lesser Antilles forearc block ?[END_REF]. Bathymetry data are 500 m contour lines [START_REF] Smith | Global sea floor topography from satellite altimetry and ship depth soundings[END_REF].

Figure 1: Réseau permanent d'accéléromètres du Rap (numériseurs Titan Agecodagis, voir http ://www-rap.obs.ujf-grenoble.fr/) en Guadeloupe : conditions de site type sol (carrés) ou rocher (triangles) (voir (Bengoubou-Valérius et al., 2008) pour plus de détails). Les épicentres de quelques séismes instrumentaux sont indiqués par des étoiles avec dates (voir texte). Encart : context tectonique des Petites Antilles. CA = plaque Caraïbe, NA = plaque Nord-Américaine, SA = plaque Sud-Américaine. Courbe noir avec triangles : fosse frontale du prisme d'accrétion. Lignes noires : système de failles crustales d'après [START_REF] Feuillet | Arc parallel extension and localization of volcanic complexes in Guadeloupe, Lesser Antilles[END_REF]. Vecteur d'Euler pour la plaque NA d'après [START_REF] Lopez | Is there a northern Lesser Antilles forearc block ?[END_REF]. Données de bathymetrie courbes de niveaux 500 m [START_REF] Smith | Global sea floor topography from satellite altimetry and ship depth soundings[END_REF].

thin few tens of minutes after an event, thanks to observatory permanent duty. Macroseismic intensities are determined later, as a result of detailed investigations in the field. However, in the case of strong felt earthquake, the first need of local authorities is to get practical information on event location and maximum possible effects in the living areas. If this information can be delivered rapidly, it may be used to evaluate and focus assistance in the most concerned zones. On November 21, 2004, the occurrence of Les Saintes event, M w = 6.3 and thousands of aftershocks in few days [START_REF] Ipgp | Bilan mensuel de l'activité volcanique de la Soufrière de Guadeloupe et de la sismicité régionale[END_REF][START_REF] Bertil | Séisme des Saintes : rapport de synthèse[END_REF]Beauducel et al., 2005a;[START_REF] Courboulex | Ground-motion simulations of the 2004 mw 6.4 les saintes, guadeloupe, earthquake using ten smaller events[END_REF][START_REF] Bazin | The 2004-2005 les saintes (french west indies) seismic aftershock sequence observed with ocean bottom seismometers[END_REF] offered an exceptional new strong-motion database thanks to the french permanent accelerometric network [START_REF] Pequegnat | The French accelerometric network (RAP) and national data centre (RAP-NDC)[END_REF](Pequegnat et al., ) installed in 2002(Pequegnat et al., -2004. . Combined with collected testimonies and official intensity estimations for largest events, this provided a unique opportunity to establish a first local ground motion model adapted to the observatory needs.

In this paper, we present the modeling strategy, dataset, results and applications of our empirical model. This work has been previously described in an internal report [START_REF] Beauducel | Loi d'atténuation B-Cube pour l'évaluation rapide des intensités sismiques probables dans l'Archipel de Guadeloupe[END_REF], named B 3 (from initials of the three original authors), and is presently used in Guadeloupe and Martinique seismological observatories to produce automatic reports.

Methodology

Our goal is to produce a predictive model of macroseismic intensities with a final uncertainty of about one intensity level, paying special attention to the maximum values that will be published after each earthquake. To be usable in an operative way, the model must be applicable to a wide range of magnitudes and hypocentral distances, and ideally independently from its tectonic context or depth.

Due to insular configuration of Lesser Antilles, most of epicenters occur offshore : it concerns 95% of M ≥ 2.5 detected events (OVSG-IPGP database). Classical macroseismic intensity models cannot be used because they are based on maximum intensity at epicenter, I 0 (see for instance [START_REF] Pasolini | The attenuation of seismic intensity in Italy, part I : Theoretical and empirical backgrounds[END_REF]; [START_REF] Sorensen | Attenuation of Macroseismic Intensity : A New Relation for the Marmara Sea Region, Northwest Turkey[END_REF]), a meaningless parameter for offshore events. Morever, we do not have sufficient intensity data to well constrain a predictive model for intensities. We have then proceeded by combining, first, a ground motion predictive equation (GMPE) constrained by peak ground accelerations (PGA) local data, and second, applying a forward empirical relation between intensities and accelerations.

Many empirical relations to predict earthquake ground motions have been developed for engineering purposes (see [START_REF] Abrahamson | Overview (of modern attenuation relationships)[END_REF]; [START_REF] Douglas | Earthquake ground motion estimation using strong-motion records : a review of equations for the estimation of peak ground acceleration and response spectral ordinates[END_REF]; [START_REF] Strasser | Sigma : Issues, insights, and challenges[END_REF]; [START_REF] Bommer | On the selection of ground-motion prediction equations for seismic hazard analysis[END_REF] for a short review). Due to the necessary high precision for these specific applications (like building damages study), models are developed using very selected datasets for specific applicability ranges of site conditions, magnitude and depth. Moreover, none of them are valid for magnitudes lower than 4.

Furthermore, in a recent study [START_REF] Douglas | A preliminary investigation of strong-motion data from the french antilles[END_REF], shows that ground motions observed on Guadeloupe and Martinique are poorly estimated by commonly-used GMPE, having smaller and more variable amplitudes than expected.

In this work, we do not intend to produce a new GMPE for the engineering community ; we need a more general model with certainly higher uncertainty, but applicable over a wide range of earthquakes to be used in an operative way. In the following we check results and residuals of our obtained PGA model as an intermediate stage, but in order to validate the choices made to produce automatic reports, we emphasize tests of the final intensity model performance in terms of medians across full range of intensity and distance applicability and beyond.

Intermediate PGA predictive equation

Formulation and dataset

Due to the limited database and model purpose, we use one of the simplest form of GMPEs with only 3 parameters [START_REF] Berge-Thierry | New empirical response spectral attenuation laws for moderate european earthquakes[END_REF] :

log(PGA) = aM + bR -log(R) + c (1)
where PGA is the horizontal acceleration peak (in g), M is the magnitude, R is the hypocentral distance (in km), and a, b, c are constant parameters. This functional form implies many hypothesis. In particular, a radial distribution of ground motion around a point source, neglecting geological heterogeneities, tectonic origin, source extension and radiation pattern. [START_REF] Fukushima | Scaling relations for strong ground motion prediction models with m2 terms[END_REF] also points out that a linear log(D)/M formulation is not verified for magnitudes ≥ 6.5 for which a M 2 term should be necessary. This concerns magnitudes out of our study range, but we will keep in mind that accelerations should be underestimated at long distance for large magnitudes.

To inverse the three parameters, we use seismic data recorded at 14 strong-motion permanent stations in Guadeloupe (see Fig. 1), with mixed site conditions, rock and soil (details about the seismic stations can be found in [START_REF] Bengoubou-Valérius | CDSA : A New Seismological Data Center for the French Lesser Antilles[END_REF]), in the period from November 21 to December 28, 2004. The dataset includes about 400 earthquakes associated to 1430 triggers of 3-component acceleration waveforms. These events correspond to Les Saintes main shock M w = 6.3 and mostly the associated aftershocks, but also some regional events that we voluntarily kept in the database.

Locations and magnitudes come from the seismic catalog of Guadeloupe observatory (OVSG-IPGP). Magnitudes were computed using the classical formula of duration magnitude from [START_REF] Lee | HYPO71 (revised) : a computer program for determining hypocenter, magnitude, and first motion pattern of local earthquakes[END_REF] for events M d ≤ 4.5 [START_REF] Feuillard | Macrosismicité de la Guadeloupe et de la Martinique[END_REF][START_REF] Clément | Compilation et validation du catalogue de sismicité des observatoires IPGP des Antilles[END_REF], and we imposed the moment magnitude from worldwide networks for greater events. This allows to overcome the problem of duration magnitude saturation for magnitude greater than 4.5. The consistency of magnitude scale (M d versus M w ) has been checked by [START_REF] Bengoubou-Valérius | CDSA : A New Seismological Data Center for the French Lesser Antilles[END_REF].

For each event, a value of PGA is calculated as the maximum amplitude of horizontal acceleration signals, using the modulus of a complex vector defined by the two horizontal and orthogonal components x(t) and y(t). The PGA dataset is presented in Fig. 2 and3. Magnitudes range from 1.1 to 6.3, hypocentral distances from 2 to 450 km, and PGA from 16 µg to 0.36 g.

Best model determination and residuals

To calculate the 3 parameters in Eq. ( 1), we minimized a misfit function using the L 2 -norm. Due to the inhomogeneous dataset (magnitudes follow a power-law and there is more shortdistances values), we applied a simple weighting function by multiplying the misfit by the magnitude and a power of the hypocentral distance. This gave more weight for large magnitudes and long distances.

The inversion scheme yields the following parameters : a = 0.61755, b = -0.0030746, and c = -3.3968. It produced an RMS residual on log (PGA) of 0.47 (a factor of 3 in PGA, see Fig. 4). This value is higher than classical published GMPE results (around 0.3, see [START_REF] Strasser | Sigma : Issues, insights, and challenges[END_REF]), and it confirms the observation of [START_REF] Douglas | A preliminary investigation of strong-motion data from the french antilles[END_REF] about abnormal data variability in Lesser Antilles. But interestingly, this factor corresponds to the average ratio between rock and soil conditions in the observed PGA [START_REF] Bengoubou-Valérius | CDSA : A New Seismological Data Center for the French Lesser Antilles[END_REF]. This might also reflects the wide range of magnitudes and distances in a too simple functional form. In order to follow some of the key considerations used to develop GMPEs [START_REF] Bommer | On the selection of ground-motion prediction equations for seismic hazard analysis[END_REF], we checked medians and sigmas of PGA residuals (Fig. 4) : it shows a very consistent distribution in the full magnitude range (from 2 to 6), while we observe a significant PGA underestimation (median around +0.5 so a factor 3 in amplitude) for D < 15 km. Figure 2: Presentation of the seismic events dataset : 1430 triggers recorded by 14 permanent accelerometer stations from 2004.11.21 to 2004.12.28. Data are presented in 3-D view for different X-Y combinations of magnitudes, hypocentral distances (in km) and PGA values (in g). In the 3-D graph the view angle has been chosen to highlight the planar characteristic of the dataset.

Figure 2: Données sismiques utilisées : 1430 déclenchements enregistrés par 14 stations accélérométriques permanentes entre le 2004.11.21 et le 2004.12.28. Les données sont présentées pour différentes combinaisons X-Y de magnitudes, distances hypocentrales (en km) et valeurs de PGA (en g). Sur le graphe 3D l'angle de vue a été ajusté pour mettre en évidence l'aspect planaire du jeu de données. Predicted PGA (g)

Hypocentral distance (km) Eq. ( 1) with found parameters is represented as an abacus in Fig. 5 showing calculated PGA as a function of hypocentral distance (from 3 to 500 km) and magnitudes 1 to 8.

Note that we voluntarily limited the minimum hypocentral distance for each magnitude, as we do not take into account the near fault saturation term. It is reasonable to assume that this minimum hypocentral distance is greater than rupture size. Earthquake magnitude reflects the seismic moment which is proportional to the total displacement averaged over the fault surface [START_REF] Aki | Earthquake mechanism[END_REF][START_REF] Kanamori | The energy release in great earthquakes[END_REF]. Many authors propose a simple formula to express the relationship between magnitude and fault length or rupture area [START_REF] Liebermann | Source dimensions of small earthquakes as determined from the size of the aftershock zone[END_REF][START_REF] Mark | Application of linear statistical models of earthquake magnitude versus fault length in estimating maximum expectable earthquakes[END_REF][START_REF] Wyss | Estimating maximum expectable magnitude of earthquakes from fault dimensions[END_REF][START_REF] Wells | New empirical relationships among magnitude, rupture length, rupture width, rupture area, and surface displacement[END_REF]).

Here we use Wyss's formula [START_REF] Wyss | Estimating maximum expectable magnitude of earthquakes from fault dimensions[END_REF] :

M = log(A) + 4.15 (2)
where M is the magnitude and A the rupture surface. We decide to restrict the attenuation law of Eq. (1) to the domain R > L, where L ≈ A 1 2 is an estimation of the fault characteristic size.

Examples of predicted and observed PGA

Fig. 6 shows representative events with observed PGA compared to our model prediction. We do not limit examples to the events from the dataset which reflects the previous residual analysis (Fig. 4), but present events in the period 2004 to 2007 with various depths, in crustal or subduction context, and for which sufficient triggers were available. As seen in Fig. 6, most of PGA values are predicted within the model uncertainty. Medians of log(PGA) residuals are equal to +0.15, +0.28, +0.10, +0.19, +0.24, and -0.01 for Fig. 6a to f events, respectively. We denote, for these 6 particular examples, a light tendency for PGA underestimation, which seems independent from magnitude. This is consistent with Fig. 4 residual analysis. The only significant PGA misfit appears for one soil condition station in the near field (≈ 15 km) for Les Saintes aftershocks (Fig. 6b andc), that is systematically underestimated by a factor of about 10.

We also compare these results with two published GMPE adapted to shallow crustal events : [START_REF] Sadigh | Attenuation relationships for shallow crustal earthquakes based on california strong motion data[END_REF] and [START_REF] Ambraseys | The prediction of earthquake peak ground acceleration in europe[END_REF]. [START_REF] Sadigh | Attenuation relationships for shallow crustal earthquakes based on california strong motion data[END_REF] model is very similar to our PGA model for magnitudes ≥ 5.0 (Fig. 6a, c andf) but has poor fitting for lower magnitudes (Fig. 6b, d ande) with a systematic overestimation. [START_REF] Ambraseys | The prediction of earthquake peak ground acceleration in europe[END_REF] model has a global poor fitting with overestimation of PGA, particularly for M < 5.0.

Macroseismic intensities

Formulation

Although we know that spectral frequency content of ground acceleration and peak velocity have important implications on the building damages, establishing a direct relation between a single PGA value and macroseismic intensity has proved its efficiency in many cases [START_REF] Murphy | The correlation of peak ground acceleration amplitude with seismic intensity and other physical parameters[END_REF][START_REF] Chiaruttini | The correlation of peak ground horizontal acceleration with magnitude, distance, and seismic intensity for Friuli and Ancona, Italy, and the Alpide belt[END_REF][START_REF] Margottini | Intensity versus ground motion : a new approach using italian data[END_REF][START_REF] Wald | Relationships between peak ground acceleration, peak ground velocity, and modified Mercalli intensity in California[END_REF]. For the Lesser Antilles, we follow the suggestion by [START_REF] Feuillard | Macrosismicité de la Guadeloupe et de la Martinique[END_REF] who studied the historical and instrumental seismicity using the simple empirical relation of [START_REF] Gutenberg | Earthquake magnitude, intensity, energy, and acceleration[END_REF] :

I = 3 log(PGA) + 3/2 ( 3
)
where I is the mean intensity (MSK scale ), PGA is maximum acceleration (in cm.s -2 ≈ mg). Combining equations (1), ( 2) and (3) made the final empirical model formulation (hereafter called B 3 prediction) :

I = 1.85265M -0.0092238R -3 log(R) + 0.3096 R > 10 M-4.15 2 (4)
Note that following MSK scale, intensity must be an integer value. In this paper, we decided arbitrarily to round I to the nearest and smallest integer (e.g., I = 6.0 to 6.9 correspond to intensity of VI).

The resulting model for intensities is presented as right Yaxis in Fig. 5. Following Eq. ( 3), the 0.47 uncertainty on our predicted log(PGA), would imply an uncertainty on I of ±1.4, on which we should add the uncertainty of Eq. (3) itself, which is unknown.

Intensity model residuals

We test our model on a database of 20 recent earthquakes for which we have intensity reports (a total of 254 observations) as well as instrumental magnitudes and hypocenter locations. Events are from various origins with magnitudes 1.6 to 7.4, distances from 4 to 500 km, and observed intensities from I to VIII. This wide panel of event characteristics allows to check our model applicability.

We present in Fig. 7 the intensity residuals versus observed intensity and hypocentral distance. Global standard deviation equal 0.8, with a near zero median value. Residuals are also well distributed over the intensity and distance ranges. Since this database is not statistically sufficient, we will keep uncertainty on intensities deduced from the PGA residuals, i.e., σ = 1.4 corresponding to 68% confidence interval. We also checked that maximum observed intensity for each event is strictly below this probability level (see Fig. 7 solid circles).

Examples of simulated and observed intensities

In Fig. 8, we detail eight examples of the most significant events with observed and predicted intensities (see epicenters in Fig. 1).

Fig. 8a shows the October 10, 1974 "Antigua" earthquake [START_REF] Tomblin | Reconnaissance report of the antigua, west indies, earthquake of october 8, 1974[END_REF][START_REF] Mccann | A large normal-fault earthquake in the overriding wedge of the Lesser Antilles subduction zone : The earthquake of 8 October 1974[END_REF], M s = 7.4, a shallow 30 km-depth with normal-fault mechanism, M s from NEIC USGS, location and MSK intensities from [START_REF] Mccann | A large normal-fault earthquake in the overriding wedge of the Lesser Antilles subduction zone : The earthquake of 8 October 1974[END_REF]. Maximum intensities and distance of observations vary from VIII at 45 km in Antigua to II at 400 km in Virgin Islands. All the observations (9 sites) are within the B 3 prediction uncertainty limits. Median of intensity residuals equals -0.6, sigma is 0.5. This is an unexpected positive result since the model is extrapolated for magnitudes larger than Les saintes (M w = 6.3) ; so this magnitude 7.4 is formally out of our interval of validity. Note also that near-field intensities (at 45 km) seem correctly fitted by the model while this hypocentral distance is very close to our limit defined by Eq. ( 2), which gives L = 42 km.

Fig. 8b shows the March 10, 1976 earthquake, a magnitude M b = 5.9, 56 km-depth on subduction interface north of Guadeloupe (M b from USGS-NEIC, location and MSK intensities from [START_REF] Feuillard | Macrosismicité de la Guadeloupe et de la Martinique[END_REF]). Maximum intensities and distances of observations vary from V in Le Moule (Guadeloupe) at 85-km, to II in Martinique at 150 km distance. Most of the 22 observed intensities are underestimated (median of residuals is +0.4) but still within one sigma uncertainty (RMS equals 0.5).

Fig. 8c shows the January 30, 1982 earthquake, a magnitude M w = 6.0, 63 km-depth on subduction interface north of Guadeloupe (M w and location from Global CMT Project, MSK intensities from [START_REF] Feuillard | Macrosismicité de la Guadeloupe et de la Martinique[END_REF]). Maximum intensities and distances of observations vary from V in various urban districts of Guadeloupe and Antigua at 90 km distance, to II in Barbuda (130 km). Most of the 34 observed intensities are within the B 3 uncertainty limits, with a zero median and RMS on intensity residuals equal to 0.7.

Fig. 8d shows the March 16, 1985 "Redonda" earthquake [START_REF] Girardin | Réseau régional sismique de l'arc des Petites Antilles ; sismicité superficielle (1981-1988[END_REF][START_REF] Feuillet | Active faulting induced by slip partitioning in Montserrat and link with volcanic activity : New insights from the 2009 GWADASEIS marine cruise data[END_REF], a magnitude M w = 6.3, 10 km-depth normal-fault (M w and location from Global CMT Project, MSK intensities from [START_REF] Feuillard | Macrosismicité de la Guadeloupe et de la Martinique[END_REF]). Maximum intensities and distances of observations vary from VI at 30 km in Montserrat to II at 300 km in Martinique. We added a supposed intensity of VII-VIII (light gray dashed rectangle) because important cliff collapses have been observed in the Redonda island, at 10 km-distance from epicenter. All the 23 observed intensities are within the B 3 uncertainty limits (RMS = 0.7) with zero median. Note a very local amplification effect that occurred in the region of Pointe-à-Pitre (Guadeloupe) with an intensity of V to VI at 120 km from the hypocenter.

Fig. 8e shows the November 21, 2004 Les Saintes main shock earthquake of magnitude M w = 6.3, M w from Global CMT Project, location from [START_REF] Bazin | The 2004-2005 les saintes (french west indies) seismic aftershock sequence observed with ocean bottom seismometers[END_REF], EMS98 intensities (see definition in [START_REF] Grunthal | European Macroseismic Scale 1998 (EMS-98)[END_REF]) from an official survey by the BCSF [START_REF] Cara | Séisme des Saintes (Guadeloupe) du 21 novembre 2004, note préliminaire[END_REF]. Maximum intensities and distances of observations vary from VIII at 20 km in Les Saintes to IV at 140 km in Martinique, and correspond to detailed studies carried on by BCSF in 33 different urban districts. All the 29 observed intensities are within the B 3 uncertainty limits (RMS = 0.6 , median = -0.9).

Fig. 8f shows the largest Les Saintes aftershock, on February 14, 2005 of magnitude M w = 5.8, located south of Terre-de- Fig. 8h shows the November 29, 2007 Martinique intermediate-depth (152 km) intraslab earthquake of magnitude M w = 7.4, M w and location from [START_REF] Bouin | Source parameters of the Mw 7.4 Martinique intermediate-depth earthquake (Leeward islands) of November 29, 2007 : A slab-pull event with horizontal fault plane[END_REF] and Global CMT Project, with EMS98 intensities from an official survey by the BCSF [START_REF] Schlupp | Séisme de Martinique du 29 novembre 2007[END_REF]. Maximum intensities and distances of observations vary from VII at 150 km in Martinique to II at 400 km in St-Barthelemy, and correspond to detailed studies carried on by BCSF in 70 different urban districts in Guadeloupe and Martinique, plus other islands reports. Most of the 74 observed intensities are within the B 3 uncertainty limits (RMS = 0.83 , median = -0.1), but we note three underestimated intensities at long distances : V in Saint-Vincent (250 km) and Trinidad (500 km), and IV in Anguilla (443 km). This may be due to local site amplifications because of low frequency content of the seismic waves.

These eight examples confirm that B 3 model seems able to predict average intensities within a global residual of σ = 1.4 degree in the MSK scale, for events of magnitudes up to 7.4 in Lesser Antilles context with various hypocentral distances. This value corresponds to 68% of confidence interval and gives a convincing maximum possible intensity even when local site effects are observed.

Automatic intensity report

These good results and the apparent robustness of the B 3 model made us confident for the release of semi-automatic theoretical intensity report at the Guadeloupe and Martinique observatories. For each located event, maximum intensity is computed for all towns of Lesser Antilles islands. If at least one location reaches an intensity of II, it means that the event has been potentially felt and an automatic report is produced, waiting for seismologist validation.

This simulation allows 1) to confirm that inhabitants may have (or not) felt the event when intensity interval varies from II to III in a town ; and 2) to publish immediately and blindly (without any testimonies) the information of a possible felt earthquake when the predicted maximum intensity reaches IV, which means 68% of confidence for an intensity between I-II and IV.

The report (see an example in Fig. 9) includes a synthetic text resuming the date, location and type of event, the maximum intensity prediction value and corresponding town name and distance. To better take into account potential site effects and increase the precision of the result, the average prediction is given together with the upper limit value (I + σ = I + 1.4) for potential site effects, and MSK intensities are indicated in half-unit values, i.e., I = 6.0 to 6.4 is "VI", and I = 6.5 to 6.9 is "VI-VII". The exhaustive list of urban districts for which theoretical intensity reaches at least II is given. Note that it includes all islands in the Lesser Antilles, while our model has been mainly checked with Guadeloupe and Martinique intensities. This may constitutes a future extension of our study.

The report also includes a location map that presents the islands and towns, earthquake epicenter and theoretical isoseist curves using a shaded color map. A detailed table legend explains the MSK scale and corresponding name, color, PGA interval, potential damages and human perception.

Discussion and conclusions

We propose a simple empirical model for macroseismic intensities prediction for observatory operational purpose. The model is based on intermediate PGA model that has been adjusted using a shallow crustal normal-fault sequence of events. The functional form is only 3-parameters dependent which implies many assumptions and simplifications, but makes it also extremely robust with an uncertainty higher than usual GMPE (a factor of 3). This can be explained also by the fact that we do not select specific site conditions in the database, mixing rock and soil stations. The obtained PGA model has strong potential limits and may not be very useful for engineering purposes, but it exhibits a better fit than previous existing GMPE for Lesser Antilles. Its application domain should be limited to crustal events, magnitude range up to 6.3, and distance range up to 100-200 km.

The deduced intensity model is tested on a wider range of magnitudes, distances and source types of earthquakes. We suggest that the B 3 model is able to correctly predict intensities within ±1.4 (1σ), for magnitudes up to 7.4 and hypocentral distance up to 300 km. At longer distances, we observe a clear underestimation of intensities. A major result of our work is that the final equation seems to exhibit a larger applicability range than intermediate PGA predictive equation. In particular, greater magnitudes and other types of earthquakes such as those located in the subduction slab are well modeled within the given uncertainties.

This model is currently used to produce automatic reports in Guadeloupe (since January 2005) and Martinique (since September 2008) observatories in order to anticipate potentially felt events immediately after the location and magnitude calculation. On a total amount of about 10, 000 located events in Guadeloupe, a third has been potentially felt (minimal intensity of II) and has produced an automatic report. Following the observatory convention, only 200 reports were effectively sent as a public communiqué, when the minimum theoretical intensity reached IV or in case of lower intensity (II or III) when immediate testimonies were received from inhabitants.

During more than 5 years of continuous seismic monitoring and thanks to inhabitants testimonies, the B 3 model is daily controlled by observatory team : comparisons between observations and predicted intensities exhibit an average uncertainty less than ±1 unit in the MSK scale.

The reports were also used for seismic hazards awareness and education of the public and local authorities. Particularly, explaining the fundamental difference between magnitude and intensity of an earthquake, the MSK scale, the uncertainty of prediction due to the law's empirical aspect and simplicity, and the potential site condition effects, thus earthquake-resistant construction advices.

Figure 3 :

 3 Figure 3: Histograms of the seismic dataset : number of trigger versus magnitude, hypocentral distance (in km) and PGA value (in g).

Figure 3 :

 3 Figure3: Histogrammes des données : nombre de déclenchements en fonction de la magnitude, de la distance hypocentrale (en km) et du PGA (en g).

Figure 4 :

 4 Figure4: PGA model residuals distribution : difference between observed and calculated log(PGA) versus magnitude and hypocentral distance (km). Horizontal segments with error bars stand for median values and associated standard deviations for each corresponding intervals. Global sigma of PGA residuals equals 0.47 (a factor of 3).

Figure 4 :

 4 Figure 4: Résidus du modèle PGA : différence entre log(PGA) observé et calculé en fonction de la magnitude et de la distance hypocentrale. Les segments horizontaux avec barres d'erreur représentent les valeurs médianes et l'écart-type associé pour chaque intervalle. L'écart-type global des résidus est égal à 0.47 (un facteur 3).

Figure 5 :

 5 Figure5: The B 3 prediction model : simulation of PGA (left Y-axis and log dotted horizontal grid) and MSK intensities (right Y-axis with roman numerals and solid thin horizontal lines) versus hypocentral distance (X-axis) and magnitude (thick solid lines and numbers). Uncertainty interval (σ = factor 3 in PGA) is indicated as gray vertical error bar. MSK is defined in[START_REF] Medvedev | Seismic intensity scale version 1964[END_REF].

Figure 5 :

 5 Figure 5: Modèle B 3 : simulation des PGA (axe des ordonnées à gauche et grille logarithmique horizontale) et intensités MSK (axe des ordonnées de droite avec chiffres romains et lignes horizontales) en fonction de la distance hypocentrale (abscisses) et de la magnitude (lignes épaisses et chiffres en gras). L'intervalle d'incertitude (σ = facteur 3 sur le PGA) est indiqué par une barre d'erreur verticale grise. L'échelle MSK est définie dans Medvedev et al. (1967).

Figure 6 :

 6 Figure 6: PGA prediction (solid thick curve) and uncertainty (solid thin curves) for various magnitudes and depths. The recorded PGA values are presented in squares in the case of soil site condition and in triangles for rock site condition. Are also shown other attenuation laws : Sadigh et al. (1997) (dashed curve) and Ambraseys (1995) (dotted curve). Inset maps show epicenter location for each event (star). a) Les Saintes main shock, b) and c) are aftershocks. d) is a shallow subduction event located east of Guadeloupe. e) is a deep subduction event located north of Guadeloupe. f ) is a shallow crustal event located south of Antigua.

Figure 6 :

 6 Figure 6: Prédiction des PGA (courbe épaisse) et incertitudes (courbes fines) pour différentes magnitudes et profondeurs. Les PGA observés sont présentés par des carrés pour les conditions de site type sol, et par des triangles pour les conditions de site type rocher. Sont aussi indiquées d'autres lois d'atténuation : Sadigh et al. (1997) (courbe tiretée) et Ambraseys (1995) (courbe pointillées). Les cartes en encart indiquent la position de l'épicentre pour chaque événement (étoile). a) Choc principal du séisme des Saintes, b) et c) sont des répliques. d) Séisme superficiel sur le plan de subduction à l'est de la Guadeloupe. e) Séisme profond de subduction au nord de la Guadeloupe. f ) Séisme crustal superficiel au sud d'Antigua.

Figure 7 :

 7 Figure 7: Intensity model residuals distribution for 20 earthquakes (254 macroseismic observations) : difference between observed and calculated intensities versus observed intensity and hypocentral distance (km). Solid black circles indicate the maximum observed intensity for each event. Thick black line = median value (-0.2). Dotted lines = standard deviation (±0.85). Dashed lines = PGA standard deviation equivalent (±1.4).

Figure 7 :

 7 Figure 7: Distribution des résidus du modèle prédictif des intensités pour 20 séismes (254 observations macrosismiques) : différence entre intensité observée et calculée en fonction des intensités observées et de la distance hypocentrale (km). Les disques noirs indiquent l'intensité observée maximum pour chaque événement. Ligne épaisse = valeur médiane (-0.2). Lignes pointillées = écart-type (±0.85). Lignes tiretées = écart-type PGA équivalent (±1.4).

Figure 8 :

 8 Figure 8: Examples of B 3 predicted intensities for 8 instrumental earthquakes : observed macroseismic intensities (MSK) are presented versus hypocentral distance (km). Gray thick curve stands for the predicted intensity for given magnitude, dashed curves indicate uncertainties, solid rectangles represent intensity observations. a), d), e), f ) g) are shallow crustal earthquakes, b) and c) are ≈ 60 km-depth subduction slab interface, h) is an intermediate depth intraslab subduction.

Figure 8 :

 8 Figure8: Exemples des intensités prédites B 3 pour 8 séismes instrumentaux : intensités macrosismiques observées (MSK) en fonction de la distance hypocentrale (km). La courbe épaisse grise réprésente l'intensité prédite pour une magnitude donnée, les courbes tiretées indiquent l'incertitude, les rectangles pleins représentent les intensités observées. a), d), e), f ) g) sont des séismes crustaux, b) et c) sont des séismes à l'interface de subduction à ≈ 60 km de profondeur, h) est un séisme de subduction intraslab de profondeur intermédiaire.

Acknowledgments

FB thanks Pascal Bernard and Nathalie Feuillet for useful discussions, Victor Huerfano and Fabrice Cotton for very constructive comments that helped us to greatly improve the initial manuscript. Acceleration data come from the French national strong motion permanent network Réseau Accélérométrique Permanent (RAP), available at http://www-rap.obs.ujfgrenoble.fr/. Stations have been installed and maintained from 2002 to 2004 thanks to the effort of technicians from the observatory of Guadeloupe (OVSG-IPGP) : Alberto Tarchini † , Christian Lambert, Laurent Mercier, Alejandro Gonzalez, and Thierry Kitou. Authors warmly thank Guadeloupe inhabitants for their collaboration in collecting testimonies. This is an IPGP contribution #3222.

S I M U L A T I O N Rapport préliminaire de séisme concernant la Guadeloupe et Îles du Nord

Observatoire Volcanologique et Sismologique de Guadeloupe -IPGP Le Houëlmont -97113 Gourbeyre -Guadeloupe (FWI) Tél: +590 (0)590 99 11 33 -Fax: +590 (0)590 99 11 34 -infos@ovsg.univ-ag.fr -www.ipgp.fr Gourbeyre, le 21 novembre 2004 10:00 locales Figure 9: Example of observatory communiqué published in case of potential felt event : here an a posteriori simulation of Les Saintes main shock (data from OVSG-IPGP). The B 3 model is used to estimate the probable macroseismic intensities. Are indicated the exhaustive list of towns or neighborough islands where the event could be felt, with maximum intensity in case of site amplification (after [START_REF] Beauducel | Loi d'atténuation B-Cube pour l'évaluation rapide des intensités sismiques probables dans l'Archipel de Guadeloupe[END_REF]). Title : « Preliminary report of earthquake in Guadeloupe and Northern Islands » Main text : « A strong earthquake (magnitude 6.3 on the Richter Scale) was recorded on Sunday, November 21, 2004 at 07:41 (local time) and identified with tectonic origin. The epicenter was located 14 km south of Terre-de-Haut, 10 km-depth (hypocentral distance of about 18 km). This earthquake may have generated, in the closest areas from epicenter, an average ground acceleration of 160 mg (*), corresponding to macroseismic intensity of VIII (severe potential damages). However, depending on the type of soil, the intensities may have reached locally an intensity of IX-X (potential destructions). » Towns list : « Average (and maximum) probable intensities : » Table : « Human perception / Potential damages / Accelerations (mg) / MSK Intensity » Note : « (*) "mg" is a unit of acceleration equal to one thousandth of Earth's gravity » Figure 9: Exemple de communiqué publié par l'observatoire en cas de séisme ressenti : ici une simulation a posteriori du séisme des Saintes (données OVSG-IPGP).

Le modèle B 3 est utilisé pour estimer les intensités moyennes probables. Est mentionnée la liste exhaustive des communes ou îles voisines ayant pu ressentir le séisme, avec les intensités maximales (en cas d'effets de sites) (d'après [START_REF] Beauducel | Loi d'atténuation B-Cube pour l'évaluation rapide des intensités sismiques probables dans l'Archipel de Guadeloupe[END_REF]).