M Mandjes
email: m.r.h.mandjes@uva.nl

P G Taylor
email: p.taylor@ms.unimelb.edu.au

Koen De

K De Turck
email: koen.deturck@l2s.centralesupelec.fr

The Markov-modulated Erlang loss system

Keywords: . Markov-modulation, loss systems, Erlang loss formula, scaling limits

published or not. The documents may come

Introduction

The Erlang loss model, or M/M/C/C queue, dates back to Erlang's original paper [START_REF] Erlang | Solution of some problems in the theory of probabilities of significance in automatic telephone exchanges[END_REF] in 1917. It models a setting in which calls arrive at a telephone exchange according to a Poisson process with rate λ, call holding times are exponentially distributed with mean µ, and C lines are available. The blocking probability turns out to have the clean, explicit form embodied in the expression

p bl = 1 C! λ µ C C =0 1 ! λ µ . (1)
Equation ( 1) is often referred to as the Erlang-B formula. Since Erlang's time, the model has been broadly applied for dimensioning circuit-switched telephone networks. Arguably, [START_REF] Anderson | A functional central limit theorem for a Markov-modulated infinite-server queue[END_REF] is the most frequently used formula in communications engineering.

The model has been extended in many directions. In the first place, Erlang himself analyzed the model in which customers who find all servers busy are not lost but rather join a queue, leading to the Erlang-C formula. Later, it was found that p bl depends on the service time distribution only through its mean; for nice accounts of this remarkable insensitivity property, see [START_REF] Miyazawa | Insensitivity and product form decomposability of reallocatable GSMP[END_REF][START_REF] Whittle | Systems in Stochastic Equilibrium[END_REF].

A further generalization is the loss network [START_REF] Kelly | Loss networks[END_REF], which was defined about 25 years ago to study the probabilistic properties of a network with R user classes (characterized by their Poisson arrival processes and exponential service times) and K links, with C k circuits being available at link k. On arrival, a user of class r seizes a specified number A rk circuits from each link k if this number of circuits is available, or is rejected and 'lost' from the system if there is a link at which the required A r circuits are not free. Loss networks also have the insensitivity property with respect to the distribution of the holding times of each of the user classes [START_REF] Kelly | Loss networks[END_REF].

Even with the advent of modern packet-switched data networks, the loss model (together with its many variants) has many uses in communications engineering via concepts such as effective bandwidth [START_REF] Elwalid | Effective bandwidth of general Markovian traffic sources and admission control of high speed networks[END_REF], and also in a broad range of other resource allocation problems.

One obvious variant of the Erlang loss model is the Markov-modulated M/M/C/C queue. The point of difference between this model and the standard Erlang loss model is that the arrival and service rates vary over time. More specifically, the system has Poisson arrivals and exponential service time distributions, but with the special feature that the arrival rates and service rates are determined by an external, independently-evolving, continuous-time Markov chain. This queue considered is still a loss system, with customers arriving to find all servers busy being rejected from the system and lost. However, despite the fact that M. Neuts analyzed the related Ph/M/C/K + C model in his celebrated monograph [16, pp. 92-94], the Markov-modulated loss queue does not appear to have received much attention in the literature.

The evolution of the Markov-modulated M/M/C/C queue is uniquely described by the bivariate continuous-time Markov chain (M (t), J(t)) t 0 , where M (t) denotes the number of customers in the system, and J(t) the state of the modulating Markov chain (also referred to as the background process), at time t. In the literature on matrix-analytic methods, the process (M (t)) t 0 is usually referred to as the level process, whereas (J(t)) t 0 is called the phase process. The process (J(t)) t 0 is assumed to be irreducible, taking values in some finite state space S := {1, . . . , d}. We use Q = (q ij ) i,j∈S to denote the generator of the modulating Markov chain. As usual, q ij 0 for i = j and Qe = 0, with e denoting a (d × 1) column vector of ones. The (1 × d) vector α denotes the stationary distribution of the Markov chain (J(t)) t 0 : it is the unique non-negative vector satisfying αQ = 0 and αe = 1.

The vectors λ = (λ j ) j∈S and µ = (µ j ) j∈S contain the arrival rates and per-customer service rates. At time t the arrival rate is λ J(t) , whereas the total service rate at time t is µ J(t) M (t), proportional to the number of customers currently present. The number of arrivals between 0 and t is Poisson with the random parameter t 0 λ J(s) ds.

(

) 2
Under the assumption that at least one λ j and one µ j are positive, the irreducibility of Q leads to the fact that the whole queue is modelled by an irreducible finite state Markov chain. As a consequence, the stationary distribution

π k,j := lim t→∞ P(M (t) = k, J(t) = j) (3)
exists for k ∈ {0, . . . , C} and j ∈ S. We shall use (M, J) to denote random variables distributed according to this stationary distribution.

In this paper we provide an explicit analysis of the stationary distribution (3) of the Markovmodulated M/M/C/C queue, which we encode by its parameters (λ, µ, Q, C). In addition, we provide remarkably simple, yet accurate approximations in particular scaling regimes. More specifically, the main contributions of our work are the following.

• We derive an explicit expression for the probabilities π k ≡ (π k,j ) j=1,...,d . In particular, we identify the blocking probability in closed form, which can be considered the true Markovmodulated counterpart of the classical formula (1).

• We consider the scaling (N λ, µ, N f Q, C N ) for N large; here C N := N + βN γ , with the system's offered load (αλ)/(αµ), β > 0 a given constant (typically referred to as the 'hedge'), and γ := max{1 -f /2, 1/2}. In this central-limit type of scaling, we show that an appropriately centered and scaled version of the stationary number of customers converges to a truncated Normal random variable. This finding also provides us with a compact, closed-form approximation of the blocking probability in this 'fast' regime. Importantly, the parameters of the truncated Normal random variable critically rely on whether f < 1, f = 1, or f > 1.

The proof relies on a variant of the methodology developed in [START_REF] Blom | Analysis of Markov-modulated infinite-server queues in the central-limit regime[END_REF].

• In addition, we study a 'slow' regime of the type (λ, µ, εQ, C), with ε ↓ 0. Relying on the theory of singularly perturbed Markov chains, we propose explicit approximations in this case.

• The paper concludes with a series of representative experiments. In both the fast and slow regime, the proposed approximations perform remarkably well, despite their simple form.

The organization of our paper is as follows. Section 2 presents the explicit results for the non-scaled model (λ, µ, Q, C). Results for the fast and slow regimes are covered in Section 3. The numerical experiments can be found in Section 4.

Exact analysis

The objective of this section is to identify the blocking probability in the Markov-modulated counterpart of the M/M/C/C queue characterized by the parameters (λ, µ, Q, C). We do so by providing an expression for the vector of probabilities π C that M = C while J runs over the various values in S. In addition, we identify the full distribution π k ≡ (π k,j ) for k ∈ {0, . . . , C} and j ∈ S.

In the sequel we frequently work with probability generating functions. For k ∈ {0, . . . , C}, we write π k for the 1 × d vector whose entries are π k,j , with j ∈ S. Then

Π(z) := C k=0 π k z k
is the vector-valued probability generating function of the level process, defined for all z ∈ C.

We first derive the equations that determine the stationary distribution of (M (t), J(t)). Denoting by 1 A the indicator of the event A, it is immediate that the vectors π k , with k = 0, . . . , C, satisfy the equations

π k Q -1 {k<C} diag{λ} -kdiag{µ} + π k-1 1 {k>0} diag{λ} + π k+1 1 {k<C} (k + 1)diag{µ} = 0, (4)
with diag{η} the diagonal matrix whose entries are given by the vector η. Multiplying Equation (4) by z k and summing over k = 0, . . . , C, we derive the fact that Π(z) satisfies

Π(z)Q + (z -1) Π(z) -π C z C diag{λ} - dΠ(z) dz diag{µ} = 0. (5)
By substituting z = 1 in Equation ( 5), we obtain the equation

Π(1)Q = 0, (6)
which makes sense, because Π(1) = C k=0 π k is the marginal stationary distribution α of the phase, which must satisfy αQ = 0. Combining this with the obvious identity Π(1)e = 1, Equation (6) uniquely defines Π(1) = α.

Differentiating Equation [START_REF] Blom | Markov-modulated infinite-server queues with general service times[END_REF] with respect to z and putting z = 1 we see that

Π (1) [Q -diag{µ}] + (Π(1) -π C ) diag{λ} = 0. (7)
Because Q -diag{µ} is an irreducible generator with at least one row sum strictly less than zero, it must be non-singular (see, for example, [17, pp. 8 and 31]). It follows that we can rewrite Equation [START_REF] Auria | Stochastic decomposition of the M/G/∞ queue in a random environment[END_REF] as

Π (1) = [Π(1) -π C ] diag{λ} [diag{µ} -Q] -1 . (8)
We can also derive a similar relationship for the higher derivatives of Π(z), again evaluated at z = 1. Repeated differentiation of [START_REF] Blom | Markov-modulated infinite-server queues with general service times[END_REF] shows us that, for k = 1, . . . , C -1,

Π (k) (z)Q + k Π (k-1) (z) - C! (C -k + 1)! z C-k+1 π C diag{λ} -Π (k) (z)diag{µ} + (z -1) Π (k) (z) - C! (C -k)! z C-k π C diag{λ} -Π (k+1) (z)diag{µ} = 0,
where we have written Π (k) (z) for the k-th derivative of Π(z). Inserting z = 1 we see that

Π (k) (1) [Q -k • diag{µ}] + k Π (k-1) (1) - C! (C -k + 1)! π C diag{λ} = 0. (9)
The matrix Q -k • diag{µ} is nonsingular for the same reason that Q -diag{µ} is nonsingular, and so we can write

Π (k) (1) = k Π (k-1) (1) - C! (C -k + 1)! π C diag{λ} [k • diag{µ} -Q] -1 . (10)
We are now in a position to determine the vector π C . It follows from Equations ( 8) and [START_REF] Huang | Large deviations for Markov-modulated diffusion processes with rapid switching[END_REF] that

Π (k) (1) = k! Π(1) C i=1 diag{λ} [i • diag{µ} -Q] -1 -k! π C k =1 C -1 k i= diag{λ} [i • diag{µ} -Q] -1 (11)
where the matrix products multiply on the right as i runs through successive values.

As an immediate consequence of the fact that the level process attains values in {0, . . . , C}, observe that, for any value of z,

Π (C) (z) = d C dz C C k=0 π k z k = C! π C .
Substituting this into Equation [START_REF] Jagerman | Some properties of the Erlang loss function[END_REF], we see that

π C C =0 C C i= +1 diag{λ} [i • diag{µ} -Q] -1 = α C i=1 diag{λ} [i • diag{µ} -Q] -1 , (12)
with the empty product taken to be the identity matrix. We thus arrive at the following theorem.

Theorem 2.1. In the Markov-modulated M/M/ C / C queue with parameters (λ, µ, Q, C), the vector π C satisfies equation [START_REF] Kelly | Loss networks[END_REF]. If the matrix

A(λ, µ, Q, C) := C =0 C C i= +1 diag{λ} [i • diag{µ} -Q] -1 (13)
is nonsingular, then π C is given by the expression

π C = α C i=1 diag{λ} [i • diag{µ} -Q] -1 C =0 C C i= +1 diag{λ} [i • diag{µ} -Q] -1 -1 . ( 14
)
Remark 2.1. The expression ( 14) is the Markov-modulated analogue of the celebrated Erlang loss formula. The fraction of time that all servers are busy and the background process is in phase j is given by (π C ) j and the fraction of time that all servers are busy irrespective of the state of the background process is π C e. In all instances we considered, the matrix A(λ, µ, Q, C) turned out to be nonsingular, and we conjecture that this is true for all (λ, µ, Q, C). However, a proof of this conjecture has eluded us to date.

Expression ( 14) also provides us with an explicit formula for the customer blocking probability p bl , defined as the probability that an arbitrary arriving customer finds all servers busy. It is obtained by weighting the probabilities π C,j by the fractions λ j /α λ.

Corollary 2.1. The blocking probability p bl is given by

p bl = π C λ α λ = d j=1 π C,j λ j d j=1 α j λ j .
In fact, the above analysis not only yields the vector π C , but actually all π k for k ∈ {0, . . . , C}. To this end, realize that

πk := Π (k) (1) = d k dz k C =0 π z z↑1 = C =k ! ( -k)! π ;
having identified the vectors πk , it is trivial to obtain π k . Directly from Equation [START_REF] Jagerman | Some properties of the Erlang loss function[END_REF] and Theorem 2.1,

Π (k) (1) = k! π C C =0 C C i= +1 diag{λ} [i • diag{µ} -Q] -1 - k! π C k-1 =0 C k i= +1 diag{λ} [i • diag{µ} -Q] -1 .
In the special case where λ i ≡ λ and µ i ≡ µ, for all i ∈ S, we can verify that

πk = 1 B(C) C =k 1 ( -k)! λ µ , B(C) := C =0 1 ! λ µ .
In particular, as desired, πC = C! π C = (λ/µ) C /B(C), demonstrating that the classical Erlang loss formula (1) is indeed a special case of our model.

We can find the moments of the stationary number of clients in the system M , jointly with the state J of the background process. It takes an elementary computation to verify that

EM 1 {J=1} , . . . , EM 1 {J=d} = (α -π C ) diag{λ} [diag{µ} -Q] -1 = C-1 k=0 π k diag{λ} [diag{µ} -Q] -1 ,
and likewise

EM (M -1)1 {J=1} , . . . , EM (M -1)1 {J=d} = 2 Π (1) (1) -Cπ C diag{λ} [2 • diag{µ} -Q] -1 = 2 C-1 k=0 k π k diag{λ} [2 • diag{µ} -Q] -1 .
These expressions facilitate the computation of Var M .

Remark 2.2. The naïve way to numerically compute the process's equilibrium distribution amounts to solving the (C + 1)d-dimensional system of balance equations, which (in its basic form) takes O(C 3 d 3 ) operations. This can already be substantially be reduced by exploiting the system's band structure; from the results in e.g. [START_REF] Thorson | Gaussian elimination on a banded matrix[END_REF] it immediately follows that in our model the computational cost drops to

(C + 1)d(d + 1) 2 /4 = O(Cd 3
). Careful counting reveals that evaluation of ( 14) has complexity O(Cd 3 ) as well. This means that such methods may be prohibitively slow when C or d is large (where the value of d has more impact than the value of C).

Asymptotic analysis

As argued in Remark 2.2, the evaluation of the quantities derived in the previous section may be time consuming in specific regimes, in particular when C and/or d is large. This motivates interest in considering scalings under which more explicit results can be derived. In this section we consider two scalings, which we refer to as the fast regime and the slow regime. The fast scaling, which directly relates to the scalings proposed in [START_REF] Anderson | A functional central limit theorem for a Markov-modulated infinite-server queue[END_REF][START_REF] Blom | Analysis of Markov-modulated infinite-server queues in the central-limit regime[END_REF][START_REF] Blom | Markov-modulated infinite-server queues with general service times[END_REF], involves blowing up the arrival rates and the transition rates of the background process, but, importantly, at different rates. In the slow regime, the transitions of the background process are made increasingly rare. For convenience, we assume henceforth that λ j > 0 and µ j > 0 for all j ∈ S.

The fast regime

We proceed by introducing the scaling of the fast regime in more detail; it involves adaptation of (i) the arrival rates λ, (ii) the transition rate matrix Q, and (iii) the capacity C. We first define by the load imposed on the system: with

λ ∞ := α λ = i∈S α i λ i and µ ∞ := α µ = i∈S α i µ i , this offered load is defined as := λ ∞ µ ∞ .
The scaling regime is now defined as follows:

• We linearly scale the arrival rates by a factor N : λ → N λ (and leave the service rates unaltered).

• The jumps of the modulating process are sped up by a factor N f , for some f > 0, which amounts to Q → N f Q.

• In, for example, [START_REF] Blom | Analysis of Markov-modulated infinite-server queues in the central-limit regime[END_REF] it was pointed out that in the corresponding infinite-server queue the variance of the stationary number of customers essentially grows like N 2γ , with the parameter

γ defined as 1 2 max{2 -f, 1} ∈ [ 1 2 , 1]
. This motivates us to consider a scaling of the capacity

C → C N , with C N := N + βN γ ,
for β > 0. The parameter β is usually referred to as the hedge: the higher the hedge, the lower the blocking probability.

As both the arrival process and the jump process of the modulating Markov chain (J(t)) t 0 are accelerated (albeit at different rates), this regime can be thought of as being fast.

Denote the stationary number of clients in the resulting scaled system by the random variable M (N ) ∈ {0, . . . , C N }. In self-evident notation, the objective of this section is to study the stationary number of clients in the scaled version (N λ, µ, N f Q, C N ) of our original system (λ, µ, Q, C), when N grows large.

We now further motivate the relevance of this scaling regime, by pointing at a link with the phenomenon of overdispersion. As was mentioned in the introduction, the number of arrivals between times 0 and t has a mixed Poisson distribution, where the random parameter is given by ( 2). In the asymptotic regime introduced above, more explicit statements can be made: as shown in [START_REF] Lu | A functional central limit theorem for Markov additive arrival processes and its applications to queueing systems[END_REF], the number of arrivals up to time t in the scaled system, denoted by A (N ) (t), obeys a (functional) central limit theorem with scale function N γ . More specifically, as N grows large,

A (N ) (t) -N λ ∞ t N γ
converges to a zero-mean Normally distributed random variable. It entails the fact that the mean of the number of arrivals grows essentially linearly in N , whereas the variance behaves as N 2γ . We observe that the case where f > 1 (and therefore γ = 1 2 ) corresponds to the 'traditional Poisson regime', in which the mean and the variance of the number of arrivals roughly match. The situation where f < 1 (in which γ > 1 2 ), however, exhibits overdispersion [START_REF] Bassamboo | Capacity sizing under parameter uncertainty: safety staffing principles revisited[END_REF]: the variance grows faster than the mean. The phenomenon of overdispersion has been observed in various operational contexts such as call center data [START_REF] Brown | Statistical analysis of a telephone call center: A queueing-science perspective[END_REF] and hospital arrivals [START_REF] Kim | Are call center and hospital arrivals well modeled by non-homogeneous Poisson processes[END_REF].

The boundary case when f = 1 is more subtle and has to be dealt with separately: the variance behaves linearly, as is the case for f > 1, but the system is more variable than in the Poisson regime. This is expressed by the fact that the mean and variance are both essentially proportional to N , but the proportionality constant corresponding to the variance is larger than that corresponding to the mean.

Before we can state our main scaling result, we have to introduce new notation. Define by D the deviation matrix corresponding to the background process J(•); its (i, j)-th entry is given by

D ij = ∞ 0 (e Qt ) ij -α j dt,
for i, j = 1, . . . , d. This allows us to introduce σ 2 := τ 2 1 {f 1} + 1 {f 1} , with

τ 2 := U µ ∞ , where U := α (diag{λ} -diag{µ}) D (diag{λ} -diag{µ}) e; (15)
below, in Corollary 3.1, we show that U plays a central role in an asymptotically exact approximation of p (N ) bl . Our main result in this section shows that an appropriately centered and normalized version of M (N ) converges in distribution to a truncated Normal random variable; the variance of this Normal distribution critically depends on whether f is smaller or larger than 1. Here and in the sequel, let N(µ, σ 2 ) denote a random variable that has a Normal distribution with mean µ and variance σ 2 . Proposition 3.1. Consider the scaling (N λ, µ, N f Q, C N ). The random variable

M (N ) := M (N ) -N N γ converges, as N → ∞, to (N(0, σ 2 ) | N(0, σ 2 ) β).
With π

(N ) k
denoting the counterpart of π k in the N -scaled model, define

P (N ) (s) := E e -s M (N ) 1 {J=1} , . . . , E e -s M (N ) 1 {J=d} = C N k=0 exp -s k -N N γ π (N ) k .
In order to prove Proposition 3.1, by Lévy's convergence theorem [START_REF] Williams | Probability with Martingales[END_REF]Thm. 18.1] it suffices to establish the convergence

P (N ) (s)e → E(exp(-sN(0, σ 2 )) | N(0, σ 2 ) β), (16)
as N → ∞. From ( 16) it then follows after a routine calculation that we should prove P (N ) (s)e → e σ 2 s 2 /2 P(N(0, 1) σs + β/σ) P(N(0, 1) β/σ) .

Several techniques can be used to prove this convergence. Perhaps the most straightforward among these is a variation of that presented in [START_REF] Blom | Analysis of Markov-modulated infinite-server queues in the central-limit regime[END_REF], of which we have included a sketch in Appendix A.

The following corollary of the derivation in Appendix A provides us with the asymptotics of the blocking probability, with this probability in the N -scaled model denoted by

p (N ) bl := π (N ) C λ α λ = d j=1 π (N ) C N ,j λ j d j=1 α j λ j .
Let φ N (•) and Φ N (•) denote the density and cumulative distribution function of a standard Normal random variable, respectively.

Corollary 3.1. Consider the scaling (N λ, µ, N f Q, C N ). As N → ∞, N 1-γ p (N ) bl → µ ∞ λ ∞ 0 -∞ e -r 2 σ 2 /2-rβ dr -1 = σ φ N (β/σ) Φ N (β/σ) =: b. ( 18
)
The above corollary is the Markov-modulated counterpart of the Halfin-Whitt-type result for the non-modulated loss system, as has been presented in e.g. [START_REF] Whitt | Stochastic-process Limits[END_REF]Thm. 10.4.3]. The non-modulated result [START_REF] Whitt | Stochastic-process Limits[END_REF]Thm. 10.4.3] states that, under essentially the same scaling as the one that we consider, the blocking probability is inversely proportional to √ N (where the proportionality constant was explicitly given). Corollary 3.1 entails that in the Markov-modulated case the blocking probability is inversely proportional to N 1-γ , rather than √ N , where it is noted that

N 1-γ = √ N for f 1.
Now we study the asymptotics of the blocking probability more closely, To this end, consider the f > 1, f < 1, and f = 1 cases separately. We denote by

g N ∼ h N that g N /h N → 1 as N → ∞. • For f > 1, Corollary 3.1 yields (realizing that γ = 1 2 ), p (N ) bl ∼ b √ N = 1 √ N 1 √ φ N (β/ √ ) Φ N (β/ √ ) ,
in line with the (classical) asymptotics for the non-modulated M/M/C/C queue under a Halfin-Whitt-type scaling [START_REF] Jagerman | Some properties of the Erlang loss function[END_REF] and [START_REF] Whitt | Stochastic-process Limits[END_REF]Thm. 10.4.3]. For the non-scaled model (λ, µ, Q, C), this leads to the approximation

p bl ≈ 1 √ φ N ((C -)/ √ ) Φ N ((C -)/ √ ) . (19)
• For f < 1, we have γ = 1 -f /2, and therefore Corollary 3.1 yields

p (N ) bl ∼ b √ N f = 1 √ N f √ U µ ∞ λ ∞ φ N (β µ ∞ /U ) Φ N (β µ ∞ /U ) ,
with U as in [START_REF] Miyazawa | Insensitivity and product form decomposability of reallocatable GSMP[END_REF]. For the non-scaled model (λ, µ, Q, C), we thus obtain the approximation

p bl ≈ √ U µ ∞ λ ∞ φ N ((C -) µ ∞ /U ) Φ N ((C -) µ ∞ /U ) . (20)
• For f = 1, again due to Corollary 3.1, we see that

p (N ) bl ∼ b √ N = 1 √ N 1 √ U λ ∞ -1 φ N (β/ τ 2 + ) Φ N (β/ τ 2 + ) .
We observe the √ N scaling as in the case f > 1, but with a different proportionality constant. This yields for the non-scaled model

p bl ≈ U/µ ∞ + • φ N C - U/µ ∞ + Φ N C - U/µ ∞ + . ( 21
)
In Section 4 we get back to these approximations, also proposing a single approximation that works for all three regimes simultaneously, across all values of f > 0.

The slow regime

A crucial characteristic of the regime discussed in Section 3.1 is that the modulating Markov chain jumps relatively fast. In this section we consider the situation in which the opposite applies: we study the scaling (λ, µ, εQ, C) for ε ↓ 0. Note that for ε = 0, the system decomposes into d different ergodic classes, each constituting an (unmodulated) Erlang loss model, while for any positive ε, the system is ergodic. As a consequence, the regime ε ↓ 0 can be analyzed relying on the theory of singularly perturbed Markov chains.

To this end, let Q(ε) denote the family of generator matrices of the process (M (t), J(t)) t 0 , indexed by the 'rarity parameter' ε. For notational convenience, we swap levels and phases compared with the rest of the paper: the levels represent the background state and the phases represent the number of customers present. Thus, Q(ε) can be decomposed as follows: with ⊗ as usual denoting the Kronecker product,

Q(ε) =     Q1 . . . Qd     + ε • (Q ⊗ I) =: Q (0) + ε • Q (1)
,

where Qi represents the (C + 1) × (C + 1)-dimensional generator matrix of a (non-modulated) Erlang loss model with parameters (λ i , µ i , C). We let ν i denote the corresponding the stationary distribution vector: for j = 0, . . . , C,

ν i,j = 1 j! λ i µ i j C =0 1 ! λ i µ i .
The idea now is to determine the stationary distribution π(ε), making use of the results of [2, Section 3]. As established there, π(ε) is analytic around ε = 0, and as such admits a Taylor-series expansion:

π(ε) = ∞ k=0 ε k π (k) , (22)
for appropriate (vector-valued) coefficients π (k) . Before pointing out the main result of this subsection, we first introduce a number of matrices; D is the (d × d) deviation matrix corresponding to Q, as before:

• V and W are d × (C + 1)d and (C + 1)d × d matrices, respectively, defined by

V :=     ν T 1 . . . ν T d     ; W :=     e . . . e     = I ⊗ e,
• and D represents the deviation matrix for Q(0), with Di the deviation matrix of Qi ,

D =     D1 . . . Dd     .
Then the results of [START_REF] Avrachenkov | The first Laurent series coefficients for singularly perturbed stochastic matrices[END_REF] imply that the vectors π (k) have the geometric structure given by

π (k) = π (0) Q (1) D(I + Q (1) W DV ) k .
This expression can be simplified since, after some manipulations, it turns out that DQ (1) W DV vanishes in our specific case. Indeed, as Q (1) = (Q ⊗ I) and W = I ⊗ e, we have

DQ (1) W =     D1 . . . Dd     • (Q ⊗ e) = 0,
where the last step is due to the fact that Di e = 0. We obtain the following result.

Proposition 3.2. Consider the scaling (λ, µ, εQ, C). As ε ↓ 0, π(ε) admits the Taylor-series expansion (22), where

π (k) = π (0) Q (1) D k .
We now analyze the complexity of computing π (k) relying on Prop. 3.2. Observe that (in selfevident notation) π (0) i = α i ν i , which is in line with the intuition that when the switching between different background states is very slow, the stationary distribution is effectively reached within each regime. Note that, due to the fact that the matrices Qi are tridiagonal, computing products of the form v Di is an operation that has a computational complexity of O(C), and hence computing an additional term π (k) takes O(Cd) operations. The multiplication with Q (1) takes O(Cd 2 ) operations, and is therefore the bottleneck. It means that when K terms of the expansion (22) are to be evaluated, the complexity is O(d 3 +KCd 2 ). Our experiments (reported in the next section) indicate that the evaluation of π (k) through Prop. 3.2 typically compares favorably with the complexity of the approaches discussed in Remark 2.2 (which were O(Cd 3 ) or higher).

The above results lead to the following approximation in the slow regime. Corollary 3.2. Consider the scaling (λ, µ, εQ, C). As ε ↓ 0, with π(ε) admitting the Taylor-series expansion (22), and π (k) as in Prop. 3.2, the following Taylor expansion applies to p bl :

p bl = 1 αλ d j=1 π C,j (ε)λ j .

Numerical experiments

In this section we report on the output of a set of numerical experiments. Section 4.1 studies the impact of the model's parameters. In Section 4.2 we perform extensive accuracy tests for a broad range of parameter settings.

Figure 1: log p bl as a function δ ∈ [0, 1], for three values of q. The top curve corresponds to q = 1 10 , the middle one to q = 1, and the bottom one to q = 10.

Impact of model parameters

In this subsection we consider the following experiments:

• We first demonstrate by means of a concrete example that modulation can result in an increase as well as a decrease of the blocking probability.

• Next, we present an example that is indicative of the potential errors that would result from replacing a Markov-modulated loss system by its non-modulated counterpart, relying on the expressions obtained in Section 2. This example underscores the relevance of the model studied in this paper; naïvely using the 'classical' Erlang loss model typically leads to substantial errors.

• For the fast regime, we then study the numerical performance of approximations based on the asymptotic results derived in Section 3. Importantly, we derive a new approximation of the blocking probability that is valid across all f > 0, with the situations f < 1, f > 1, and f = 1 not needing to be distinguished, which outperforms ( 19)-( 20).

• We conclude the section with experiments for the slow regime, showing a high degree of accuracy for small ε.

A. Modulation can affect blocking probability both ways. Following the intuition that modulation increases the variance and more variance means worse performance, one might be led to think that adding modulation (while keeping constant) will always increase the blocking probability. We show that this is not the case and that the effect can go both directions. To this end, consider a system with C = 5 and a three-state background chain with generator Q defined as follows:

Q =    -10 5 5 1 -2 1 1 8 -9    .
• If we take the arrival and service rates constant (thus effectively ending up in an unmodulated system) at λ 1 = λ 2 = λ 3 = 2.0727 and µ 1 = µ 2 = µ 3 = 1, we find an offered load of = 2.0727, and a blocking probability of 0.0409.

• With the same C and Q but with λ 1 = 12, λ 2 = 4, λ 3 = 8, µ 1 = 3, µ 2 = 5, and µ 3 = 1, we find that is still 2.0727, but the blocking probability drops to 0.031.

• Alternatively, with the same C and Q, but with λ 1 = 0.6820, λ 2 = 2.0727, λ 3 = 3, and µ 1 = µ 2 = µ 3 = 1, we find again the same , but the blocking probability is now 0.0422, which is higher than in the unmodulated case.

B. Impact of burstiness on blocking probability. In this example, with d = 2, we assume q 12 = q 21 =: q, so that α 1 = α 2 = 1 2 , and µ 1 = µ 2 = 1. The number of servers C is set to 100. We create burstiness by assuming heterogeneity in the arrival rates: λ 1 = 70 + 50 δ and λ 2 = 70 -50 δ, with δ ∈ [0, 1]. Observe that (i) the system is underloaded on average, in the sense that := (α 1 1 + α 2 2 )/C = 0.7 < 1 (with i := λ i /µ i ), irrespective the value of δ, (ii) for δ = 0 the system is an ordinary (nonmodulated) Erlang-loss system operating in underload (in that /C < 1), (iii) for δ = 1 the system alternates between 'overloaded periods' (as 1 /C = 120/100 > 1) and 'underloaded periods' (as When evaluating the performance of this system, one could naïvely use the classical (i.e., nonmodulated) Erlang-loss formula, neglecting the burstiness in the arrival process. This approach evidently gives the correct result for δ = 0, but becomes increasingly questionable when δ grows. Fig. 1 shows log p bl as a function of δ ∈ [0, 1], for three values of q. We observe that the blocking probability significantly grows with δ. This effect is not that strong for large q, due to the fact that the modulating process is relatively fast: the periods with overload are so short that hardly any blocking occurs before the modulating process jumps back to the underloaded state. For small q, however, the increase of the blocking probability is rather pronounced.

C. Fast regime, convergence to truncated Normal distribution. In the second series of experiments we study the convergence of the random variable M (N ) to (N(0, σ2 ) | N(0, σ 2 ) β) as N tends to ∞. We choose β = 0.5, and

λ = 1.2 0.6 , µ = 0.6 1.8 , Q = -1 1 2 -2 .
The QQ-plots in Fig. 2 illustrate the convergence as N becomes large, which is particularly fast for larger values of f .

m := E[N(0, σ 2 ) | N(0, σ 2 ) β] = -σ φ N (η) Φ N (η) , v := Var[N(0, σ 2 ) | N(0, σ 2 ) β] = σ 2 1 - φ N (η) Φ N (η) η + φ N (η) Φ N (η) .
We use the same values for Q, λ and µ as in the previous example. The numerical experiments show that v N := Var M (N ) indeed converges to v, but that this convergence is rather slow: compare N -2γ Var M (N ) in the top panels of Fig. 3 (the solid curves) with its theoretical limit (the blue step function, with the isolated point corresponding to f = 1).

As we concluded above, the approximation Var M (N ) ≈ N 2γ v typically performs poorly. Interestingly, however, there is a highly-accurate alternative. The idea is to replace σ 2 by its 'prelimit counterpart'

σ 2 N := 1 µ ∞ N 2-f -2γ U + N 1-2γ λ ∞ = N 2-f -2γ U µ ∞ + N 1-2γ , (23)
as in Equation ( 27) in Appendix A; likewise, we also replace η by η N := β/σ N . This leads to the approximation:

Var M (N ) ≈ N 2γ v N , with v N := σ 2 1 - φ N (η N ) Φ N (η N ) η N + φ N (η N ) Φ N (η N ) .
The plots in the bottom panels of Fig. 3 show that the approximation is remarkably accurate: the solid and dashed curves are close. Important from an application perspective is that the approximations are conservative.

E. Fast regime, approximation of the blocking probability. Here we see behavior that is very similar to that observed for the variance of M (N ) , in the sense that (i) the convergence ( 18) turns out to be slow, and (ii) replacing σ 2 by σ 2 N (as given by Equation ( 23)) drastically improves the performance. As a result, we have the accurate approximation

p (N ) bl ≈ N γ-1 b N , with b N := σ N φ N (η N ) Φ N (η N ) .
The results displayed in Fig. 4 again show a nearly perfect fit for f > 0.5, and still reasonably accurate performance for f ∈ (0.1, 0.5); again, we have used the same values for Q, λ and µ.

Based on the above observations, we therefore propose the replacement of the approximations ( 19) and ( 20) for the blocking probability in the non-scaled model (λ, µ, Q, C) for f > 1 and f < 1 by the approximation ( 21) that corresponds to f = 1 across all values of f > 0. Throughout our experiments, we have observed that this approximation is typically highly-accurate, and when there is a bias it tends to be conservative.

p bl ≈ U/µ ∞ + • φ N C - U/µ ∞ + Φ N C - U/µ ∞ + . (24)
F. Slow regime, convergence to the limiting variance and approximation of the blocking probability. This last example considers the approximation proposed for the slow regime. We conduct an experiment showing the good fit even a lower-order approximation can provide. With d = 2, we assume q 12 = q 21 = 1, so that α 1 = α 2 = 1 2 . Also, we take µ = (1, 2), and λ = (2, 1), and as a consequence the offered load is equal to 1. We consider the scaling (nλ, µ, εQ, 1.2 n) and let ε be small; we take n = 500 and n = 5 000. In Figs. 5-6 we consider the blocking probability p bl and the variance of the stationary number of customers Var M , by truncating the series (22) after the quadratic term.

First observe that the blocking probability is decreasing in ε, which can be understood as follows. The parameters are chosen such that the system is a temporarily underloaded (overloaded) when in state 1 (state 2) respectively. When ε increases, the background process jumps more frequently between the states, thus reducing the burstiness, and hence the blocking probability. A similar reasoning applies to the variance. Figs. 56show that for relatively small ε the approximation is accurate, even though just three terms, a constant, a linear, and aquadratic term have beeen included. The fit can be made more accurate for larger ε by taking into account a quadratic term.

Validation of approximation

In this subsection we discuss the performance, under a broad range of possible parameter values, of the approximation (23) of the variance of M (N ) , and the approximation (24) of p bl . Recall that these approximations rely on our result that the M (N ) converges to a truncated Normal random variable. As this convergence can be made arbitrarily slow (just as this can be done in a conventional central limit theorem setting), it is clear that it is possible to construct cases in which the approximations lose accuracy. More specifically, when (i) the transition rates of the modulating Markov chain are low, or when (ii) the arrival rates (and/or service rates) during different states are highly heterogeneous, we anticipate that the accuracy of our approximation is likely to degrade. We shall see that this intuition is confirmed by the experiments below. In view of the above observations, we performed the following extensive set of numerical experiments. We chose thoroughly to examine the case where the modulating Markov chain has two states, fixing the hedge β to some values that we will detail below, and fixing f to 1 (as any other value of f can be obtained by rescaling the parameters). This brings us to six free parameters, namely λ 1 , λ 2 , µ 1 , µ 2 , q 1 := q 12 , and q 2 := q 21 .

We bring this down to four by imposing the constraint that the average arrival rate α 1 λ 1 + α 2 λ 2 , as well as the average service rate α 1 µ 1 + α 2 µ 2 , must be equal to 1. As we scale up time and space afterwards, these constraints in fact do not restrict the parameter space, it just makes it so that for the same scaling parameter N , and identical β and f , we have the same service capacity (because the loads are the same), and the same 'relaxation rate' µ ∞ (which governs the rate of convergence to equilibrium in the model's infinite-server counterpart; see for example [START_REF] Blom | Analysis of Markov-modulated infinite-server queues in the central-limit regime[END_REF]). This leaves us with four (positive) parameters: q 1 , q 2 , and the ratios r 1 := λ 1 /λ 2 and r 2 := µ 1 /µ 2 . In order to explore this four dimensional parameter space, we randomly generated these parameters in the following way. Let ν 1 , . . . , ν 4 denote four standard normal random variables and c 1 , . . . , c 4 four constants. Then the random parameters Q 1 , Q 2 , R 1 , R 2 were given by

Q 1 = 10 c 1 +ν 1 ; Q 2 = 10 c 2 +ν 2 ; R 1 = 10 c 3 +ν 3 ; R 2 = 10 c 4 +ν 4 .
Hence, the parameters essentially adhere to a log-normal law. Choosing the constants c i appropriately allow us to explore more extreme areas of the parameter space.

For each row of the two tables in Appendix B, we generated 50 000 points in parameter space. We list the constants c i , β, and four characterizations of the relative error of (A) the standard deviation of the system content and (B) the blocking probability (indicated in the following by a superscript '(s)' and '(b)', respectively): with ∈ {(s), (b)},

• p 5 and p 10 denote the fraction of samples that have a relative error of more than 5% resp. 10%, and

• E and σ denote the average relative error and its standard deviation.

The first table shows the results for scaling parameter N = 500, the second for N = 5 000. The general conclusion is that the approximations perform well with our procedure being set up to cover a broad range of parameter values (some of which rather extreme), typically a relatively low fraction of scenarios leads to substantial inaccuracies. In addition, the tables show that even if the low value of N brings already quite satisfactory results, they improve significantly when N is increased to 5 000. We also observe an asymmetric effect when the mean of the Q 2 samples is shifted: we see an improvement when the values are on average higher, and a marked deterioration when the values are lower. This is as could have been expected, as the convergence of the background chain will be respectively higher and lower in these cases. For the other two parameters, shifts in both directions lead to loss of accuracy. The quality of the approximation of the standard deviation generally improves when β increases, whereas the opposite is true when the blocking probability is concerned. This can be explained from the fact that for higher β the blocking probability requires convergence of the tail of the distribution.

Concluding remarks

In this paper, we have presented a full analysis of the Markov-modulated counterpart of the classical Erlang loss system. Our first main result is the derivation of an explicit formula describing the stationary behavior, which also leads to a closed-form expression for the blocking probability. As a second contribution we provide simple, yet highly-accurate approximations which are provably correct under specific scalings.

The approximations relate to a fast regime (in which arrivals as well as the background process' transitions are sped up) and a slow regime (in which the background process jumps slowly). In the fast regime, the steady-state distribution tends to that of a truncated Normal random variable. Potential topics for future research are the following:

• Large deviations of the blocking probability under the scaling (N λ, µ, N f Q, N C), where the capacity is blown up proportionally to the arrival rate, whereas the transitions of the background process can be relatively slow (f < 1) or fast (f > 1). We anticipate that the blocking probability will decay essentially exponentially. Recent work on similar models [START_REF] Huang | Large deviations for Markov-modulated diffusion processes with rapid switching[END_REF] has shown that the case f = 1 is delicate and should be handled separately.

• The development of sound staffing rules for the Markov-modulated loss system. A similar issue could be studied for the model in which customers finding all servers busy join a waiting line.

A Sketch of the proof of Proposition 3.1

First observe that

(P (N ) ) (s) = - 1 N γ C N k=0 k exp -s k -N N γ π (N ) k + N 1-γ P (N ) (s),
or, alternatively,

C N k=0 k exp -s k -N N γ π (N ) k = N P (N ) (s) -N γ (P (N ) ) (s). (25)
Multiplying the balance equations of the scaled model by exp(-sN -γ (k -N )) and summing over k ∈ {0, . . . , C N } yields the identity

C N k=0 exp -s k -N N γ π (N ) k N diag{λ}1 {k<C N } + k • diag{µ} -N f Q = C N k=0 exp -s k -N N γ π (N ) k-1 N diag{λ}1 {k>0} + C N k=0 exp -s k -N N γ π (N ) k+1 diag{µ}1 {k<C N } . (26)
be concluded that c 1 = 1 due to the lack of the non-homogeneous term; evidently, we cannot do this here. The next (standard) step is to try g(s)h(s) as solution to the non-homogeneous limiting differential equation, with the function h(•) to be identified. To this end, we first define

κ := lim N →∞ N 1-γ (κ N diag{λ}e)/µ ∞ ,
whose value we determine below. We obtain the equation

g (s)h(s) + g(s)h (s) = sg(s)h(s)σ 2 + κe -sβ ,
such that h (s) = κe -sβ /g(s), and hence

h(s) = s -∞ κ c 1 e -r 2 σ 2 /2-rβ dr + c 2 .
We thus have found the solution

ψ(s) = c 1 e s 2 σ 2 /2 • κ c 1 s -∞ e -r 2 σ 2 /2-rβ dr + c 2 .
The next goal is to identify the unknown constants. To this end, notice that

s -∞ e -r 2 σ 2 /2-rβ dr = sσ -∞ 1 σ exp - 1 2 u + β σ 2 e β 2 /(2σ 2 ) du,
which can be interpreted as P (N(0, 1) < σs + β/σ) σ -1 e β 2 /(2σ 2 ) . Using the standard asymptotic relation P(N(0, 1) > x) xe x 2 /2 √ 2π → 1, this quantity behaves, for s → -∞, as

1 √ 2π 1 σs + β/σ e -(σs+β/σ) 2 /2 1 σ e β 2 /(2σ 2 ) = 1 √ 2π 1 σ 2 s + β e -βs e -s 2 σ 2 /2 .
We conclude that

lim s→-∞ 1 s log e s 2 σ 2 /2 • 1 c 1 s -∞ e -r 2 σ 2 /2-rβ dr = β.
Because M (N ) β, we also know that, for s < 0, ψ(s) e -sβ , and hence

-lim s→-∞ 1 s log ψ(s) β,
and therefore c 2 = 0. Now from ψ(0) = 1 it follows that

κ = 0 -∞ e -r 2 σ 2 /2-rβ dr -1 = σ φ N (β/σ) Φ N (β/σ) .
where we recall that φ N (•) and Φ N (•) denote the density and cumulative distribution function of a standard Normal random variable, respectively. It takes elementary calculus to verify that the resulting expression for ψ(s) coincides with ( 16).

B Output of numerical experiments

[c 1 , c 2 , c 3 , c 4 ] β p (s) 5 p (s) 10
E (s) σ (s) p (b) 5 p (b) 10
E (b) σ (b) [0, 0, 0, 0] 0.5 7

Figure 2 :

 2 Figure 2: QQ-plots of M (N ) versus (N(0, σ 2 ) | N(0, σ 2 )β). The dashed curve corresponds to N = 100, and the solid one to N = 500; to assess the convergence, the diagonal (in blue) has also been depicted.

Figure 3 :

 3 Figure 3: The solid curves are N -2γ Var M (N ) as a function of f , for four values of β (the same in the top graphs and the bottom graphs). In the top graph N -2γ Var M (N ) can be compared with v (in blue), and in the bottom graphs with v N (dashed curves). The red curves correspond to N = 500 and the black curves to N = 5000.

Approximation 4 . 1 .

 41 Consider the non-scaled model (λ, µ, Q, C). Then, with C > ,

Figure 4 :

 4 Figure 4: The solid curves are N 1-γ p (N ) bl as a function of f , for four values of β (the same in the top graphs and the bottom graphs). In the top graph N 1-γ p (N ) bl can be compared with b (in blue), and in the bottom graphs with b N (dashed curves). The red curves correspond to N = 500 and the black curves to N = 5000.

Figure 5 :

 5 Figure 5: p bl as a function of ε ∈ (0, 1), with the quadratic quasi-stationary approximation indicated in green.

Figure 6 :

 6 Figure 6: The variance Var M as a function of ε ∈ (0, 1), with the quadratic quasi-stationary approximation indicated in green.

Table 1 :

 1 Results for N = 500. [c 1 , c 2 , c 3 , c 4 ]

			.75%	3.23% 1.95 (-02) 6.61 (-02)	8.34%	3.10%	2.74 (-02) 4.49 (-02)
	[0, 0, 0, 0]	1.0	6.67%	2.66% 1.74 (-02) 7.38 (-02)	6.92%	3.04%	1.33 (-02) 4.28 (-02)
	[0, 0, 0, 0]	1.5	5.95%	2.31% 1.60 (-02) 4.26 (-02)	6.31%	2.78% -1.03 (-02) 6.35 (-02)
	[0, 0, 0, 0]	2.0	4.84%	2.03% 1.46 (-02) 4.69 (-02) 83.16%	3.89% -4.98 (-02) 9.44 (-02)
	[0, 2, 0, 0]	0.5	0.54%	0.14% 7.59 (-03) 8.70 (-03)	0.51%	0.10%	1.71 (-02) 6.32 (-03)
	[0, 2, 0, 0]	1.0	0.33%	0.08% 7.38 (-03) 9.77 (-03)	0.61%	0.23%	1.73 (-03) 1.01 (-02)
	[0, 2, 0, 0]	1.5	0.49%	0.15% 7.65 (-03) 9.46 (-03)	0.26%	0.07% -2.67 (-02) 9.15 (-03)
	[0, 2, 0, 0]	2.0	0.29%	0.08% 6.57 (-03) 5.70 (-03) 97.69%	0.29% -7.23 (-02) 2.01 (-02)
	[0, 0, 2, 0]	0.5 14.39%	6.77% 3.15 (-02) 1.02 (-01) 12.73%	4.67%	3.24 (-02) 4.65 (-02)
	[0, 0, 2, 0]	1.0 12.28%	5.82% 2.76 (-02) 8.40 (-02) 10.65%	4.20%	1.90 (-02) 4.80 (-02)
	[0, 0, 2, 0]	1.5 11.28%	5.38% 2.58 (-02) 7.93 (-02)	9.51%	4.22% -2.60 (-03) 5.78 (-02)
	[0, 0, 2, 0]	2.0 10.18%	4.88% 2.51 (-02) 9.95 (-02) 71.63%	4.54% -3.75 (-02) 7.85 (-02)
	[0, 0, 0, 2]	0.5 15.55%	7.62% 3.33 (-02) 8.02 (-02) 11.70%	4.67%	3.17 (-02) 4.81 (-02)
	[0, 0, 0, 2]	1.0 14.11%	6.15% 2.92 (-02) 7.18 (-02)	9.38%	4.00%	1.71 (-02) 4.99 (-02)
	[0, 0, 0, 2]	1.5 11.73%	5.42% 2.73 (-02) 7.75 (-02)	8.25%	4.03% -5.36 (-03) 7.01 (-02)
	[0, 0, 0, 2]	2.0 10.85%	4.60% 2.51 (-02) 6.94 (-02) 76.54%	5.25% -4.26 (-02) 9.36 (-02)
	[0, -2, 0, 0]	0.5 21.90% 13.95% 7.31 (-02) 3.89 (-01) 28.24% 18.46%	8.40 (-02) 2.15 (-01)
	[0, -2, 0, 0]	1.0 20.78% 13.27% 6.93 (-02) 3.74 (-01) 26.20% 18.44%	7.93 (-02) 2.36 (-01)
	[0, -2, 0, 0]	1.5 18.80% 12.13% 6.29 (-02) 4.20 (-01) 23.98% 16.92%	5.89 (-02) 2.96 (-01)
	[0, -2, 0, 0]	2.0 17.73% 11.64% 6.40 (-02) 2.84 (-01) 86.37% 23.07%	5.68 (-02) 5.83 (-01)
	[0, 0, -2, 0]	0.5 14.21%	6.49% 3.05 (-02) 1.24 (-01) 12.43%	4.31%	3.19 (-02) 4.82 (-02)
	[0, 0, -2, 0]	1.0 12.50%	5.66% 2.65 (-02) 7.28 (-02) 10.94%	4.08%	1.89 (-02) 4.55 (-02)
	[0, 0, -2, 0]	1.5 11.18%	5.18% 2.58 (-02) 9.08 (-02)	9.67%	4.65% -1.47 (-03) 6.18 (-02)
	[0, 0, -2, 0]	2.0 10.37%	4.97% 2.41 (-02) 8.16 (-02) 72.37%	5.23% -3.51 (-02) 9.33 (-02)
	[0, 0, 0, -2]	0.5 16.73%	7.47% 3.44 (-02) 1.16 (-01) 12.01%	4.68%	3.23 (-02) 5.80 (-02)
	[0, 0, 0, -2]	1.0 14.46%	6.61% 3.00 (-02) 7.67 (-02)	9.23%	3.85%	1.72 (-02) 5.71 (-02)
	[0, 0, 0, -2]	1.5 12.36%	5.67% 2.80 (-02) 7.25 (-02)	8.18%	4.10% -4.64 (-03) 8.30 (-02)
	[0, 0, 0, -2]	2.0 10.58%	4.67% 2.58 (-02) 1.12 (-01) 77.73%	5.41% -4.19 (-02) 8.85 (-02)

Table 2 :

 2 Results for N = 5 000.

/C = 20/100 < 1). We conclude that the burstiness in the arrival process increases when δ goes from 0 to 1; δ can therefore be interpreted as a 'burstiness parameter'.

Acknowledgements

The authors thank S. Asmussen and the anonymous referees for useful comments.

M. Mandjes' research is partly funded by the NWO Gravitation project NETWORKS, grant number 024.002.003. P.G. Taylor's research is supported by the Australian Research Council (ARC) Laureate Fellowship FL130100039 and the ARC Centre of Excellence for the Mathematical and Statistical Frontiers (ACEMS). This research was partly performed when K. De Turck was a Postdoctoral Fellow of Fonds Wetenschappelijk Onderzoek / Research Foundation -Flanders.

Now let us study (26) in greater detail. Writing κ N := π (N )

C N , the first term on the left hand side can be rewritten as N P (N ) (s) -κ N e -sβ diag{λ};

the first term on the right hand side is precisely equal to this expression, but multiplied by a factor exp(-sN -γ ). Using (25), the second term on the left hand side reads

and the second term on the right hand side is equal to this expression multiplied by exp(sN -γ ).

We thus arrive at the system of differential equations

The next steps are now exactly as in [START_REF] Blom | Analysis of Markov-modulated infinite-server queues in the central-limit regime[END_REF], and we therefore omit details here. We postmultiply the equation by the fundamental matrix F := D +A (where A := eα), use QF = I -A, bring P (N ) (s)A to the right hand side, to obtain an expression for P (N ) (s) in terms of itself and its derivative. Then we use the Taylor-series expansion

to iterate the resulting identity, which we then postmultiply with eN f . Mimicking the reasoning in [START_REF] Blom | Analysis of Markov-modulated infinite-server queues in the central-limit regime[END_REF], one obtains, with ψ (N ) (s) := P (N ) (s)e,

where the first term on the right hand side cancels (as an immediate consequence of the definition of ). After dividing by µ ∞ s, we obtain the one-dimensional ordinary differential equation,

where U , defined in [START_REF] Miyazawa | Insensitivity and product form decomposability of reallocatable GSMP[END_REF], is given by

Now we consider the limiting solution ψ(•) when N → ∞, first looking at the homogeneous version of the differential equation ( 27) which is precisely the same as the equation in [START_REF] Blom | Analysis of Markov-modulated infinite-server queues in the central-limit regime[END_REF]). By distinguishing between the cases f < 1, f > 1, and f = 1, the reasoning presented in [START_REF] Blom | Analysis of Markov-modulated infinite-server queues in the central-limit regime[END_REF] yields the solution g(s) := c 1 exp(s 2 σ 2 /2), with σ 2 as defined above, for some free constant c 1 . In [START_REF] Blom | Analysis of Markov-modulated infinite-server queues in the central-limit regime[END_REF] it could


Mathias Quillot

Cassandre Ollivier

Richard Dufour

Vincent Labatut

Exploring Temporal Analysis of Tweet Content from Cultural Events

Keywords: Word embedding, Temporal representation, Statistical analysis, Cultural events

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

Introduction

Social networks have become a new way of communicating and sharing information and views that can be accessed by billions of people. These social interactions may take the form of short text messages, such as the Twitter platform, where users have the possibility to instantly send a message (here a tweet) containing around 140 characters.

Thanks to its ease of use, Twitter is currently an essential platform for the exchange of messages. For particular events (news, concerts, festivals, presidential elections, etc.), users are increasingly inclined to express themselves through these short messages. Although this platform has become a formidable object of study for a variety of domains ranging from sociology [START_REF] Java | Why we twitter: understanding microblogging usage and communities[END_REF][START_REF] Romero | Differences in the mechanics of information diffusion across topics: idioms, political hashtags, and complex contagion on Twitter[END_REF][START_REF] Murthy | Towards a sociological understanding of social media: Theorizing twitter[END_REF] to automatic information extraction [START_REF] Endarnoto | Traffic condition information extraction & visualization from social media twitter for android mobile application[END_REF][START_REF] Sriram | Short text classification in twitter to improve information filtering[END_REF][START_REF] Nakov | Semeval-2016 task 4: Sentiment analysis in twitter[END_REF], the short format of the messages and the large size of the corpora often both make them difficult to analyze. In this article, we analyze messages exchanged through the Twitter platform in the context of cultural events, with a particular focus on festivals. More precisely, we seek to account for shared content, through the words contained in tweets. The major difficulty of this type of analysis lies in the duration of the considered events: although a festival takes place over a defined period (ranging from a few days to several weeks), the activity of users on social networks can intervene at any time (before, during, or after the festival).

We assume that it is difficult to reveal all the information conveyed through the discussions (the tweets) in a global way without taking into account the temporal aspect of the messages (i.e. their emission date). A global model would have a tendency to reveal only the frequently shared information, ignoring the uncommon ones that could nonetheless be important over a particular period of time. This is problematic when these models are used as input for information retrieval tasks, such as automatic event extraction, automatic summarization, etc. Based on this observation we propose a preliminary work close to those initiated in [START_REF] Basile | Analysing word meaning over time by exploiting temporal random indexing[END_REF][START_REF] Hamilton | Diachronic word embeddings reveal statistical laws of semantic change[END_REF] that seeks to compare two word embedding-based models: one ignoring the temporal aspect of the messages, using the state-of-the-art Word2Vec model [START_REF] Mikolov | Efficient estimation of word representations in vector space[END_REF], and one taking advantage of the emission date of tweets, using the Temporal Embedding approach.

The rest of this article is organized as follows. Section 2 presents the different methods used to analyze the impact of time in the events analysis. We presentthe experimental protocol as well as the results obtained in Section 3. Finally, we conclude and give some perspectives in Section 4.

Methods

We seek to highlight the interest of taking into account the temporal information conveyed by messages emitted through social networks, in a context of cultural events analysis. In Subsection 2.1, the two compared word embedding representations are described: one considers the complete set of messages regardless of their emission date (the Word2Vec neural network method), while the second one takes into account the chronology of the documents (Temporal embedding approach). In Subsection 2.2, we describe the methods used to compare the word embedding models, which include both subjective and objective tools. The goal here is to identify the points on which the models differ or converge. Figure 1 summarizes our overall framework, and is detailed in the rest of this section. Our different models are specified there in bold.

We note N the total number of unique words in the corpus. When performing the analysis, one can focus on a list of n words of interest corresponding to a subset of the corpus lexicon: this allows the end user to adopt either a verification-directed approach. If the user has some a priori knowledge and would like to check certain assumptions regarding the corpus, or a exploratory approach, consisting in using the whole lexicon as the word list (in this case, n = N ).

Visual comparison Statistical comparison

Visual comparison

Corpus tweets

Word2Vec

Word Filtering

Month

Word Embedding Representations

We briefly describe the classic Word2Vec model, before explaining how the Temporal Embedding model is extracted. Two distinct temporal resolutions are considered in this work (weeks vs. months).

Word2Vec Neural Network. Word2Vec models [START_REF] Mikolov | Efficient estimation of word representations in vector space[END_REF] are based on the hypothesis that semantically similar words tend to have similar contextual distributions.

Concretely, this context is a window whose size is expressed in a number of words. This is centered on the word of interest. For our experiments we use the CBOW method, it seeks to predict the word reference given a context. For instance, a context of 2 words, the CBOW neural network model takes an input taking the form of a sequence of 4 words w i-2 , w i-1 , w i+1 and w i+2 , and outputs a word w i . We only use the hidden layer of the neural networks, which means each word is represented by a vector. The length of this vector is specified by the user as a parameter d which is 200 by default in literature. We use the complete word vocabulary to train the CBOW model, the method outputs an N × d matrix. More information about Word2Vec models can be found in [START_REF] Mikolov | Efficient estimation of word representations in vector space[END_REF].

Temporal Embedding. Instead of globally taking all the corpus words into account for this method, as for the Word2Vec model, we focus on a predefined list of n words of interest (represented by the Word filtering box in Figure 1). For each word in this list we count its number of occurrences by time unit where the time unit is either one month or one week, depending on the considered model. In order to avoid zero values, which can be frequent when considering weekly occurrences, we smooth the numbers of occurrences through a moving average. This results in an n × m occurrence matrix called Temporal matrix, where m is the period covered by the corpus expressed in time units. We then perform a Principal Component Analysis (PCA) on this matrix. This provides us with the temporal embedding: another n × m matrix, whose columns are the obtained components ordered by increasing informativeness. The interest of the PCA is to get a more compact representation of the temporal embedding by focusing on the first few components.

Model Comparison

Comparing models is usually a difficult problem, the most frequent solution being to compare their performance on a targeted task (for example, in speech recognition, the best models are those that allow the lowest word error rate).

In this paper, we investigate the interest of including temporal information in word representations from cultural events. To correctly evaluate the impact of each model, an objective ground truth would be necessary, i.e. knowing which words clearly represent a particular event and should be associated with it (and therefore which words have no interest). This would reveal the model that best represents a targeted cultural event. Since no ground truth is available, we first perform an objective comparison relying on statistical tests. These results are generally considered more reliable (and more easily comparable because experiments are reproducible). We seek to know with this first evaluation how close or different the models are. Then, we perform a more detailed analysis through a subjective comparison based on a visual human interpretation in order to investigate the contribution of the temporal information to cultural event representation.

Objective comparison. Two statistical tests are used to compare the models globally: Wilcoxon-Mann-Whitney and Kendall's τ . These non-parametric tests check the hypothesis whether two samples originate from the same distribution. They are complementary, in the sense the former can be considered as a medianbased version of the t-test, whereas the latter focuses on rank correlation. We apply them to the Rescaled Correlation and Cosine Distance matrices which were previously computed for the models when extracting the plots and dendrograms.

The tests allow us to compare the way the pairs of words are ordered in the different models based on their distances. In other words, they compare the models based on the relative positions of the words in the model spaces.

Besides this global comparison, we also perform a local one by focusing on each word separately. For a given word, we perform the same Kendall's and Wilcoxon-Mann-Whitney tests as before, on the distances between this word of interest and the other considered words. When comparing two models, we ultimately identify two groups of words: those whose relative positions are significantly different in these models (using a significance level of .05), and the others, whose relative positions are supposedly similar in both models.

Subjective comparison. The second comparison is visually performed by humans, based on graphical representations of the models. We consider two of them: 1) the projection of the words in a 2D space, and 2) dendrograms. The former representation allows us to identify opposition between words, whereas the second one focuses on their associations.

The 2D space representation is obtained by considering the first 2 components (i.e. the most informative ones) of a PCA. Since the Temporal Embedding model includes a PCA, no additional process is required. For the Word2Vec approach, some additional processing is needed: we extract the n rows corresponding to the word list from the N ×d matrix (this step is represented as the Word filtering box in Figure 1), resulting in the n × d so-called Filtered Matrix, on which a PCA is performed. The subjective comparison is conducted by checking how the words are spatially separated by the plot axes and how this differs from one model to the other. Put differently, if an axis separates two pairs of words and these words are away from this axis, we will consider that they are in opposition. If an opposition is present in both models, we consider this as a similarity between the models. On the contrary, if an opposition is found in one model only, we consider it as a difference between the models.

The dendrograms allow us to identify how the models gather words and organize them hierarchically. We obtain them using the standard hierarchical clustering algorithm available in the R Language. Note that this implementation requires a dissimilarity matrix as its input. Moreover, we use the complete linkage approach in order to favor compact clusters with small diameters. In the case of the temporal embedding, we first build an n × n correlation matrix based on the (week/month) temporal matrix (i.e. without the PCA). The correlation between two words is obtained by processing Pearson's coefficient for the two rows associated to these words in the matrix. When the correlation is negative, we set it to zero: this is a common practice when dealing with temporal series because in this context they are generally considered as noise. The dissimilarity is then obtained through the following rescaling: 1 -Cor, which in our case produces values ranging from 0 (similar) to 1 (dissimilar). This gives us an n × n matrix, which is called Rescaled Correlation Matrix in Figure 1. For the Word2Vec model, we build an n × n Distance Matrix based on the previously computed n×d filtered matrix. Each element of this distance matrix corresponds to the Cosine distance between two words of the list. Let us note w 1 and w 2 the respective word embeddings of these two words, then the distance is given by: d(w 1 , w 2 ) = 1 -Sim(w 1 , w 2 ), where Sim(w 1 , w 2 ) is the classic Cosine similarity between the words. Here the Cosine approach seems more appropriate than the Correlation used on the temporal data because we know that the rows of the considered matrix are inter-dependent by construction. After having generated the dendrograms, we make the visual comparison by checking if two words which are contiguous in one dendrogram are also placed together in another dendrogram.

Experiments

In this section, we briefly present the data analyzed in this study (Subsection 3.1), before describing and discussing the obtained results (Subsection 3.2).

Corpus and List of Words

We use the corpus provided for the MC2 CLEF 2017 lab1 which contains 70 million tweets [START_REF] Goeuriot | CLEF 2017 MC2 search and time line tasks overview[END_REF]. These were automatically retrieved from Twitter using a predefined set of keywords related to cultural festivals in the world. They cover a period ranging from May 2015 to November 2016 and are composed by 134 different languages [START_REF] Goeuriot | Overview of the CLEF 2016 cultural micro-blog contextualization workshop[END_REF][START_REF] Ermakova | Cultural micro-blog contextualization 2016 workshop overview: data and pilot tasks[END_REF].

We focus on a manually curated list of words of interest. It was originally designed by cultural sociologists to focus on festivals. We extended the list in the following way. Firstly, we added certain cities of interest based on various general and specialized sources: Wikipedia's List of the world's most liveable (31 cities) [START_REF]Wikipedia: List of the world's most liveable cities[END_REF], BFM Business's List of the top 20 European cities [START_REF]BFM Business: List of the top 20 european cities[END_REF], and festival cities from Wikipedia's List of theatre festivals (30 cities) [START_REF]Wikipedia: List of theatre festivals[END_REF], Red Bulletin's List of the top 15 music festivals (12 cities) [START_REF]Red Bulletin: List of the top 15 music festivals[END_REF], Sky Scanner's List of the top 10 music festivals [START_REF]Sky Scanner: List of the top 10 music festivals[END_REF], and Temps de Vivre's List of top cinema festivals [START_REF]Temps de Vivre: List of the top cinema festivals[END_REF]. Second, we added other words related to the concept of festival in general: "theater", "music", "film". Third, we added some commercial brands also related to festivals, such as apple or deezer (33 words). In total, the list contains 119 different words.

Results

In this subsection we compare two word embedding representations under the form of 3 distinct models2 . One Word2Vec model and two Temporal Embedding models based on two different time units (weeks and months respectively). We first discuss the outcome of the statistical tests before presenting the visual comparison of the plots and dendrograms (see Section 2).

Statistical methods. We apply both the Wilcoxon-Mann-Whitney and Kendall's tests on all 3 models: Month vs. Week Temporal Embeddings, Month Temporal Embedding vs. Word2Vec, and Week Temporal Embedding vs. Word2Vec. All tests return a p-value smaller than 10 -15 : for these implementations of the tests, this means that they always reject model independence. Kendall's τ , which is an association measure ranging from -1 to +1, is 0.68 when comparing the Month vs. Week Temporal Embeddings: this corresponds to a strong correlation between these models. For the Word2Vec vs. Month and Week Temporal models, we get τ = -0.06 for both comparisons, which means that Word2Vec is almost independent from our temporal models. According to these tests, the information encoded in the temporal models is not the same as the one conveyed by Word2Vec. Thus, they can be considered as complementary and are likely to lead to different results depending on the task at hand.

After the global analysis we switch to individual words to compare the models in terms of which words have significantly different relative positions in the two considered models or a similar position in both models. Based on Kendall's 2D representation. We now switch to the visual comparison starting with the 2D plot based on the 2 main components obtained from the PCA. The semantics of these components are not available for the Word2Vec model so we do not discuss them. We only consider the positions of word pairs in the graphical representation, and stress the presence of oppositions. We first compare the Month and Week Temporal Embedding models whose PCA are shown in Figures 2 and3, respectively. Globally, they seem to present the same oppositions. For instance, "Deezer" vs. "Venise" (in red in all the figures), and "Vancouver" vs. "Valdez" (in blue). This comparison is in line with our statistical results and seems to indicate that it is not necessary to use the week as a time unit because the month-based model requires less data and captures roughly the same information. We then compare both temporal models with the Word2Vec one, whose PCA is shown in Figure 4. Visually, the oppositions seem to differ more than previously. For instance, "deezer" and "Venise" are not opposed anymore, while "Vancouver" and "Valdez" are less opposed. However we also see some oppositions such as "Madrid" vs. "free" (in green), which are common to all 3 models. Let us focus on an illustrative example: the Apple Music Festival, which takes place in London. Both Temporal Models tend to group the 3 concerned words "London", "Apple" and "music" (represented in yellow in the figures), whereas Word2Vec does not. This means that temporal models tend to gather words from a same cultural event. If we examine finely the Word2Vec representation, we can observe it groups words by semantic category: there are several clusters of cities, whereas "Film", "Filmmakers", "Hollywood" and "movie" are together (cinematic items), and so are "musique", "opera", "theatre", "jazz" and "Bach" (musical items). This is of course consistent of how Word2Vec is supposed to work. In conclusion of this visual comparison, we can state that not only are the temporal and Word2Vec models different as shown by the statistical tests, but the chronological information encoded by the former also allows to identify relevant groups of words relatively to what we know of the studied corpus.

Dendograms. The dendrograms of the Month and Week Temporal models, which are represented in Figures 5 and6, respectively, look strongly similar. In particular, we get the same direct connection between "Kyoto" and "Albuquerque" (in blue in all the figures), or "Frankfurt" and "Teatro" (in purple). Note that they nevertheless differ on some pairs of words such as "Berne" and "Teatro" (in green), which are directly connected in the month model whereas "Berne" is connected to "hiphop" in the week model. Like before, there are visible differences between the representations of the temporal models and that of the Word2Vec, represented in Figure 7. Focusing on the same pairs of words, we see that "Kyoto" is not connected with "Albuquerque" anymore, and neither are "Frankfurt" and "Teatro". Moreover, "Berne" is neither connected with "Teatro" or "hiphop", but rather with "MontNaeba" (a Japanese mountain).

In the case of Word2Vec the words are grouped by semantic similarity. The hierarchical nature of the groupings seem to connect them according to hyper/hyponymy relationships. For instance, "rock" is connected to "heavy.metal" (in orange in the figure), the latter is a subgenre of the former. Both are close to "Jazz" in the dendrogram, itself connected to "reggae" (in red) which can be considered as musical style strongly influenced by certain forms of jazz. Elsewhere, "Tokyo" is connected to "Kyoto" (both Japanese cities, Kyoto is in blue), "Santa Maria" to "Los Angeles" (both Californian cities, in Cyan) or "Paris" and "Cannes" (both French cities, in yellow).

Unlike Word2Vec, the temporal models lead to clusters of words which are consistent with the fact that the corpus is related to festivals. For instance in the dendrograms of both temporal models, we observe that "Cannes" (French city hosting a Cinema festival) and "Hollywood" (Californian center of the movie industry) are directly connected (both in orange in the figures), they are close to "Film" and "movie" (themselves directly connected aas well and represented in red). We also have "Avignon" (French city) connected to "Santarcangelo di Romagna" (Italian city): both cities (in cyan in the figures) host a renown theater festival. An other example is "London", directly connected to "Music" and close to "Apple" (all in yellow): this would represent the previously mentioned Apple Music Festival in London. Generally speaking in the dendrograms, festival cities are associated to words which are semantically related to the nature of their festival.

The visual comparison of the dendrograms confirm and complete our previous observations. The Word2Vec model is definitely different from the temporal ones which are relatively similar to each other. The context-based approach of Word2Vec does not capture the temporal information conveyed by the other models. In both cases, words are grouped according to their relative semantics. But in the case of Word2Vec these groups are built on semantic proximity. However with the Temporal Embeddings words from the same groups are indirectly related through a festival. In conclusion, the latter models seem more appropriate to a festival-oriented analysis of this corpus.

Conclusion

In this preliminary work, our objective was to study how considering temporal information affects the word embedding-based modeling of text corpora. We built two Temporal Embedding models and one Word2Vec model based on a corpus of tweets focusing on cultural events. We studied and compared them objectively through statistical tests and visually through PCA plots and dendrograms. It turns out that both temporal models appear to be highly similar according to their PCA plots and dendrograms whereas they seem different from the Word2Vec model. The statistical tests conclude that when taken globally, they are all significantly different from on an other. However, when considering each word separately we show that they differ only on a small proportion of words albeit with a large magnitude. These first results show that temporal information is not completely captured on Word2Vec model and Temporal Embedding is worth considering, at least when dealing with analyzing temporally messages related to cultural events. However, a finer analysis is indispensable to correctly characterize the contribution of time information. To this end, we must carry on several steps. Firstly, we could adopt a more exploratory approach by expanding our analysis to the whole lexicon instead of focusing on a predefined list of words. Secondly, we could compare our models built on the considered corpus with out-of-the-box Word2Vec models in order to analyze the contribution of the corpus regarding how words are represented.

Fig. 2 .

 2 Fig. 2. Two first principal components of the Month Temporal Embedding model.

Fig. 3 .

 3 Fig. 3. Two first principal components of the Week Temporal Embedding model.

Fig. 4 .

 4 Fig. 4. Two first principal components of the Word2Vec model.

Fig. 5 .

 5 Fig. 5. Dendrogram of the Month Temporal Embedding model.

Fig. 7 .

 7 Fig. 7. Dendrogram of the Word2Vec model.

  Dendrogram of the Week Temporal Embedding model.

	Angra dos Reis Children Theatres Annecy Wight Teatro Frankfurt Fukuoka Athens Sydney Tbilisi deezer Barcelona trance Lisbonne Glastonbury Grahamstown Vienne movie Film Hollywood Cannes Filmmakers Paris Thespian Festival Santa Maria São Paulo Miami Auckland Bogotá dogg Garden Melbourne Curitiba spotify Berlin Byron Bay Valdez electronic Madrid Belgrade Adelaide Albuquerque Kyoto Amsterdam Genève kpop Performance Writing Music London apple Zurich early Toronto San Sebastián Lyon Perth Niagara-on-the-Lake Tehran Rio de Janeiro Deauville Nevada Venise Florianopolis Budapest Helsinki Boom Oslo Rennes Edinburgh Bergen Stockholm Mont Naeba Williamstown Santarcangelo di Romagna Avignon Costa Mesa Hamburg Calgary Londrina Vancouver jardin Teatri Salt Lake City Chicago Los Angeles celtic Tokyo musicians folk rock bach heavy metal reggae Salzburg rave Dublin Munich opera musique jazz Theatre hiphop Lincoln international International Theatre historic music Reggio Calabria Atlanta Sibiu New York City Moscow Brussels free Manchester showcasing Liverpool Berne Copenhague Luxembourg
	Fig. 6. Florianopolis Londrina Angra.dos.Reis Rio.de.Janeiro Curitiba São.Paulo Teatri Reggio.Calabria Santarcangelo.di.Romagna Valdez Costa.Mesa San.Sebastián Santa.Maria Los.Angeles Barcelona Madrid Bogotá Teatro Venise Deauville Avignon Lyon Rennes Paris Cannes Berne Mont.Naeba Genève jardin musique Lisbonne deezer kpop apple bach celtic rock heavy.metal jazz reggae hiphop electronic rave trance spotify dogg Glastonbury Wight Boom Amsterdam Miami Fukuoka Tokyo Kyoto Sibiu Tbilisi Budapest Frankfurt Vienna Hamburg Copenhagen Oslo Helsinki Belgrade Annecy Toronto Zurich Berlin Munich Tehran Moscow Stockholm Athens Luxembourg Brussels Children.Theatres Filmmakers Performance.Writing Music Garden Thespian.Festival Film Theatre International.Theatre Bergen Grahamstown Williamstown Hollywood Lincoln New.York.City Salt.Lake.City Vancouver Calgary Atlanta Chicago Auckland Byron.Bay Edinburgh Sydney Adelaide Melbourne Perth London Dublin Manchester Liverpool Albuquerque Nevada Niagara.on.the.Lake Salzburg opera showcasing musicians historic.music folk free early international movie

http://mc2.talne.eu/

http://tac.talne.eu

Acknowledgments

This work was funded by the GaFes project supported by the French National Research Agency (ANR) under contract ANR-14-CE24-0022.