Manaf Ahmed

Véronique Maume-Deschamps

Pierre Ribereau

Céline Vial

Spatial risk measure for Gaussian processes

Keywords: Risk measures, Spatial dependence, Gaussian processes

come

Introduction

A heat wave is a prolonged period where maxima temperatures are unusually high with respect to the usual ones. Most of the times, these heat waves have also a huge spatial component. For example, in 1936, an extremely severe heat wave hits North America. Many states record high temperatures set during this canicule stood until the canicule of 2012. In 2003, a major heat wave hit Europe (cf [START_REF] Chase | Was the 2003 european summer heat wave unusual in a global context?[END_REF], [START_REF] Herrera | A review of the european summer heat wave of 2003[END_REF]), specially France, leaving over 70,000 deaths (around 15,000 only in France). In France, this climatic event was exceptional due to its intensity since some cities registered eight consecutive days with temperature greater than 40 • , but it was also exceptional due to its spatial extent, covering almost all the country. In probability, this means that the underlying spatial process has a strong spatial dependence even at long distance.

On the other hand, a "classical" storm type in the south of France is a cevenol event. These storms are a particular kind of rainfall, hiting usually the Cevennes in France. They are characterized by extremely heavy and localized rainfalls that lead to severe floods. For example, on september 2002, the Gard department was hitted by an exceptional storm. Some locations received more than 700mm of rain in 24h. This event caused the death of 23 persons. Another example, on june 15th 2010 Draguignan was severly flooded (cf [START_REF] Payrastre | TTO2D Team of HYMEX[END_REF], [START_REF] Martin | Les inondations du 15 juin 2010 dans le centre var : réflexion sur un épisode exceptionnel[END_REF]), leaving 27 dead and 1 billion Euros of damages.

In both situations, one of the main characteristic of the event is its spatial dependence: very strong even for large distances for the heat waves and strong at short distances and weak at larger for the cevenol events. When trying to detect the dangerousness of a region using risk measure, the notion of spatial dependence must taken into account.

Risk measures has been widely studied in literature in the univariate setting (i.e. for random variables). The axiomatic formulation of univariate risk measures has been presented in [START_REF] Artzner | Coherent measures of risk[END_REF]. [START_REF] Föllmer | Spatial risk measures and their local specification: The locally lawinvariant case[END_REF] is concerned with financial products on a network, collections of risk measures indexed by the network are considered. The extension of definitions of risk measures to spatial processes may take various forms. As mentioned previously, one of the main characteristics of these phenomenon is the spatial dependence, for example, in the two previous examples (canicule and cevenol event), after normalisation of the marginals, the risk measure should not be the same since the phenomenons are spatially very different. In [START_REF] Keef | Spatial risk assessment for extreme river flows[END_REF], the authors proposed to evaluate the risk on a region by a probability P(S < s) where S is an integrated damage function. In [START_REF]Spatial risk measures and applications to max-stable processes[END_REF] or [START_REF] Koch | Tools and models for the study of some spatial and network risks:application to climate extremes and contagion in france[END_REF] this idea is developed to define spatial risk measures taking into account the spatial dependence. In these works, the sensitivity of the risk measures with respect to spatial dependence and space is studied. In the same spirit as [START_REF] Artzner | Coherent measures of risk[END_REF], the author propose a set of axioms that a risk measure in the spatial context should verify. The author focus on max-stable processes.

We consider similar spatial risk measures well suited for Gaussian processes. Gaussian processes are relevant for some climate models (temperatures e.g.). As an example, in [START_REF] Leblois | Space-time simulation of intermittent rainfall with prescribed advection field: Adaptation of the turning band method[END_REF], a spatial rainfall generator based on a Gaussian process is proposed.

Our main contributions concern the risk measure based on the excess damage function D + X,u = (X -u) + over a threshold u ∈ R + . We calculate the risk measure for Gaussian processes; we study its properties with respect to the parameters of the model (with a focus on the dependence parameter). Moreover, we study the axiomatic properties of a class of risk measures. This paper is organized as follows: in Section 2 we consider quite general spatial risk measures and develop the axiomatic setting of [START_REF] Koch | Tools and models for the study of some spatial and network risks:application to climate extremes and contagion in france[END_REF]. Section 3 is devoted to the study of the risk measure with damage function (X -u) + and Gaussian processes. We propose explicit forms of this risk measure and derive its behavior. We present in Section 4 a simulation study in order to evaluate this spatial risk measure. Finally, we compute the risk measure on the air pollution in Northern Italy model proposed in [START_REF] Bande | Spatio-temporal modelling for pm10 in piemonte[END_REF] in Section 5 and concluding remarks are discussed in Section 6.

Spatial risk measures.

Considering a process X we will define a risk measure associated to X on a region A ⊂ R 2 of the space. It will be a non negative quantity which represents an average damage or cost due to X on A. Throughout the paper, X is a spatial process:

X = (X(s), s ∈ S) with S ⊂ R 2 .
denotes the euclidian norm on R 2 .

2.1. Normalized loss function. A damage function D represents the relationship between the aggregate losses (e.g economic, health) and the environmental (climate) indicator (e.g air pollution, temperature levels), some economic interpretations may be found in [START_REF] Bosello | Estimating a climate change damage function through general equilibrium modeling[END_REF]. More precisely, we will focus on the excess damage function: let u > 0 be fixed threshold for s ∈ S, (2.2) D + X,u (s) = (X(s) -u) + . For example, when considering air pollutants (like in [START_REF] Bande | Spatio-temporal modelling for pm10 in piemonte[END_REF]), u is a regulatory level which is determined by experts.

2.2. Spatial risk measures. As already mentioned, in spatial contexts the spatial dependency is an important characteristic. Considering the risk measure as the expectation of a normalized loss will not take into account the spatial dependency, it is nevertheless useful to quantify the magnitude of risk with respect to different thresholds u.

We shall consider the spatial risk measure composed from two components: the expectation and variance of the normalized loss,

R(A, D X ) = {E[L(A, D X )], Var L(A, D X ) }, (2.3) =: {R 0 (A, D X ), R 1 (A, D X )}
For stationary processes, the expectation component gives informations on the severity of the phenomenon, while the variance component is impacted by the dependence structure.

Let us remark that

(2.4) R 1 (A, D X ) = 1 |A| 2 A×A Cov D X (s), D X (t) dsdt.
2.3. Axiomatic properties of spatial risk measures. In [START_REF] Artzner | Coherent measures of risk[END_REF], [START_REF]Higher moment coherent risk measures[END_REF], [START_REF] Tsanakas | Risk measures and theories of choice[END_REF] and others, axioms and the behavior of univariate risk measures are presented, while [START_REF]Spatial risk measures and applications to max-stable processes[END_REF] provides an axiomatic setting of risk measures in a spatial context.

In this section we will present a set of spatial axiomatic properties describing the behavior of a real valued spatial risk measure R * (A, D). Axioms 1. and 4. below have been introduced in [START_REF]Spatial risk measures and applications to max-stable processes[END_REF], and studied for some max-stable processes.

Definition 2. Let A ⊂ R 2 be a region of the space.

(1) Spatial invariance under translation

Let A + v ⊂ R 2 be the region A translated by a vector v ∈ R 2 . Then for v ∈ R 2 , R * (A + v, D) = R * (A, D). (2) Spatial anti-monotoncity Let A 1 ,A 2 ⊂ R 2 , two regions such that |A 1 | ≤ |A 2 |, then R * (A 2 , D) ≤ R * (A 1 , D). (3) Spatial sub-additivity Let A 1 ,A 2 ⊂ R 2 be two regions disjointed, then R * (A 1 ∪ A 2 , D) ≤ R * (A 1 , D) + R * (A 2 , D). (4) Spatial super sub-additivity Let A 1 ,A 2 ⊂ R 2 be two regions disjointed, then R * (A 1 ∪ A 2 , D) ≤ min i=1,2 [R * (A i , D)]. (5) Spatial homogeneity Let λ > 0 and A ⊂ R 2 then R * (λA, D) = λ k R * (A, D), that is R * is homogenous of order k, where λA is the set {λx, x ∈ A}.
Remark that in [START_REF]Spatial risk measures and applications to max-stable processes[END_REF], the following damage functions are considered for some max-stable processes: D X (s) = 1 {X(s)>u} , D X (s) = X(s) β . The author proves the invariance by translation in this context, he also proves the monotonicity and super sub-additivity in the case where A 1 , A 2 are either disks or squares. We shall prove monotonicity for Gaussian processes and the damage function D X (s) = D + X,u (s) = (X(s) -u) + in the case where A 1 , A 2 are either disks or squares (see Section 3).

Theorem 2.1. Let X be a stationary process and D X be a positive damage function of X. The risk measure R 1 (•, D X ) is invariant by translation and sub-additive.

Proof. The invariance by translation follows directly from the stationarity. On one other hand, consider

A 1 ,A 2 ⊂ R 2 two disjointed regions. R 1 (A 1 ∪ A 2 , D X ) = Var L(A 1 ∪ A 2 , D X ) = 1 (|A 1 | + |A 2 |) 2 |A 1 | 2 R 1 (A 1 , D X ) + |A 2 | 2 R 1 (A 2 , D X ) + 2Cov   A 1 D X (s)ds, A 2 D X (s)ds     . ≤ 1 (|A 1 | + |A 2 |) 2 |A 1 | 2 R 1 (A 1 , D X ) + |A 2 | 2 R 1 (A 2 , D X ) + 2|A 1 ||A 2 | R 1 (A 1 , D X ) R 1 (A 2 , D X ) by using the Cauchy-Schwarz inequality, ≤ R 1 (A 1 , D X ) + R 1 (A 2 , D X ).
So that we have proved the sub-additivity.

Risk measures for spatial Gaussian processes.

We consider the excess damage function D + X,u = (X -u) + , and X an isotropic standard spatial Gaussian process on S ⊂ R 2 , with auto-correlation function ρ. Then for the fixed threshold u and A ⊂ S, we have

(3.1) L(A, D + X,u ) = 1 |A| A X(s) -u + ds, and R 1 (A, D + X,u ) = Var(L(A, D + X,u )).
In what follows, ϕ is the density of the univariate standard normal distribution, Φ is the survival function of the standard normal distribution, u, v, w is the total probability of a truncated bivariate standard normal distribution with correlation w, that is

(3.2) u, v, w = 1 2π(1 -w 2 ) 1/2 ∞ u ∞ v e -1 2(1-w 2 ) [x 2 -2wxy+y 2 ] dxdy.
In this Section, we first give explicit forms for the risk measure, then we will study the behavior of R 1 (λA, D + X,u ) with respect to λ.

3.1.

Explicit forms for R(A, D + X,u ). We are interested in the explicit calculation of the expectation and variance of L(A, D + X,u ).

Proposition 3.1. Consider X := {X(s), s ∈ S} an isotropic standard Gaussian process with auto-correlation function ρ. Let u ∈ R + be a fixed threshold. We have:

(3.3) R 0 (A, D + X,u ) = ϕ(u) -uΦ(u), and
(3.4) R 1 (A, D + X,u ) = |A| 2 A×A G(τ s,t , u) dsdt,
with τ s,t = s -t and for any h, s ∈ S (3.5)

G(h, u) :=Cov D + X,u (s), D + X,u (s + h) ; G(h, u) = ρ(h) + u 2 u, u, ρ(h) -2uϕ(u)Φ u(1 -ρ(h)) (1 -ρ 2 (h)) 1/2 + 1 -ρ 2 (h) 1/2 ϕ u (1 + ρ(h)) 1/2 2 -ϕ(u) -uΦ(u) 2 .
Proof. Let X be an isotropic standard Gaussian process and u ∈ R + ,

E L(A, D + X,u ) = 1 |A| A E (X(s) -u + ds = 1 |A| A ∞ u xϕ(x)dx -u ∞ u ϕ(x)dx ds = 1 |A| A (ϕ(u) -uΦ(u))ds =ϕ(u) -uΦ(u). (3.6)
On one other hand,

Var L(A, D + X,u ) = 1 |A| 2 A×A Cov D + X,u (s), D + X,u (t) dsdt.
We calculate Cov D + X,u (s), D + X,u (t) by using the results from [START_REF] Rosenbaum | Moments of a truncated bivariate normal distribution[END_REF] on moments m 10 , m 11 of truncated bivariate normal distributions. Let f X 1 ,X 2 be the density function of the random vector (X 1 , X 2 ).

E D + X,u (s)D + X,u (t) = ∞ u ∞ u xy -2ux + u 2 f X(s),X(t) x, y dxdy = u, u, ρ(τ s,t ) m 11 -2u u, u, ρ(τ s,t ) m 10 + u 2 u, u, ρ(τ s,t ) , (3.7)
with

u, v, ρ m 10 = 1 2π(1 -ρ 2 ) 1/2 ∞ u ∞ v x exp - 1 2(1 -ρ 2 ) x 2 + 2ρxy + y 2 dxdy, = ϕ(u)Φ v -ρu (1 -ρ 2 ) 1/2 + ρϕ(v)Φ u -ρv (1 -ρ 2 ) 1/2 ; and u, v, ρ m 11 = 1 2π(1 -ρ 2 ) 1/2 ∞ u ∞ v xy exp - 1 2(1 -ρ 2 ) x 2 + 2ρxy + y 2 dxdy, = ρ u, v, ρ + ρuϕ(u)Φ v -ρu (1 -ρ 2 ) 1/2 + ρvϕ(v)Φ u -ρv (1 -ρ 2 ) 1/2 + (1 -ρ 2 ) 1/2 √ 2π ϕ (u 2 -2ρuv + v 2 ) 1/2 (1 -ρ 2 ) 1/2 . For v = u, we have, u, u, ρ m 10 = (1 + ρ)ϕ(u)Φ u(1 -ρ) (1 -ρ 2 ) 1/2 and u, u, ρ m 11 = ρ u, u, ρ +2ρuϕ(u)Φ u(1 -ρ) (1 -ρ 2 ) 1/2 + (1 -ρ 2 ) 1/2 √ 2π ϕ (2u 2 (1 -ρ)) 1/2 (1 -ρ 2 ) 1/2 .
Finally we get,

E D + X,u (s)D + X,u (t) = (1 -ρ 2 (τ s,t )) 1/2 √ 2π ϕ (2u 2 (1 -ρ(τ s,t ))) 1/2 (1 -ρ 2 (τ s,t )) 1/2 +2uρ(τ s,t )ϕ(u)Φ u(1 -ρ(τ s,t )) (1 -ρ 2 (τ s,t )) 1/2 +ρ(τ s,t ) u, u, ρ(τ s,t ) -2u(1 + ρ(τ s,t ))ϕ(u)Φ u(1 -ρ(τ s,t )) (1 -ρ 2 (τ s,t )) 1/2 +u 2 u, u, ρ(τ s,t ) = u, u, ρ(τ s,t ) ρ(τ s,t ) + u 2 -2uϕ(u)Φ u(1 -ρ(τ s,t )) (1 -ρ 2 (τ s,t )) 1/2 + 1 -ρ 2 (τ s,t ) 1/2 ϕ u (1 + ρ(τ s,t )) 1/2 2 .
The result follows.

Corollary 3.2. Let Y := {Y (s), s ∈ S} be an isotropic Gaussian process with mean µ and variance σ 2 . Let X = Y -µ σ an isotropic and standard Gaussian process. The spatial risk measure R(A, D + Y,u ) statisfies

(3.8) R(A, D + Y,u ) = σE[L(A, D + X,u 0 )], σ 2 Var L(A, D + X,u 0 ) , with u 0 = (u -µ)/σ.
Proof. From the definition of D + Y,u , we have:

E[L(A, D + Y,u )] = 1 |A| A E(Y (s) -u) + ds = 1 |A| A E(µ + σX(s) -u) + ds = σ |A| A E X(s) -( u -µ σ ) + ds =σE[L(A, D + X,u 0 )].
(3.9)

On one other hand,

Var L(A, D + Y,u ) = 1 |A| 2 A×A E D + Y,u (s)D + Y,u (t) -E D + Y,u (s) E D + Y,u (t) dsdt = 1 |A| 2 A×A E (Y (s) -u) + (Y (t) -u) + -E (Y (s) -u) + E (Y (t) -u) + dsdt = 1 |A| 2 A×A σ 2 E (X(s) -u 0 ) + (X(t) -u 0 ) + -σ 2 E (X(s) -u 0 ) + E (X(t) -u 0 ) + dsdt = σ 2 |A| 2 A×A E D + Xs,u 0 (s)D + Xt,u 0 (s) -E D + Xs,u 0 (s) E D + Xt,u 0 (s) ds.dt Therefore, (3.10) Var L(A, D + Y,u ) = σ 2 Var L(A, D + X,u 0 ) .
Corollary 3.2 implies that without loss generality, we may calculate the risk measure for an isotropic standard Gaussian process, expressions for an isotropic non standard Gaussian process will follow. Furthermore from these results we can see that R 0 (A, D + Y,u ) does not depend on the region A but only on the characteristics of the underlying Gaussian process. Then in the following study of the risk measure we will focus on the component R 1 (A, D + Y,u ). The following Theorem is useful to compute the risk measure because it reduces to a one dimension integration. Theorem 3.3. Let X := {X(s), s ∈ S} be an isotropic standard Gaussian process. If the region A is either a disk or a square, the expression Var L(A, D + X,u ) reduces to a one dimensional integration. When A is a disk of radius R

(3.11) Var L(A, D + X,u ) = 2R h=0 G(h, u)f disk (h, R)dh,
where

(3.12) f disk (h, R) = 2h R 2 2 π arccos h 2R - h πR 1 - h 2 4R 2 ,
and G is defined in Equation (3.5).

When A is a square of side R

(3.13) Var L(A, D + X,u ) = √ 2R h=0 G(h, u)f square (h, R)dh,
where, for h ∈ [0, R]

f square (h, R) = 2πh R 2 - 8h 2 R 3 + 2h 3 R 4 and for h ∈ [R, √ 2R], (3.14) f square (h, R) = 2h R 2 -2-b+3 √ b -1+ b + 1 √ b -1 +2arcsin( 2 -b b )- 4 b 1 -(2-b) 2 b 2
,

where b = h 2 R 2 .
Proof. The strategy of proof is the one adopted in [START_REF]Spatial risk measures and applications to max-stable processes[END_REF] for some max-stable processes.

Let S and T be two independent random variables uniformly distributed on A. For any function γ defined on R + , we have

E γ(||S -T ||) = 1 |A| 2 A×A γ(||s -t||)dsdt.
Using [START_REF] Moltchanov | Distance distributions in random networks[END_REF] if A is a square of side R,

(3.15) E γ(||S -T ||) = √ 2R h=0 γ(h)f square (h, R)dh,
with f square given by Equation (3.14). If A is a disk of radius R then

(3.16) E γ(||S -T ||) = 2R h=0 γ(h)f disk (h, R)dh.
Moreover, by (3.4) In what follows, we write our results for square regions A, but the results hold for disks as well.

Var L(A, D + X,u ) = E G(||S -T ||, u) .

3.2.

Behavior of R 1 (λA, D + X,u ) with respect to λ. The following expression of R 1 (λA, D + X,u ) is a keystone to understand its behavior. Lemma 3.4. Let λ ≥ 0 and A be a square of side R, then

(3.17) R 1 (λA, D + X,u ) = √ 2R h=0 f square (h, R)G(λh, u)dh.
Proof. Theorem 3.3 gives:

R 1 (λA, D + X,u ) = Var L(λA, D + X,u ) = √ 2λR h=0 f square (h, λR)G(h, u)dh. = √ 2R h=0 f square (λh, λR)G(λh, u)λdh.
Remark that f square (λh, λR) = λ -1 f square (h, R). Thus,

R 1 (λA, D + X,u ) = √ 2R h=0 f square (h, R)G(λh, u)dh.
The same calculations would give the same result if A is a disk of radius R (by replacing f square by f disk ).

Lemma 3.4 gives the following two results on the behavior of the mapping λ → R 1 (λA, D + X,u ).

Corollary 3.5. Let X be an isotropic standard Gaussian process on S ⊂ R 2 with auto-correlation function ρ. Let A ⊂ S be either a disk or a square. The mapping λ → R 1 (λA, D + X,u ) is non-increasing if and only if h → ρ(h), h > 0 is non-increasing and non-negative .

Proof. It suffices to remark that by its definition, for any h > 0, the function λ → G(λh, u) is non-increasing provided the auto-correlation function is non-negative and non-increasing. Corollary 3.6. Let X := {X(s), s ∈ S} be an isotropic standard Gaussian process with auto-correlation function satisfying decreasing to 0 as h goes to infinity. Then, for A either a disk or a square, we have

(3.18) lim λ→∞ R 1 (λA, D + X,u ) = 0.
Proof. Let A be a square of side R,

(3.19) R 1 (λA, D + X,u ) = √ 2R h=0 f square (h, R)G(λh, u)dh,
the monotonic convergence theorem gives:

(3.20) lim λ→∞ R 1 (λA, D + X,u ) = √ 2R h=0 f square (h, R) lim λ→∞ G(λh, u)dh.
Since ρ(h) goes to 0 as h goes to infinity

lim λ→∞ G(λh, u) = u 2 u, u, 0 -u 2 Φ 2 (u)
and the result follows.

We finish this section with the remark that Lemma 3.4 implies the antimonotonicity for regions A 1 , A 2 which are either a disk or a square. Property 3.7. Let X := {X(s), s ∈ S} be an isotropic standard Gaussian process with non-negative and non-increasing auto-correlation function, let A 1 , A 2 be either squares or disks such that

|A 1 | ≤ |A 2 | then R 1 (λA 2 , D + X,u ) ≤ R 1 (λA 1 , D + X,u ). Proof.
Let us do the proof in the square case. By invariance by translation, we may assume A 1 = λA 2 for some λ ≤ 1. Equation (3.17) gives the result.

We now perform some simulation study for various shapes of auto-correlation functions.

Simulation study

In this section, we study the behavior of the proposed spatial risk measure R(A, D + X,u ), through some simulations.

4.1. Analysis of G(h, u) and R 1 (λA, D + X,u ). We begin this simulation section with the study of the covariance damage function G which plays a central role in the behavior of R(A, D + X,u ). Following [START_REF] Abrahamsen | A review of gaussian random fields and correlation functions[END_REF], we consider five Gaussian models depending on the choice of the correlation structure, more precisely,

(1) Spherical correlation function:

ρ sph θ (h) = 1 -1.5 h θ + 0.5 h θ 3 1 {h>θ} .
(2) Cubic correlation function :

ρ cub θ (h) = 1 -7 h θ + 35 2 h θ
where Γ is the gamma function, K κ is the modified Bessel function of second kind and order κ > 0, κ is a smoothness parameter and θ is a scaling parameter.

In order to emphasize the dependence of the damage covariance function G to the correlation parameter we will denote it by G θ (h, u) for any triplet (h, u, θ). X,u (•) and D + X,u (• + h) with respect to the distance h, when the correlation length is set to θ = 0.50 and the threshold to u = Φ -1 (0.75), where Φ -1 is the quantile function of the standard normal distribution. It shows that G θ (h, u) tends to 0 as h tends to infinity with different decreasing speed. This is obviously the expected behavior, because the process (D + X,u (s), s ∈ S) is (spatially) asymptotically independent. Whereas, for spherical and cubic correlation functions, G θ (h, u) = 0 as soon as h > θ, which means that the process (D + X,u (s), s ∈ S) is θ-independent (independent at distance larger that θ).

In order to study the behavior of the damage covariance function with respect to θ, we set the threshold u = Φ -1 (0.75) and the distance h = 0.30. In Figure 1.(b) we remark that G θ (h, u) is increasing with θ. Finally, we study the behavior of the damage covariance function with respect to the threshold u = Φ -1 (p), p ∈ [0, 1]. We set θ = 0.50 and h = 0.30. Remark (see Figure 1.(c)) that even if h is small, G θ (h, Φ -1 (p)) goes to zero as p goes to 1, so that it will be difficult to approximate correctly the covariance damage function when u is large. Figure 2. focusses on the behavior of R 1 (λA, D + X,u ) with respect to (λ, u, θ), when A is a square of side R = 1. In order to see the influence of the homothety rate λ, we set u = Φ -1 (0.75) and θ = 0.50.

To tackle the behavior with respect to θ we choose λ = 1 and u = Φ -1 (0.75).

To study the behavior of the variance with respect to the threshold u = Φ -1 (p), p ∈ [0, 1], we set λ = 1 and θ = 0.50. 4.2. Numerical computation. We generated isotropic standard spatial Gaussian processes X on S = R 2 with different non-negative correlation functions (exponential, Gaussian, spherical, cubic and Matérn with κ = 1) for θ = 0.5. The process X is simulated on a (15 × 15) irregular grid with n = 125 locations over A = [0, 1] 2 .We set the threshold u = Φ -1 (p), for p := {0.75, 0.85, 0.95}. This section is devoted to a numerical study of the computation of R 1 (λA, D + X,u ), where A = [0, 1] 2 . We compare the computation of R 1 (λA, D + X,u ) by the one dimensional integration using (3.4) with the intuitive Monte-Carlo computation (M1). The (M1) computation is obtained by generating a m = 1000 sample of X on the grid. That is, (4.1)

L j (A, D + X,u ) = 1 |A| 1 n -1 2 n-1 i=1 (X(s ij ) -u) + j = 1, ..., m. (4.2) R M 1 0 (λA, D + X,u ) = 1 m m j=1 L j (A, D + X,u ) and (4.3) R M 1 1 (λA, D + X,u ) = Var(L j (A, D + X,u )).
Boxplots in Figure 3 represent the relative errors over 100 (M1) computations with respect to the one dimensional integration. Because exponential, Gaussian and matern correlation models have relatively simple forms, the relative errors are expected to be smaller compared to spherical and cubic ones. For cubic and spherical models, the discontinuity at h = θ induces more instability in the simulations.

Piemonte case study

We terminate this paper with the computation of the risk measure R 1 (λA, D + X,u ) on pollution in Piemonte data. The air pollution is measured by the concentration in P M 10 , particulate matter with an aerodynamic diameter less than 10µm. The observed values of P M 10 are frequently larger than the legal level fixed by the European directive 2008/50/EC (see [START_REF] Cameletti | Spatio-temporal modeling of particulate matter concentration through the spde approach[END_REF] for details). The data has been fitted and analyzed in [START_REF] Bande | Spatio-temporal modelling for pm10 in piemonte[END_REF] The random variable L(A, D + Y,log u ) is the average over the square A of the values of Y that exceed the legal threshold log u. This is a quantity of interest for health public policies. Our study shows that the standard deviation of L(A, D + Y,log u ) is large with respect to its expectation. This means that the dependence structure of the underlying process highly impacts the random variable L(A, D + Y,log u ).

Conclusion

We have proposed a spatial risk measure R(A, D + X,u ). It takes into account the spatial dependence over a region. We showed that some proposed axioms are valid for any stationary processes. Properties such as anti-monotonicity is verified for isotropic Gaussian processes and A a disk or a square (the same result holds for some max-stable processes, see [START_REF]Spatial risk measures and applications to max-stable processes[END_REF]). A simulation study emphasized the behavior of the risk measure with respect to the various parameters. Finally, the computation on pollution data showed the interest of using the variance of L(A, D + Y,log u ) as a spatial risk measure in concrete cases.

Definition 1 .

 1 (Normalized loss function) Consider a damage function D : R d → R + . For any set A ∈ B(R d ) the normalized aggregate loss function on A is (2.1) L(A, D) = 1 |A| A D(s) ds, where |A| stands for the volume (or the Lebesgue measure) of A. The quantity A D(s)ds represents the aggregated loss over the region A. Therefore the function L(A, D) is the proportion of loss on A. In our context, D will be a function of the process X, denoted D X .

Using ( 3 .

 3 [START_REF] Moltchanov | Distance distributions in random networks[END_REF] and(3.16) with the function γ(h) = G(h, u) we obtain the result.

Figure 1 .

 1 Figure 1. Behavior of G θ (h, u) with respect to the threshold u, the correlation length θ and the distance h. Five nonnegative correlation functions (exponential, Gaussian, spherical, cubic and Matérn with κ = 1) have been examined. The graphs (a), (b) and (c) show the behavior of G • (•, •) with respect to: (a) the distance h, when u = Φ -1 (0.75) and θ = 0.50; (b) the correlation length θ, when u = Φ -1 (0.75) and h = 0.30; (c) the threshold u = Φ -1 (p), p ∈ [0, 1], when θ = 0.50 and h = 0.30.

Figure 1 .

 1 Figure 1.(a) shows the behavior the spatial covariance between two damage functions D + X,u (•) and D + X,u (• + h) with respect to the distance h, when the correlation length is set to θ = 0.50 and the threshold to u = Φ -1 (0.75), where Φ -1 is the quantile function of the standard normal distribution. It shows that G θ (h, u) tends to 0 as h tends to infinity with different decreasing speed. This is obviously the expected behavior, because the process (D + X,u (s), s ∈ S) is (spatially) asymptotically independent. Whereas, for spherical and cubic correlation functions, G θ (h, u) = 0 as soon as h > θ, which means that the process (D + X,u (s), s ∈ S) is θ-independent (independent at distance larger that θ).In order to study the behavior of the damage covariance function with respect to θ, we set the threshold u = Φ -1 (0.75) and the distance h = 0.30. In Figure1.(b) we remark that G θ (h, u) is increasing with θ. Finally, we study the behavior of the damage covariance function with respect to the threshold u = Φ -1 (p), p ∈ [0, 1]. We set θ = 0.50 and h = 0.30. Remark (see Figure1.(c)) that even if h is small, G θ (h, Φ -1 (p)) goes to zero as p goes to 1, so that it will be difficult to approximate correctly the covariance damage function when u is large.

Figure 2 .

 2 Figure 2. The behavior of R 1 (λA, D + X,u ) for A = [0, 1] 2 . Exponential, Gaussian, spherical, cubic and Matérn with κ = 1 non-negative correlation functions. The graphs (a), (b) and (c) show the behavior of R 1 (λA, D + X,u ) for a fixed h = 0.30 with respect to : (a) λ, when u = Φ -1 (0.75) and θ = 0.50; (b) θ, when u = Φ -1 (0.75) and λ = 1; (c) u = Φ -1 (p), p ∈ [0, 1], when λ = 1 and θ = 0.50.

Figure 3 .

 3 Figure 3. The boxplots represent the relative errors of R 1 (λA, D + X,u ) between the one dimensional integration computation and the M1 method for different thresholds u = Φ -1 (p), p := {0.75, 0.85, 0.95} and five correlation functions (exponential, Gaussian, spherical, cubic and Matérn with κ = 1) for correlation length θ = 0.5 over A = [0, 1] 2 .

  . The data contains the daily concentration of P M 10 during the winter season 2005-March 2006. The authors considered 24 monitoring stations for estimating the parameters of this model and 10 stations for validation. The log of P M 10 has been fitted on an isotropic Gaussian process with

Figure 4 .

 4 Figure 4. [3]. Locations of the 24 P M 10 monitoring sites (red dots) and 10 validation stations (blue squares) in northern Italy between Alps and Appenises (Piemonte region).

Acknowledgements: This work was supported by the LABEX MILYON (ANR-10-LABX-0070) of Université de Lyon, within the program "Investissements d'Avenir" (ANR-11-IDEX-0007) operated by the French National Research Agency (ANR).


Maxime Baelde

Christophe Biernacki

Raphaël Greff

A MIXTURE MODEL-BASED REAL-TIME AUDIO SOURCES CLASSIFICATION METHOD

Keywords: real-time, audio identification, statistical learning, mixture models, sound classification

Recent research on machine learning focuses on audio source identification in complex environments. They rely on extracting features from audio signals and use machine learning techniques to model the sound classes. However, such techniques are often not optimized for a real-time implementation and in multi-source conditions. We propose a new real-time audio single-source classification method based on a dictionary of sound models (that can be extended to a multi-source setting). The sound spectrums are modeled with mixture models and form a dictionary. The classification is based on a comparison with all the elements of the dictionary by computing likelihoods and the best match is used as a result. We found that this technique outperforms classic methods within a temporal horizon of 0.5s per decision (achieved 6% of errors on a database composed of 50 classes). Future works will focus on the multi-sources classification and reduce the computational load.

INTRODUCTION

Audio source classification is a vast and trendy topic in machine learning, and can be divided into three groups. Music Information Retrieval (MIR) aims at recognizing musical instruments [START_REF] Joder | Temporal Integration for Audio Classification With Application to Musical Instrument Classification[END_REF] or musical genres [START_REF] Kour | Music Genre Classification using MFCC, SVM and BPNN[END_REF] from musics. Automatic Speech Recognition (ASR) aims at detecting and identifying speakers in audio recordings [START_REF] Palaz | Convolutional Neural Networks-based continuous speech recognition using raw speech signal[END_REF]. Environmental Sound Recognition (ESR) aims at recognizing classes of sounds that are neither music nor speech [START_REF] Chachada | Environmental sound recognition: A survey[END_REF]: for instance, airplanes or gunshots.

Typical audio source classification methods include two stages. The first one consists in extracting features from the signals using audio descriptors [START_REF] Peeters | The Timbre Toolbox: extracting audio descriptors from musical signals[END_REF]. Common features are temporal (energy, zero crossing rate,...), spectral (Mel Frequency Cepstral Coefficient, centroid,...) or harmonic (fundamental frequency, inharmonicity,...). This step extracts relevant information from the signal and can be seen as a dimension reduction. The second one uses machine learning algorithm to model the sound classes based on the previous features. Commonly used algorithms are Gaussian Mixture Models (GMM) [START_REF] Clavel | Events Detection for an Audio-Based Surveillance System[END_REF][START_REF] Choi | On using acoustic environment classification for statistical model-based speech enhancement[END_REF], Support Vector Machines (SVM) [START_REF] Lecomte | Abnormal events detection using unsupervised One-Class SVM -Application to audio surveillance and evaluation[END_REF][START_REF] Sameh | On the Use of Time-Frequency Reassignment and SVM-Based Classifier for Audio Surveillance Applications[END_REF], Hidden Markov Models (HMM) [START_REF] Heittola | Audio context recognition using audio event histograms[END_REF][START_REF] Bietti | An online em algorithm in hidden (semi-)Markov models for audio segmentation and clustering[END_REF] or Neural Networks (NN) [START_REF] Biondi | Low Cost Real Time Robust Identification of Impulsive Signals[END_REF][START_REF] Dadula | Neural network classification for detecting abnormal events in a public transport vehicle[END_REF].

More recently, research focuses on neural networks and uses them for features extraction and classification at the same time. These networks often involve convolutional layers and deep architecture (Deep Convolutional Neural Network, DCNN) to model fine details in signals [START_REF] Palaz | Convolutional Neural Networks-based continuous speech recognition using raw speech signal[END_REF][START_REF] Karol | Environmental sound classification with convolutional neural networks[END_REF].

In this study, we develop an audio single-source classification system based on a dictionary of models in a probabilistic framework, using the mixture model theory [START_REF] Mclachlan | Finite Mixture Models[END_REF]. Indeed, the future goal of this research is to identify several sound sources at the same time, which is merely impossible with the current techniques. The use of mixture models allows to deal with mixture of sounds and therefore to identify simultaneous audio sources. Each sound spectrum is modeled by a mixture model and all the models constitute a dictionary. The classification is performed by comparing an unknown signal with the elements of the dictionary by computing likelihoods, aggregating these probabilities and taking the Maximum A Posteriori (MAP).

The rest of the paper is organized as follows. The creation of the dictionary is presented in Section 2. The single-source identification procedure that uses the previous dictionary is detailed in Section 3. The experiments carried out to assess the performance of the method are presented in Section 4. Section 5 presents the results of the experiments. Finally, a brief discussion is presented in Section 6 and concludes the paper.

CREATION OF THE DICTIONARY

The classification algorithm deals with multiple classes of sounds, denoted G i , i = 1, ..., I, and a mono-channel stream. For instance, the classes can be G 1 = airplane, G 2 = gunshot. In each group G i there are multiple sounds, labeled C ij , j = 1, ..., J i . For instance in G 1 , C 11 = airplane 1, and in G 2 , C 23 = gunshot 3. As a typical signal processing step for real-time application, the sounds C ij are splitted into buffers, labeled c ijk , k = 1, ..., K ij of size T samples with a time shift of D samples (see Fig. 1). In a signal notation, c ijk = C ij [kD : kD + T ]. In real world applications, the signal does not necessarily fit correctly into the buffer. To tackle this problem, Gaussian white noise is added at the beginning and the end of each sound. D T Fig. 1. A signal (black plain line) splitted into buffers of size T samples (red empty square) with a time shift of D samples (blue empty dashed square). Gaussian white noise is added at the beginning and the end (gray full block). Magnitude

f f [n] S ijk [n] f [N ] Spectrum Model Components
Fig. 2. Modeling of the spectrums as histograms. The spectrum is represented by the blue dots, the mixture model by a red plain line and its components by green dotted lines.

The Fourier spectrum of each buffer is computed, denoted by s ijk . Only the first N bins in the spectrum are kept, because most of the information is included in the low-frequency content of the signal. As our algorithm is designed to classify multiple sounds present at the same time, the energy spectrum is used. Indeed, when two uncorrelated signals are mixed, their energy spectrums are mixed in the same proportion. The spectrums are normalized so that they sum up to N :

S ijk [n] = N s ijk [n] 2 P N p=1 s ijk [p] 2 , (2.1)
where |z| denotes the modulus of the complex number z. The square values are considered since we want to keep the additivity of the spectrums. These spectrums are considered as histograms and are modeled using a mixture model for binned data [START_REF] Mclachlan | Fitting Mixture Models to Grouped and Truncated Data via the EM Algorithm[END_REF] (see Fig. 2). For the sake of clarity, we will drop the indices (i, j, k) until the end of this section. A mixture model [START_REF] Mclachlan | Finite Mixture Models[END_REF] is a mixture of several probability density functions (pdf), the components, which represents the distribution of a random variable. Here the random variable is the frequency f at which the spectrum is computed, and the pdf is parameterized by a set of parameters ✓ = (⇡ m , µ m , 2 m ) m=1,..,M :

p (f |✓) = M X m=1 ⇡ m N f µ m , 2 m , (2.2)
where p(f |✓) means the probability of f parameterized by ✓, ⇡ m the mixing coefficients (

P m ⇡ m = 1, ⇡ m > 0) and M is the number of components. N f µ m , 2
m is the density of the univariate normal distribution of mean µ m and variance 2 m , taken at point f :

N f µ m , 2 m = 1 p 2⇡ 2 m exp 1 2 ✓ f µ m m ◆ 2 ! .
The frequency f is not known precisely: we know only how many frequencies S[n] fall into the frequency range

⇥ f [n], f[n+1]
⇥ . This framework is close to the one of McLachlan [START_REF] Mclachlan | Fitting Mixture Models to Grouped and Truncated Data via the EM Algorithm[END_REF] except that S[n] is not an integer. Thus, given that the f s are independent and identically distributed (i.i.d.), the probability that one sample falls into

⇥ f [n], f[n + 1] ⇥ is: p (S[n]|✓) = Z f [n+1] f [n] p(f |✓)df ! S[n] . (2.
3)

The Expectation-Maximization (EM) algorithm of Dempster [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF] is used to estimate the set of parameters ✓. This algorithm finds the parameters that maximize the likelihood of the model parameterized by ✓. Given that the S[n]s are i.i.d., the likelihood is:

L(✓) = p(S|✓) = N Y n=1 p (S[n]|✓), (2.4)
where

S = ⇥ S[1], ..., S[N ] ⇤ .
Starting with an initial parameter estimate ✓ (0) , the EM algorithm iterates the two following steps until convergence:

(E) Compute: Q ⇣ ✓ ✓ (p) ⌘ = E h log L(✓) S, ✓ (p) i , (M) Maximize: ✓ (p+1) = argmax ✓ Q ⇣ ✓ ✓ (p)
⌘ .

This model and the related algorithm require the number of components M in the mixture to be specified. For each buffer, the optimal number of components is chosen by minimizing the Bayesian Information Criterion (BIC) of Schwarz [START_REF] Schwarz | Estimating the Dimension of a Model[END_REF]:

BIC(M ) = 2 log L ⇣ b ✓ ⌘ + d log N, (2.5)
where b ✓ denotes the Maximum Likelihood (ML) value of ✓ and d = 3M 1 is the dimension of ✓, with M varying from 1 to 20 in this case. The estimated parameters are gathered into

n b ✓ ijk o , with k = 1, ..., K ij , j = 1, .
.., J i , i = 1, ..., I, which is named the dictionary.

IDENTIFICATION PROCEDURE

Once the dictionary is computed (this corresponds to the learning step), it can be used to identify unknown sounds. Suppose that at a time r there is a new buffer which is transformed into a normalized energy spectrum S r . The identification process consists in finding the correct label G r i by aggregating the following probabilities:

1. Compute the likelihoods of all the models of the dictionary p

⇣ S r c r ijk ⌘ = p ⇣ S r b ✓ ijk

⌘

, for k = 1, ..., K ij , j = 1, ..., J i and i = 1, ..., I using Eq. (2.4).

2. Compute the likelihood of the sound C r ij , for j = 1, ..., J i and i = 1, ..., I:

p(S r |C r ij ) = Kij X k=1 p(S r |c r ijk )p(c r ijk |C r ij ). (3.1)
3. Compute the likelihood the group G r i , for i = 1, ..., I:

p(S r |G r i ) = Ji X j=1 p(S r |C r ij )p(C r ij |G r i ). (3.2)
4. Finally compute the conditional probability of the group G r i for i = 1, ..., I:

p(G r i |S r ) = p(S r |G r i )p(G r i ) P h p(S r |G r h )p(G r h ) . (3.3)
The decision rule uses R buffers (as in [START_REF] Clavel | Events Detection for an Audio-Based Surveillance System[END_REF]) and is computed as follows. For each buffer r and each class the classifier computes a probability p(G r i |S r ). Then these probabilities are multiplied so as to aggregate the results over a larger temporal horizon:

p (G i |S) = R Y r=1 p ⇣ G r i S r ⌘ , (3.4)
where S = ⇥ S 1 , ..., S R ⇤ . The MAP estimate is finally taken on these probabilities: b

G i = argmax Gi p (G i |S).

EXPERIMENTS

Setup of the experiments

Several experiments were carried out to assess the quality of our sound sources classification algorithm. Three databases were considered for these experiments. One database (from A-Volute) was composed of 704 video games sounds divided into 9 classes. The ESC databases [START_REF] Karol | ESC: Dataset for Environmental Sound Classification[END_REF] (50 and 10) were also considered, composed of 2000 environmental sounds from 50 classes for the former and 400 sounds from 10 classes for the latter. Each sound was resampled to 44.1kHz and the mean was subtracted to the signal because it caused instability in the fitting process. Several values were considered for the window size T and the time shift D, mainly T = [512, 1024, 2048] samples and D = 512 samples. The number of kept bins N was set to T /5 which corresponded to approximately 8kHz in our setting. R was set to 10, that corresponds approximately to 0.5s. Each dataset were splitted into a training and a test set for cross-validation procedure. We used 80% of the set for training and the remainder for testing. The division was done in a vfold manner, with v = 5. The recognition rate, defined by the number of buffers correctly labeled over the total number of buffers, was used to assess the performance of the system. In a cross-validation setting, a recognition rate for each fold was computed and the cross-validation recognition rate was simply the mean of the v folds recognition rates.

Three probabilities were introduced in Section 3 and have to be specified:

• The probability of a buffer c ijk given a sound C ij was set to 1: p(c ijk |C ij ) = 1 for k = 1, ..., K ij , j = 1, ..., J i and i = 1, ..., I. • The probability of a sound C ij given a class G i was uniformly distributed over the class

G i : p(C ij |G i ) =
1/#{G i } for j = 1, ..., J i and i = 1, ..., I, where #{A} means the number of elements in A. • The probability of a class G i was set to the ratio between the number of sounds in G i and the total number of sounds:

p(G i ) = #{G i }/ P h #{G h } for i = 1, ..., I.

Comparison with other techniques

We compared our algorithm with other state-of-the-art techniques. The first one was a parametric method inspired from Clavel [START_REF] Clavel | Events Detection for an Audio-Based Surveillance System[END_REF]. Acoustic features were extracted from the signals: energy, 8 MFCC, spectral centroid, spectral spread, plus the first and second derivatives. A Principal Component Analysis step was applied on these descriptors and only the first 13 principal components were kept. The classifier was a standard GMM.

The second one was a non-parametric method: a DCNN. A neural network is a set of neurons (computation units) stacked together in multiple layers that can performs classification or regression. The input of a layer l is denoted by h l 1 , and the convolutional kernel W l . The network computes an output h l by taking an activation function g: h l = g(W l ? h l 1 ), where ? is the convolution operation. Recently, the most popular activation function is the ReLU (Rectified Linear Unit), which is simply g(X) = max (X, 0) element-wise. After a convolutional layer, a max-pooling operation consisting in merging adjacent cells by taking the maximum value is applied. We used the network developed by Piczak [START_REF] Karol | Environmental sound classification with convolutional neural networks[END_REF]. The input was a log-mel spectrogram with its delta, considered as a 2-channel images of size 60 ⇥ 41. The first convolutional layer consisted in 80 filters of size 57 ⇥ 6, with a max-pooling of size 4 ⇥ 3. The second convolutional layer consisted in 80 filters of size 1 ⇥ 3, with a max-pooling of size 1 ⇥ 3. Finally, two fully connected layers composed of 5000 units each and a softmax layer ended the network. For each previous technique, the hyperparameters were adapted so as to have comparable results.

Human listening tests

Because a human score was available for the ESC dataset [START_REF] Karol | ESC: Dataset for Environmental Sound Classification[END_REF], we carried out a listening test on the A-Volute dataset. 21 participants were selected and had to classify 10 sounds from the 9 classes of the dataset. It gave a rough estimate of the human good classification rate on this dataset.

RESULTS

Recognition rate

The experiments were carried out for all the values of T . However, only the results for T = 2048 are shown because they are the best, and are presented in Table 1. The results for the A-Volute's database are the following. The parametric method achieves 73.6% and the non-parametric method 46.6%. Our algorithm outperforms the current methods by achieving a recognition rate of 96.5%. Concerning the ESC-50 database, the parametric method achieves 44%, the nonparametric method 53.2% and our algorithm 94.0%. Finally, on the ESC-10 dataset, the parametric method achieves 73.5%, the non-parametric method 76.0% and our algorithm 96.0%. It is worth noticing that a human can achieve 91.8% on the A-Voute dataset, 94.0% on the ESC-50 database and 95.7% on the ESC-10 database [START_REF] Karol | ESC: Dataset for Environmental Sound Classification[END_REF]. Our algorithm outperforms state-of-the-art methods and humans on these databases.

Complexity and execution time

We evaluate the complexity of the previous algorithms at the identification step. Our algorithm is mainly concerned with multiplication. Indeed, by computing the log-probability, computations are just additions and multiplications. All the operations needed to compute the decision for one buffer can be resumed as a O(#{D}T /5). The parametric method relies mainly on computing exponentials and matrix operations (inversion and determinant), so the complexity is roughly O(d 3 ) where d is the dimension of the feature vector. The neural network uses convolutions which are the more expensive operations when evaluating the decision. The complexity is thus O

⇣ P L l=1 n l 1 (s 1 l s 2 l )n l (m 1 l m 2 l )
⌘ , with L the number of layers, n l 1 the number of input channels for layer l, s i l the dimensions of the input, n l the number of filters of the layer l and m i l the size of the output. Considering the settings of the experiments, the numbers of operations needed for one buffer for our algorithm are O(28• ). We also report the time needed to compute the decision. The machine used a Intel R Core TM i7-5820K CPU @3.30GHz. Using a dictionary composed of approximately 70,000 models, it takes 185ms to the system to infer the decision for one buffer. By using an NVidia GeForce R GTX 960 to compute the multiplications, this computational time can be reduced to 35ms.

DISCUSSION AND CONCLUSION

The aim of this algorithm is to perform real-time audio multisource identification. Up to now, this study examined only the single-source case. As we can see in the results, our technique outperforms standard methods (GMM), more recent algorithms (DCNN) and even humans on both benchmark [START_REF] Karol | ESC: Dataset for Environmental Sound Classification[END_REF] and industrial databases. Despite these good recognition rates, the main advantage of our algorithm is that it can theoretically handle multi-source conditions. The learning algorithm of classic techniques would have to learn every combination of sound classes, which is practically impossible because of the combinatorial and the amount of data required for training. Some research on neural networks begins to study polyphonic identification, as in [START_REF] Parascandolo | Recurrent Neural Networks for Polyphonic Sound Event Detection in Real Life Recordings[END_REF].

The energy spectrum is chosen so as to extend this technique to the multi-source setting. Indeed, we make the assumption that the additivity of two energy spectrums is preserved (uncorrelated sources). Mixture models were employed because they allow a flexible modeling of any signal. Moreover, in a multi-source setting, the signal model is a mixture of mixture, which is well defined in the mixture model framework.

The real-time constraint requires that the signals are not known by advances. Besides, we want the system to have a low latency. This is why such values of T and D were chosen. However, it is not reach yet. Indeed, the larger the dictionary the longer it takes to compute the decision. This is why future work will try to reduce the computation time, for instance by organizing the dictionary in a binary tree.

10 6 )

 106 (A-Volute database), O(63 • 10 6 ) (ESC-10) and O(120 • 10 6 ) (ESC-50), for the parametric method O(2 • 10 3 ) and for the neural network O(14 • 10 6

Table 1 .

 1 The recognition rates in % of the considered methods on the different datasets (A-Volute, ESC-50 and ESC-10).

		A-Volute ESC-50 ESC-10
	Parametric method	73.6	45.5	73.5
	Non-parametric method	46.6	53.2	76.0
	Our algorithm	96.5	94.0	96.0
	Human	91.8	81.3	95.7


domain: info sde

text:
Wilfried Yoro

Mamdouh El Tabach

Taoufik En-Najjary

Azeddine Gati
email: azeddine.gati@orange.comtijani.chahed@telecom-sudparis.eu

Tijani Chahed

Service-oriented Sharing of Energy in Wireless Access Networks Using Shapley Value

Keywords: Energy consumption, Service-oriented, Mobile networks, Shapley value

We investigate in this paper the sharing of energy consumption among service categories in the access of a wireless network. We focus on the fixed component of the energy consumption, which is known to be significantly larger than the loaddependent variable component, and propose its sharing among the service categories based on coalition game concept, the Shapley value. We consider five service categories, two large players: streaming and web browsing, and three smaller ones: download, voice and other minor services, and compare our proposal with two other sharing methods: uniform and proportional which follows the same traffic proportions. Our results, applied on a real dataset extracted from an operational network in Europe, show that our proposal is more fair both towards small services in that it reduces their shares in comparison to the uniform approach and towards larger services as it reduces their shares in comparison with the proportional one. Indeed, our Shapleybased model accommodates both short term network behavior, in which the fixed energy component is independent of the traffic load, and longer term behavior, in which it varies with the load and infrastructure. Uniform sharing accounts only for the short term, and the proportional one only for the longer term.

I. INTRODUCTION

A CCORDING to Cisco [1], overall IP traffic will grow at a compound annual growth rate of 23 percent from 2014 to 2019. In order to face the Internet's growth, Internet providers upgrade their networks so as to keep up or improve the users perceived Quality of Experience. This leads to an increase in power consumption, resulting in turn in two main challenges: economical, as operators margin is decreasing and ecological, in a context aiming at reducing greenhouse gas emissions. Thus the research in green networking has fostered a number of investigations. For instance, Debaillie et al. [START_REF] Debaillie | A flexible and future-proof power model for cellular base stations[END_REF] and Vu et al. [START_REF] Vu | Statistical power consumption assessment of mobile access networks at country level[END_REF] propose models for assessing the energy consumption of the mobile access network. Vu et al. [START_REF] Vu | Energy modeling and optimization of radio access network using stochastic deployment models[END_REF] and Jada et al. [START_REF] Jada | Power efficiency model for mobile access network[END_REF] introduce models for optimizing the energy consumption of the mobile access. Bianzino et al. [START_REF] Bianzino | A survey of green networking research[END_REF] and McLauchlan et al. [START_REF] Mclauchlan | A survey of green energy technology and policy[END_REF] discuss the main trends in the fields of green networking and green energy.

If modeling the energy consumption of network elements is important for optimization sake, investigations are more and more focusing on the energy impact of services on the network. In fact, knowing the energy consumption induced by services should help in the eco-design of applications. Indeed, in order to reduce the energy consumption of the network, a promising way consists in optimizing the energy consumed by the network elements to offer a given service. To do so, one needs a model to assess the energy consumption induced by services so as to measure the energy efficiency gain related to green techniques. For example, the model can be implemented in a Life Cycle Assessment (LCA) software. A LCA software allows users to evaluate, compare and track the environmental performance of products, using notably steps suggested by ISO14040 which specifies the principles and the framework applicable to realize LCA, Ines et al. [START_REF] Ines | AHP multicriteria decision making for ranking life cycle assessment software[END_REF].

To date, the existing models assess mostly the energy consumption for some specific applications (we will report on some of them in section II). We propose in this work a model for sharing the global energy consumption of a mobile access network among the provided services. By "sharing", we mean the sharing of the responsibility of the energy consumption between the different service providers using the network. In effect, traffic in the network originates from different service providers, such as Google, Youtube, voice calls, etc. We group in this paper services into five categories: streaming, web, download, voice and other minor services. We call each service category a player. We propose to share the responsibility of the network's energy consumption between these players. The players' responsibility may represent their contribution in the network's carbon footprint for example.

We focus specifically on the fixed component of energy consumed in a mobile access network and its sharing between different service categories. We consider several service categories representing players of different sizes, large and small, in terms of traffic loads. We decompose the energy into variable versus fixed components and share the former in a manner that is proportional to the traffic proportion of each service category. As of the latter, we propose an approach based on coalition game concept, the Shapley value [START_REF] Shapley | Contributions to the theory of games II. Reading, Massachusetts[END_REF].

We consider two settings: one with a constant network infrastructure and one with an evolving one, over a larger time scale, in terms of additional and/or changing equipment so as to keep up with an increase in the traffic load. In the former setting, energy consumption as well as traffic load shall be considered as constant, whereas in the latter setting, they would both increase.

We study two variants of the sharing models: one where no service category is mandatory and one with mandatory player(s) which reflects the realistic case where some operators may be legally mandated by the state to offer a certain service, such as voice.

Unlike an equal sharing of the fixed energy that is unfair towards players with small traffic load, and a sharing proportional to the traffic load that puts too much weight on the players with large volume of traffic, our Shapley-based model is a trade-off for all the players in that it puts less weight on small players than equal sharing, and less weight on big players than proportional sharing. Our results encourage hence the transport of services with small traffic as well as the introduction of novel ones, and also acknowledge the role of services with large traffic as major drivers for network activity and increased deployment.

The remainder of this paper is organized as follows. In section II, we review some literature related to the assessment of energy consumption per service category. In section III, we describe our Shapley-based models for sharing the fixed component of the energy consumption among service categories. We discuss some implementation issues of the Shapley-based model and how we tackle these issues in section IV. In section V, we run numerical applications, comparing our Shapleybased proposal to uniform as well as proportional sharing of energy, on a real dataset taken from an operational European network transporting three main service categories: streaming, browsing and download, in addition to voice and other minor services. Eventually, section VI concludes the paper.

II. RELATED WORK.

Marquet et al. [START_REF] Marquet | Scientific methodology for telecom services energy consumption and co2 emission assessment including negative and positive impacts[END_REF] investigate the energy consumption of information and communication technology services and CO2 emission at life-cycle of the equipment, including negative and positive impacts: positive impact refers to potential gains due to dematerialization, such as physical transport substitution. Negative impact refers to CO2 emissions notably. The model for assessing the end-to-end energy consumption of a service is based on the notion of consumption rate, i.e., the energy consumed per unit of service, for example kWh/hour/user. However, no distinction is made between fixed and dynamic energy consumption components.

Jalali et al. [START_REF] Jalali | Energy consumption of photo sharing in online social networks[END_REF] discuss an energy model in order to estimate the energy consumption of cloud applications, and is applied to the case of sharing photos on Facebook. The model consists in determining the energy consumed per bit on a device, then multiplied by the traffic volume of the service. Only the load dependent power component is considered, the fixed power component is ignored since it is independent of services.

Preist et al. [START_REF] Preist | Analyzing end-to-end energy consumption for digital services[END_REF] provide a statistical model to calculate the overall energy output required for a digital service, from a Datacenter to the end user, using Monte Carlo analysis. No information is given about the nature of the energy consumption, fixed, dynamic or both.

To date, the investigations in the literature related to modeling the energy consumption of services are based on inputs that are very difficult or even impossible to measure. Marquet et al. [START_REF] Marquet | Scientific methodology for telecom services energy consumption and co2 emission assessment including negative and positive impacts[END_REF] and Jalali et al. [START_REF] Jalali | Energy consumption of photo sharing in online social networks[END_REF] cited above base their model on the energy consumed per transmitted bit of the service by the equipment implied on the path of the service flow. This approach has several limits. Firstly, it allows modeling much more the energy consumed by an application than the one consumed by a service. Secondly, it is quite impossible to measure the energy consumed per bit of service, as most of the time network equipment serve several services simultaneously. In order to overcome this complexity, the authors mostly refer to the power consumption models of the constructors which do not reflect the reality of the field. We propose to base our model on realistic inputs, i.e., the traffic volumes or proportions of the service categories.

Depending on which network segment is considered, one of these component is preponderant over the other. For example, the energy consumption in the core of the network is largely load-dependent because routers energy consumption varies significantly with utilization, while it is largely independent of the load in the RAN because the access is typically underloaded (typically ρ < 50% so as not to exceed some operating load threshold). Fig. 1 shows the power consumption of an operating 4G base station versus its traffic load. These values come from traffic and power measurement probes installed in Orange France's network. At 10% of load, the fixed power consumption represents 91% of the base station total power consumption.

III. MODELING OF THE ENERGY CONSUMPTION SHARING.

A. Description of the system.

A wireless network is composed mainly of three segments, the access, the transport and the core, [START_REF] Vereecken | Power consumption in telecommunication networks: overview and reduction strategies[END_REF], [START_REF] Hinton | Power consumption and energy efficiency in the internet[END_REF], running possibly several technologies: 2G, 3G and 4G, as shown in Fig. 2.

The Radio Access Network is the segment of the mobile network interfacing the end-users and the mobile core network. The GSM EDGE Radio Access Network is composed of the Base Transceiver Station (BTS) and the Base Station Controller (BSC). A BTS implements minimum shift keying modulation for GSM and phase shift keying modulation for EDGE. It provides essentially voice services. The controller is in charge of the radio resource management and implements The UMTS Terrestrial Radio Access Network is composed of the NodeB and the Radio Network Controller (RNC). The NodeB implements hyper phase shift keying modulation and provides voice and data services. The RNC, as with GSM, is in charge of the management of the radio resource and implements Wideband Code Division Multiple Access as resource allocation algorithm.

The eNodeB hosts both the base station and the controller functions in a single equipment, for LTE networks. Orthogonal Frequency Division Multiplex is the modulation technology and Orthogonal Frequency Division Multiple Access the resource allocation algorithm. LTE provides data services only.

Alongside with the network elements on a site, there are the equipment of the Technical Environment composed of the cooling system, the rectifiers and the backup battery.

Traffic and energy measurements are regularly made on the network for management and investigation purposes. These measures will feed our models in the numerical applications.

B. Energy consumption model

Let us first consider a radio access network with only one radio technology (homogeneous network) transporting a set N of N service categories, consuming energy E to be shared among the provided service categories. As stated earlier, the energy consumed by the access equipment is composed of a variable and a fixed components, denoted by E v and E f , respectively. Then we have

E = E v + E f (1)
Denoting by E i the energy consumption induced by service category i, with variable and fixed components E v i and E f i , respectively,

E i = E v i + E f i (2)
We first focus on the variable component of the energy consumption of service category i. Let us denote by v i the traffic volume of service category i.

E v i = ϕ i × E v (3)
where ϕ i is the share of service category i in E v , given by :

ϕ i = v i N k=1 v k (4)
As of the fixed energy component :

E f i = φ i × E f (5)
where φ i is the share of service category i in E f . Unlike the variable component of the energy consumption, that is, as shown above, proportional to the traffic proportions, we propose in this work to determine φ i using the Shapley value.

We begin by giving some introductory material on the Shapley value concept, [START_REF] Shapley | Contributions to the theory of games II. Reading, Massachusetts[END_REF], [START_REF] Myerson | Game Theory, Analysis of Conflict[END_REF]. This mathematical tool has a number of applications in telecommunications, [START_REF] Razzac | Dimensioning and profit sharing in hybrid LTE/DVB systems to offer mobile TV services[END_REF]- [START_REF] Singh | Cooperative profit sharing in coalition-based resource allocation in wireless networks[END_REF].

C. Shapley value 1) Cooperative game theory :

In game theory, a cooperative game (or coalitional game) is a game which allows grouping of players within so-called coalitions, thanks for instance to the possibility of external enforcement of cooperative behavior (e.g., through contract law). These are opposed to non-cooperative games in which there is no possibility to forge alliances. Cooperative games are typically analyzed in the framework of cooperative game theory, which focuses on predicting which coalitions will form, the joint actions that groups take and the resulting collective payoffs. It is opposed to non-cooperative game theory which focuses on predicting individual players actions and payoffs and analyzing Nash equilibriums.

In cooperative game theory, there are two types of solution concepts. The unobjectionable solutions and the equitable solutions. The former guarantees a sharing between the players such that any coalition (grouping of players) cannot increase its gain by leaving the coalition composed of all the players, called the grand coalition. Such solutions include the core. The core is the set of imputations under which no coalition has a value greater than the sum of its members' payoffs. Therefore, no coalition has incentive to leave the grand coalition and receive a larger payoff. An imputation is a solution of the cooperative game that exactly splits the total value of the grand coalition among the players, such that no player receives less than what he could get on his own. A solution is a vector of R N that represents the allocation to each player, with N the number of player. Equitable solutions take into account some consideration of equity between players. Such solutions include the Shapley value. In this paper, we are considering the sharing of the responsibility of the network's energy consumption between the players in a fair way, and hence the use of the Shapley value concept.

2) Shapley value definition and properties :

The Shapley value is a cooperative game concept used for sharing a common gain between the members of a coalition, in a fair way. The fairness of the Shapley value comes from the fact that the payoff of each player in the common gain is proportional to its contribution in the coalition's gain. Let us consider a game ξ(N, V ) with N denoting the number of players and V the characteristic function associating to each coalition of the game a value. A coalition is a set of players that cooperate so as to improve their revenue. The grand coalition consists of all the players. There are N ! possible scenarios of constructing the grand coalition.

Let S σ i denote the largest coalition not containing yet player i in the construction of the grand coalition with regard to scenario σ. We define the incremental cost vector associated to the scenario σ by:

c σ inc = (c σ inc ({1}), • • • , c σ inc ({i}), • • • , c σ inc ({N })) (6) where c σ inc ({i}) = V (S σ i ∪ {i}) -V (S σ i ), i.e.
, the marginal contribution of player i in V (S σ i ∪ {i}). The Shapley value x Shapley is the arithmetic mean of the incremental cost vectors associated to the scenarios of constructing the grand coalition,

x Shapley = 1 N ! σ c σ inc (7)
The Shapley value derived from the following axioms (φ i is the payoff of player i) :

• Efficiency Axiom: i∈N φ i (V ) = V (N )
• Symmetry Axiom: If player i and player j are such that V (S ∪ {i}) = V (S ∪ {j}) for every coalition S not containing player i and player j,

then φ i (V ) = φ j (V ). • Dummy Axiom: If player i is such that V (S) = V (S ∪ {i}) for every coalition S not containing i, then φ i (V ) = 0. • Additivity Axiom: If u and v are characteristic functions, then φ(u + v) = φ(v + u) = φ(u) + φ(v).

D. Game without mandatory players

As stated earlier, the game is characterized by ξ(N, V ), with N the number of players and V the characteristic function. The characteristic function allocates to each coalition a cost, corresponding to a fraction of the fixed energy consumption.

Let S denote a coalition of size s, with s = |S|, |.| the cardinal function. In the sequel, the payoffs of the players and values of the coalitions are normalized by the fixed energy consumption of the network, unless otherwise stated.

V (S) = s k1=1 v k1,S C s N j2=1 s k1=1 v k1,Sj 2 (8)
The value of coalition S is the ratio of its traffic volume versus the traffic volume of all the coalitions of same size as S, whose number is C s N . v k,S is the traffic volume of the k th player of coalition S.

The intuition behind the characteristic function can be seen as follows: the value of a coalition should be a function of the traffic volume and the number of services, as the fixed energy can be shared equally among the players when considering it is constant, or shared proportionally to the traffic volumes when considering it is load-dependent. Then, when considering only coalitions having the same number of members, we fix the variable related to the number of services, and so the value of a coalition should be only proportional to its traffic volume. Thus, the value of a coalition is the ratio of its traffic volume versus the traffic volume of the coalitions having the same number of members as the given coalition. This explains the characteristic function we defined. Now that the characteristic function of the game is determined, we use the Shapley value concept to compute the payoffs of the players. According to Shapley, the payoff φ i of player i is:

φ i (V, S) = 1 N ! N s=1 (N -s)!(s -1)! C s-1 N -1 j1=1 δ({i}, S) (9)
where δ({i}, S) = V (S j1,{i} )-V (S j1,{i} \{i}) is the marginal contribution of player i, in coalition S. It represents the cost gained or lost by the coalition S with the entry of player i. The computational complexity of ( 9) grows exponentially in the number of service categories, which could represent an obstacle for being implemented. We hence propose a closedform expression for the Shapley value computation, derived from [START_REF] Shapley | Contributions to the theory of games II. Reading, Massachusetts[END_REF].

Let p i denote the traffic proportion of player i: p i = vi vt . The Shapley value of i is then:

φ i (N, p i ) = ( N s=1 1 sC s N )p i + ( N s=2 (C s-2 N -1 -C s-1 N -1 )C s-2 N -2 C s-1 N -1 C s-2 N -1 sC s N )(1 -p i ) (10)
The derivation of this expression is found in Appendix A.

E. Game with a mandatory player

Let us now consider a game with a mandatory player. As stated above, this is the case for instance when an operator has the obligation, by the state, to offer a given service, notably voice, when deploying a given network infrastructure. A mandatory player is such that there can not exist any coalition without him.

Let us denote by i * the mandatory player and o a non mandatory player. The characteristic function of the game is defined as follows:

V (S) = s k1=1 v k1,S C s N j2=1 s k1=1 v k1,Sj 2 1 i * ∈S (11)
φ i * and φ o of the mandatory player i * and a non mandatory player o, respectively, are obtained by [START_REF] Shapley | Contributions to the theory of games II. Reading, Massachusetts[END_REF].

The closed-form expression of the Shapley value of the mandatory player is :

φ i * (N, p i * ) = ( N s=1 1 sC s N )p i * + ( N s=2 C s-2 N -2 C s-1 N -1 sC s N )(1 -p i * ) (12)
The derivation of this expression is found in Appendix B-A. The closed-form expression of the Shapley value of a non mandatory player is :

φ o (N, p i * , p o ) = ( N s=2 (C s-2 N -1 -C s-1 N -1 )C s-2 N -2 C s-1 N -1 C s-2 N -1 sC s N )p i * + ( N s=2 C s-2 N -2 C s-1 N -1 sC s N )p o + ( N s=3 (C s-2 N -1 -C s-1 N -1 )C s-3 N -3 C s-1 N -1 C s-2 N -1 sC s N )(1 -p i * -p o ) (13)
The derivation of this expression is found in Appendix B-B.

F. Game with only mandatory players

Let consider a game where all players are mandatory. The closed-form Shapley value of a player is given by:

φ k (N ) = 1 N (14)
The derivation of this expression is found in Appendix C. Thus, in the special case where all the players are mandatory, the fixed energy consumption is uniformly distributed among the service categories, whatever the traffic volume proportions. The uniform model is a special case of the Shapley-based model.

It is worth to note that the closed-form equations are simplification of the Shapley value equation, reducing significantly the computational complexity of the algorithm, but we lose at the same time the generality of the Shapley value equation. Indeed, since the closed-form equations come from the replacement of the characteristic function V in the Shapley equation (Eqn. ( 9)) by its analytical expressions we defined in our work, these equations are therefore specific to the game we investigate.

G. Heterogeneous radio access network model

When the network is composed of several access technologies, the sharing of the energy consumption is done per radio technology, then we deduce the overall sharing considering all the technologies.

Let T denote the set of radio technologies, F t the set of services served by the sub-network associated with radio technology t, E v t (E f t ) the variable (fixed) energy consumption of the sub-network t, v i,t the traffic volume of service category i on sub-network t. The total network energy consumption (variable and fixed) induced by the service i is :

E i = t∈T ( v i,t k∈Ft v k,t E v t + φ i,t E f t ) (15)
If θ t is the share of the variable component of the energy consumption, then :

E i = t∈T (θ t v i,t k∈Ft v k,t + (1 -θ t ) φ i,t ) E t (16)

H. Case of evolving network infrastructure

Over long periods of time, typically on the order of years, the network infrastructure needs to be upgraded in order to keep up with traffic load increase, as depicted in Fig. 3. When there are several upgrade levels to be modeled, the share of a service category in the fixed energy consumption of the network is derived from its shares per technology and per upgrade level. Let L denote the set of upgrade levels, v i,t,l the traffic volume of the service category i on the sub-network t considering the upgrade level l, E v t,l (E f t,l ) the variable (fixed) energy consumed by the sub-network t on the upgrade level l.

E i = l∈L t∈T ( v i,t,l k∈Ft v k,t,l E v t,l + φ i,t,l E f t,l ) (17)
The traffic variations impact strongly the network upgrade, and so one needs to make again a fair sharing of energy and equipment costs, that take into account this aspect. For this purpose, one can consider the variations of traffic, δv, instead of the traffic volumes, v, in the model.

E i = l∈L t∈T ( δv i,t,l k∈Ft δv k,t,l E v t,l + φ i,t,l E f t,l ) (18)
IV. IMPLEMENTATION ISSUES Fig. 4 shows the runtime (in second) of two algorithms for the computation of the Shapley values of the service categories, one using (9) -denoted by Classical -and the other using the closed-form expression (10) -denoted by Optimized.

The algorithm using the closed-form expression (10) has a runtime almost independent of the number of service categories in the network (less than 1 for up to 50 service categories, the maximum number of service categories we measure in the considered network), while the algorithm using (9) has a computational complexity growing exponentially in the number of service categories, does not converge and has some resource limitation from a certain number of service categories (depending on the hardware and software environment). This comes from that (9) computes the marginal contribution of each service category in all N ! scenarios of constructing the grand coalition, whilst the closed-form expressions we derive from (9) are simple linear functions of the traffic proportions, which can be written in the following general forms:

φ i (N, p i ) = A(N )p i + B(N ) (19)
for [START_REF] Marquet | Scientific methodology for telecom services energy consumption and co2 emission assessment including negative and positive impacts[END_REF],

φ i * (N, p i * ) = C(N )p i * + D(N ) (20)
for [START_REF] Preist | Analyzing end-to-end energy consumption for digital services[END_REF] and

φ o (N, p i * , p o ) = E(N )p i * + F (N )p o + G(N ) (21)
for [START_REF] Vereecken | Power consumption in telecommunication networks: overview and reduction strategies[END_REF], where A(N ), C(N ), E(N ), F (N ) represent the impact of traffic proportions in the sharing of the fixed energy consumption, and B(N ), D(N ), G(N ), the lower bounds of the players' share in the fixed energy consumption. For example, for N = 5, φ i (p i ) = A(5)p i + B(5) = 0.417 p i + 0.117. As depicted in Figs. 5 and6, A(N ), B(N ), C(N ), D(N ), G(N ) are asymptotically equivalent to 1/N . E(N ), F (N ) tend faster to 0. That is (10),( 12) and ( 13) become respectively φ i (N,

p i ) = 1+pi N , φ i * (N, p i * ) = 1+pi * N and φ o (N ) = 1
N for large number of services.

It is worth to note that φ i (p i ) = A(2)p i + B(2) = p i for N = 2, which means the Shapley-based model is equivalent to a proportional distribution when considering just 2 service categories, if none is considered as a mandatory service.

In addition, besides the simplification of the computational complexity, the closed-form expressions give the lower bound of the players' share in the fixed energy consumption. Considering for instance the scenario without a mandatory player, and considering 5 service categories, the minimum share of a service category is B(5) = 12%. Fig. 6 shows that the lower bound of the share of a service category depends on the number of categories defined in the model. Considering again the scenario without a mandatory player, the minimum share of a player is maximal when considering 5 service categories (12%) and minimal with 2 service categories (0%). V. NUMERICAL APPLICATIONS We now turn to the evaluation of our Shapley-based sharing model of energy between different service categories. We consider Orange France's network. The period of the study covers two years representing a mature 2G/3G network with early LTE deployments and associated traffic increase. We measure all voice and data services that are transmitted in the network with the following segmentation for the service categories: two large ones, namely streaming and web browsing, and three smaller ones: download, voice and other minor services. Fig. 7 shows their traffic proportions as taken from the real dataset. We consider just the traffic and energy consumption of the 3G sub-network (the network of NodeBs and RNCs). The variable component of the energy consumption is shared proportionally to the traffic loads, as this component is load-dependent. This implies that data services induced 90% of the UMTS Terrestrial Radio Access Network (UTRAN) variable energy consumption. These services are dominated by Over The Top (OTT) actors like Google.

A. Performance metric

In order to quantify the comparative performance between the three sharing strategies, i.e., Shapley-based, uniform and proportional, we develop next several metrics and show that the third one will be the one we will use in our performance evaluation.

Let S denote the set of strategies of the players. S = {u, p, s} with u denoting the uniform model, p the proportional model, and s the Shapley-based model. Each player has three strategies he can play. Let N denote the set of players, x k i the fixed energy share of player i when playing strategy k. Let xi denote the minimum share of player i, with regard to its strategies.

xi = min((

x k i ) k∈S ) (22)
x = (x i ) i∈N is the vector of minimum shares of the players. Let r k i denote the regret of player i, playing strategy k. It is the difference between its share when playing strategy k and its minimum share.

r k i = x k i -xi (23)
r k = (r k i ) i∈N is the regret vector of the players when strategy k is chosen.

1) Pareto-dominant strategy: A strategy k maximizes the satisfaction of the players if it is Pareto-dominant, that is, the associated regret vector, r k , is such that:

r k i ≤ r k′ i ∀ i ∈ N , ∀ k′ ∈ S (24)
and

∃ i .s.t r k i < r k′ i ( 25
)
We show that no strategy among the three we consider is Pareto-dominant. Let assume that such a strategy exists. Let k and k′ be two given strategies. k is a Pareto dominant strategy implies that

N i=1 x k i < N i=1 x k′ i =⇒ 1 < 1.
This absurd result shows the non existence of a Pareto dominant strategy in this game. We thus turn to the social welfare strategy, which minimizes the mean regret of the players.

2) Social welfare strategy: A strategy k maximizes the satisfaction of the players if it minimizes their mean regret with regard to the others strategies. However, all the strategies offer the same mean regret to the players, as shown next. Let R k denote the cumulative regret of the players, given strategy k.

R k = N i=1 r k i R k = N i=1 x k i -xi R k = N i=1 x k i - N i=1 xi R k = 1 - N i=1 xi R k = 1 - N i=1 xi
shows the cumulative regret of the players is independent of the strategy, then all the strategies have the same mean regret, hence the need to look for another order relation on S.

3) The strategy of the minimum variance: A strategy k maximizes the satisfaction of the players if it minimizes the variance of the associated regret vector, r k , with regard to the other strategies. By minimizing the variance of the regrets, the strategy k minimizes the difference between the regrets of satisfied and those of unsatisfied players. A player is satisfied when its regret is lower than the mean regret of the players, and is unsatisfied otherwise.

var(r k ) = min((var(r k ′ )) k ′ ∈S ) (26)
The idea of the "satisfaction function of the minimum variance" is as follows: Every player has a preferred sharing model, which corresponds to the model that assigns him its lowest responsibility in the fixed energy consumption. When implementing a sharing model, a player has a regret corresponding to the difference between its responsibility given the implemented sharing approach and its lowest responsibility. Given that all the players have the same mean regret whatever the sharing model, their aim is then to minimize the differences between their regrets, hence the minimization of the variance of their regrets. The sharing model that minimizes the variance of their regrets, maximizes their satisfaction.

B. Energy sharing without a mandatory service category

We now turn to the fixed component of the energy consumption and show in Fig. 8 the sharing achieved by our Shapleybased proposal along with two other approaches: uniform sharing between the different service categories, independently of their traffic loads as on the short term the fixed energy consumption is independent of the network traffic load. And a proportional sharing which follows the traffic proportions of the service categories, given that traffic increase over a larger time scale causes network upgrades that in turn augment the fixed energy consumption.

It is worth to notice in the figure that the uniform approach favors "big services" (in terms of load) while "small ones" are favored by the proportional sharing. Our Shapley-based model achieves actually a trade-off among all the players, taking into consideration the double behavior of the fixed energy as it varies or not with the traffic load according to the time scale, unlike the uniform sharing that accounts only for the short term, and the proportional approach for the longer term. Indeed "big players", namely streaming and web services, have a lower impact in the network fixed energy consumption than they would have had with a proportional approach, as well as "small players", namely voice, download and other minor services with regard to a uniform sharing. This is a good trade-off for streaming and web services as it does not penalize them a lot and acknowledges the fact that they are major drivers for network activity, and is also a good trade-off for services with small loads as it does not make them too much responsible of the fixed energy consumption and encourages their transport as well as introduction of yet new, small ones.

The trade-off offers by our Shapley-based model to all the players results in the maximization of their satisfaction for this sharing model, as illustrated in Fig. 9. Maximizing the satisfaction is equivalent to reducing the differences between the highest and the lowest regrets of the players. Based on our Shapley-based model, data services represent 85% of the UTRAN fixed energy consumption, versus 15% for voice service. We deduce that data services represent 0.91θ 3G + 0.85(1θ 3G ) of the total energy consumption (fixed and variable) of the 3G RAN. Typically θ 3G = 0.2 because the access is under-loaded (say ρ = 25%), finally data services represent 86% of the total RAN energy consumption. We consider the power consumption model of the base station in Fig. 1, i.e., P (ρ) = 0.62(1 + ρ). For ρ = 25%, θ 3G = 0.775-0.62 0.775 = 0.2.

C. Energy sharing with a mandatory service category: voice

We now turn to the case where the voice service is mandatory due to legal constraints. In this scenario, voice is not considered in the selection of the best sharing approach since it is a mandatory player, then must play whatever the sharing model.

Based on Fig. 10, the best sharing model can not be the proportional approach as the regret of big services is very high. Uniform sharing is also eliminated because it induced higher regrets for all the players with regard to the Shapley-based model. The Shapley-based sharing appears as the strategy that minimizes the difference between the players regrets, and thus maximizing their satisfaction, as depicted in Fig. 11.

It is worth to notice that the Shapley-based model takes into account the mandatory nature of the voice service by augmenting significantly its share in the fixed energy consumption (from 15% to 29%). This results in a significant reduction of the impact of data services on the total energy consumption of the network. Data services represent now 57% of the total energy consumption, that corresponds to a decrease of 29% compared to the scenario where voice is not a mandatory player. It is a good outcome that the responsibility of the mandatory player, voice service in this work, in the network's fixed energy consumption significantly increases because the operator would be mandated by the government to implement and offer it on a national basis. And so, to offer this service, the operator would need to deploy a network and dimension it in such a way so as to reach all the citizens of the given country. In this case, the network can be seen as initially deployed to transport primarily this service, and so it is natural that it would share a large part of the responsibility of the energy consumption (yet not all of it since it is sharing the infrastructure with subsequent services).

D. Case of evolving network infrastructure

We now study the case where the network infrastructure is upgraded due to a traffic increase. Fig. 12 shows real measurements for the traffic volumes of the same service categories over two periods of time corresponding to two network upgrade levels, termed levels 1 and 2 in the figure. We consider only the 3G traffic and 2 upgrade levels. The sharing of the variable component of the energy consumption follows here too the traffic proportions, and so the corresponding figure is omitted. For the fixed component however, Fig. 13 shows the new sharing based on [START_REF] Marquet | Scientific methodology for telecom services energy consumption and co2 emission assessment including negative and positive impacts[END_REF], at each level, considering both the traffic volumes and the variations of traffic, at level 2.

Considering the traffic volumes puts the weight on big services, while considering the traffic variations puts the weight on service categories whose traffic increases rapidly, taking into account the impact of traffic increase on the necessity of upgrading the network infrastructure (adding of new equipment in the network), that is itself responsible of an increase in the energy consumption, both fixed and variable.

Note that the traffic increase is correlated with the traffic volumes in the studied network. This results in a similar sharing (in terms of weights) of the fixed energy consumption among the services when considering either traffic volumes or traffic variations.

VI. CONCLUSION

We investigated in this work the sharing of the energy consumed by a wireless access network among the provided service categories. We focused on the two energy components: small, load-dependent variable one and significantly larger fixed one (since the RAN is most of the time under-loaded with traffic load ρ < 50%), load-independent over the short term but load-dependent over longer period of time (typically years). The former is to be shared among the different service categories according to their load proportions. As of the fixed component, it can be shared in a variety of ways: the simplest one being an equal share between the players -this is unfair towards players with small traffic loads. Another one is a sharing proportional to the traffic loads -this puts too much weight on the players with large volumes of traffic, which are a major driving force for the network to be operating. We proposed in this paper a third approach, based on the Shapley value, which put less weight on small players than equal sharing, and less weight on big players than proportional sharing. This is appreciable as it encourages transport and introduction of small services, and acknowledges the role of larger services as major drivers for network activity.

We considered two settings: one with constant network infrastructure and one with evolving network infrastructure over longer periods of time. Moreover, we considered the case of mandatory services wherein some service categories are legally obliged to be provided, such as voice service in several deployed operator networks.

Our next work will focus on the end-to-end path, from the content location in a datacenter for instance to the end user, and on the quantification as well as the sharing of the total energy consumption, again, among the different service categories in the overall network.

APPENDIX A CLOSED-FORM SHAPLEY VALUE OF A PLAYER IN A GAME WITHOUT MANDATORY PLAYERS

The characteristic function of the game without a mandatory player is as follows:

V (S) = s k1=1 v k1,S C s N j2=1 s k1=1 v k1,Sj 2
The value of a coalition S is the ratio of the traffic volume of that coalition and the traffic volume of all the coalitions having the same size as S, whose number is C s N , and size is s. v k,S is the traffic volume of the k th element of the coalition S. The marginal contribution of player i is the gain or loss of the coalition S to the entry of player i in the coalition. It is determined as follows:

V (S) -V (S\{i}) = s k1=1 v k1,S C s N j2=1 s k1=1 v k1,Sj 2 - s-1 k2=1 v k2,S\{i} C s-1 N j4=1 s-1 k2=1 v k2,Sj 4 Then, C s-1 N -1 j1=1 V (S j1 ) -V (S j1 \{i}) = C s-1 N -1 j1=1 s k1=1 v k1,S j 1 ,{i} C s N j2=1 s k1=1 v k1,Sj 2 - C s-1 N -1 j1=1 s-1 k2=1 v k2,Sj 1 \{i} C s-1 N j4=1 s-1 k2=1 v k2,Sj 4 C s-1
N -1 is the number of coalitions of size s containing player i, coalitions. In fact there are C s-2 N -2 coalitions of size s with both player i and a given player k. Then we have:

C s-1 N -1 j1=1 V (S j1 ) -V (S j1 \{i}) = C s-1 N -1 v i + C s-2 N -2 N k3=1 k3 =i v k3 C s-1 N -1 N k4=1 v k4 - C s-2 N -2 N k3=1 k3 =i v k3 C s-2 N -1 N k4=1 v k4 φ i (v) = 1 N ! N s=1 (N -s)!(s -1)! C s-1 N -1 v i + C s-2 N -2 N k3=1 k3 =i v k3 C s-1 N -1 N k4=1 v k4 - 1 N ! N s=1 (N -s)!(s -1)! C s-2 N -2 N k3=1 k3 =i v k3 C s-2 N -1 N k4=1 v k4 C s N = N ! s!(N -s)! =⇒ (N -s)!(s -1)! = N ! sC s N Hence, φ i (v) = 1 v T ( N s=1 1 sC s N )v i + 1 v T ( N s=2 (C s-2 N -1 -C s-1 N -1 )C s-2 N -2 C s-1 N -1 C s-2 N -1 sC s N )(v t -v i )
Let p i denote the traffic proportion of player i.

p i = v i v T φ i (N, p i ) = ( N s=1 1 sC s N )p i + ( N s=2 (C s-2 N -1 -C s-1 N -1 )C s-2 N -2 C s-1 N -1 C s-2 N -1 sC s N )(1 -p i )

APPENDIX B GAME WITH A MANDATORY PLAYER

A. Closed-form Shapley value of the mandatory player

Let i * denote the mandatory player. The value of a coalition with the mandatory player is :

V (S) =
v k5 C s-1 N -1 N k4=1 v k4 - C s-2 N -2 v i * +C s-3 N -3 N k5=1 k5 =i * ,j v k5 C s-2 N -1 N k4=1 v k4
In the case of a mandatory player, any non mandatory player can not form a coalition of less than 2 members.

φ o (v) = 1 v T ( N s=2 (C s-2 N -1 -C s-1 N -1 )C s-2 N -2 C s-1 N -1 C s-2 N -1 sC s N )v i * + 1 v T ( N s=2 C s-2 N -2 C s-1 N -1 sC s N )v o + 1 v T ( N s=3 (C s-2 N -1 -C s-1 N -1 )C s-3 N -3 C s-1 N -1 C s-2 N -1 sC s N )(v T -v i * -v o ) φ o (N, p i * , p o ) = ( N s=2 (C s-2 N -1 -C s-1 N -1 )C s-2 N -2 C s-1 N -1 C s-2 N -1 sC s N )p i * + ( N s=2 C s-2 N -2 C s-1 N -1 sC s N )p o + ( N s=3 (C s-2 N -1 -C s-1 N -1 )C s-3 N -3 C s-1 N -1 C s-2 N -1 sC s N )(1 -p i * -p o )

APPENDIX C CLOSED-FORM SHAPLEY VALUE OF A PLAYER IN A GAME WITH ONLY MANDATORY PLAYERS

Let us consider a game where all players are considered as mandatory. In this game, only the grand coalition can be formed, whose value is V (N ), the fixed energy consumption to share between services. According to [START_REF] Shapley | Contributions to the theory of games II. Reading, Massachusetts[END_REF], the payoff of a player k is:

φ k = 1 N ! (N -N )!(N -1)!
s is always equal to N and the marginal contribution of any player in the grand coalition is the value of the grand coalition as all players are mandatory. V (N ) -V (N \{k}) = 1, since V (N \{k}) = 0, ∀k.

As a reminder, payoffs and values are normalized by the fixed energy consumption, i.e., V (N ).

Hence,

φ k (N ) = 1 N
We showed that the uniform sharing is a special case of the Shapley-based sharing, when all the players of the game are mandatory. This ends the proof.

Fig. 1 .

 1 Fig. 1. Power consumption of a 4G base station.

Fig. 2 .

 2 Fig. 2. The system.

Fig. 3 .

 3 Fig. 3. Evolution of the traffic load and network infrastructure.

Fig. 4 .

 4 Fig. 4. Runtimes of the classical-based and closed-form-based Shapley value algorithms.

Fig. 5 .Fig. 6 .

 56 Fig. 5. Impact of traffic proportions in the fixed energy sharing.

Fig. 7 .

 7 Fig. 7. Traffic proportions per service category.

Fig. 8 .

 8 Fig.8. Sharing of the fixed energy using three approaches: uniform, proportional and Shapley.

Fig. 9 .

 9 Fig. 9. Players' satisfactions per sharing approaches.

Fig. 11 .

 11 Fig. 11. Players' satisfactions per sharing approaches -voice a mandatory player.

Fig. 13 .

 13 Fig. 13. Sharing of the fixed energy -case of evolving network infrastructure.

  v k2,Sj 4 are respectively the traffic volumes of the coalitions of size s and s -1. They are constant for a given coalition size s, hence we can get them out of the sum over coalitions of same size.

2 N - 2

 22 k1,S j 1 ,{i} is the sum of the traffic volumes of the coalitions of size s containing player i. v k,S j,{i} is the traffic volume of the k th element of the j th coalition of size s containing player i. Player i of course is present in all the C s-1 N -1 coalitions, while each other player is present in C s-

1 N - 1 coalitions

 11 v k2,Sj 1 \{i} is the sum of the traffic volumes of the coalitions of size s-1 not containing player i. v k,Sj \{i} is the traffic volume of the k th element of the j th coalition of size |S j | not containing player i. |.| is the cardinal function. Similarly, a given player k appears in C s-2 N -2 coalitions among the C s-1 N -1 coalitions of size s -1, derived from the C s-of size s containing i.

C s- 1 N - 1 j1=1V 2 C s- 1 N - 1 j1=1V 2 N - 2 NC s- 2 N - 2 C s- 1 N - 1 )( 1 -V 4 C s- 2 N - 2 j3=1V 4 C s- 2 N - 2 j3=1V

 112112222111422422 s)!(s -1)! (S j1,{i * } )In fact V (S\{i * }) = 0 ∀ S (S j1,{i} ) = C s-1 N -1 v i + C s-sC s N Tv i * ) φ i * (N, p i * ) = ( p i * )B. Closed-form Shapley value of a non mandatory playerLet o denote a non-mandatory player. The value of a coalition with a non mandatory player and the mandatory player is: (S j3,{i * ,o} \{o})C s-2N -2 is the number of coalitions of size s containing both players i * and o.V (S) -V (S\{o}) = (S j3,{i * ,o} ) -V (S j3,{i * ,o} \{o}) = (S j3,{i * ,o} ) -V (S j3,{i * ,o} \{o}) = C s-2 N -2 v i * +C s-2 N -2 v o + C s-

3 N 3 N

 33 s)!(s -1)! C s-2 N -2 v i * +C s-2 N -2 v o + C s-s)!(s -1)! C s-2 N -2 v i * +C s-

ACKNOWLEDGMENT Under the framework of Celtic-Plus, the studies of the SooGREEN project have been partially funded by the French Directorate General for Enterprise (DGE) and the Scientific and Technological Research Council of Turkey (TUBITAK). c 2016. This manuscript version is made available under the CC-BY-NC-ND 4.0 license http://creativecommons.org/ licenses/by-nc-nd/4.0/