Ivan Nourdin

Raghid Zeineddine
email: raghid.zeineddine@univ-lorraine.fr

An Itô's type formula for the fractional Brownian motion in Brownian time

Keywords: Fractional Brownian motion in Brownian time, change-of-variable formula in law, Malliavin calculus

, is by definition the process Z = X • Y . It is a continuous, non-Gaussian process with stationary increments, which is selfsimilar of index H/2. The main result of the present paper is an Itô's type formula for f (Z t ), when f : R → R is smooth and H ∈ [1/6, 1). When H > 1/6, the change-of-variable formula we obtain is similar to that of the classical calculus. In the critical case H = 1/6, our change-of-variable formula is in law and involves the third derivative of f as well as an extra Brownian motion independent of the pair (X, Y ). We also discuss briefly the case H < 1/6.

Introduction

If f : R + → R is C 1 then f (t) = f (0) + t 0 f ′ (s)ds for all t 0 whereas, if W is a standard Brownian motion and if f : R → R is C 2 then, by the Itô's formula,

f (W t ) = f (0) + t 0 f ′ (W s )d -W s + 1 2 t 0 f ′′ (W s )ds, t 0. (1.1)
In (1.1) the Itô integral, namely

t 0 X s d -Y s := lim n→∞ ⌊2 n t⌋-1 k=0 X k2 -n (Y (k+1)2 -n -Y k2 -n ), (1.2)
1 is of forward type. It is well-known that the additional bracket term 1 2 t 0 f ′′ (W s )ds appearing in (1.1) comes from the non-negligibility of the quadratic variation of W in the large limit; more precisely,

⌊2 n t⌋-1 k=0 (W (k+1)2 -n -W k2 -n ) 2 a.s.
-→ t as n → ∞.

(1.3)

Introducing a family {B H } H∈(0,1) of fractional Brownian motions parametrized by the Hurst parameter H may help to reinterpret (1.1) in a more dynamical way. Let us elaborate this point of view further. Recall that B 1 2 is nothing but the standard Brownian motion, whereas B 1 is the process B 1 t = tG, t 0, G ∼ N(0, 1). The extension of (1.3) to any H ∈ (0, 1) is well-known: one has

2 n(2H-1) ⌊2 n t⌋-1 k=0 (B H (k+1)2 -n -B H k2 -n ) 2 a.s.
-→ t as n → ∞.

(1.4)

Based on (1.4), it is then not difficult to prove the following two facts:

1. If H > 1 2 and f : R → R is C 2 (actually, C 1 is enough), then t 0 f ′ (B H s )d -B H s exists as a limit in probability and we have

f (B H t ) = f (0) + t 0 f ′ (B H s )d -B H s , t 0. 2. If H < 1 2 , then t 0 B H s d -B H s = -∞ a.s.,
meaning that there is no possible change-of-variable formula for f (x) = x 2 .

Thus, H = 1 2 appears to be a critical value for the change-of-variable formula involving the forward integral (1.2). This is because it is precisely the value from which the sign of 2H -1 changes in (1.4). The chain rule being (1.1) in the critical case H = 1 2 , one has a complete picture for the forward integral (1.2).

To go one step further, one may wonder what kind of change-of-variable formula one would obtain after replacing the definition (1.2) by its symmetric counterpart, namely

t 0 X s d • Y s := lim n→∞ ⌊2 n t⌋-1 k=0 1 2 X k2 -n + X (k+1)2 -n (Y (k+1)2 -n -Y k2 -n ) (1.5)
(provided the limit exists in some sense). As it turns out, it is arguably a much more difficult problem, which has been solved only recently. In this context, the crucial quantity is now the cubic variation. And this latter is known to satisfy, for any H < 1 2 ,

2 n(3H-1 2 ) ⌊2 n t⌋-1 k=0 (B H (k+1)2 -n -B H k2 -n ) 3 law → N(0, σ 2 H ) as n → ∞. (1.6)
With a lot of efforts, one can prove (see [START_REF] Cheridito | Stochastic integral of divergence type with respect to fractional Brownian motion with Hurst parameter H ∈ (0, 1 2 )[END_REF][START_REF] Gradinaru | m-order integrals and generalized Itô's formula: the case of a fractional Brownian motion with any Hurst index[END_REF] when H = 1 6 and [START_REF] Nourdin | The weak Stratonovich integral with respect to fractional Brownian motion with Hurst parameter 1/6[END_REF] when H = 1 6 ) the following three facts, which hold for any smooth enough real function f : R → R:

1. If H > 1 6 then t 0 f ′ (B H s )d • B H
s exists as a limit in probability and one has

f (B H t ) = f (0) + t 0 f ′ (B H s )d • B H s , t 0. (1.7) 2. If H = 1 6 then t 0 f ′ (B 1 6 s )d • B 1 6
s exists as a stable limit in law and one has, with W a standard Brownian motion independent of B 1 6 and κ 3 ≃ 2.322,

f (B 1 6 t ) = f (0) + t 0 f ′ (B 1 6 s )d • B 1 6 s - κ 3 12 t 0 f ′′′ (B 1 6
)dW s , t 0.

(1.8)

3. If H < 1 6 then t 0 (B H s ) 2 d • B H s does not exist in law.
(1.9)

Thus, as we see, the critical value for the symmetric integral is now H = 1 6 ; it is exactly the value of H from which the sign of 3H -1 2 changes in (1.6).

In [START_REF] Burdzy | Some path properties of iterated Brownian motion[END_REF][START_REF] Burdzy | Variation of iterated Brownian motion[END_REF] (see also [START_REF] Burdzy | Brownian motion in a Brownian crack[END_REF]), Burdzy has introduced the so-called iterated Brownian motion. This process, which can be regarded as the realization of a Brownian motion on a random fractal, is defined as

Z t = X(Y t ), t 0,
where X is a two-sided Brownian motion and Y is a standard (one-sided) Brownian motion independent of X. Note that Z is self-similar of order 1 4 and has stationary increments; hence, in some sense, Z is close to the fractional Brownian motion B 1 4 of index H = 1 4 . As is the case for B 1 4 , Z is neither a Dirichlet process nor a semimartingale or a Markov process in its own filtration. A crucial question is therefore how to define a stochastic calculus with respect to it. This issue has been tackled by Khoshnevisan and Lewis in [START_REF] Khoshnevisan | Stochastic calculus for Brownian motion on a Brownian fracture[END_REF][START_REF] Khoshnevisan | Iterated Brownian motion and its intrinsic skeletal structure[END_REF], where the authors develop a Stratonovich-type stochastic calculus with respect to Z, by extensively using techniques based on the properties of some special arrays of Brownian stopping times, as well as on excursion-theoretic arguments. See also the paper [START_REF] Nourdin | Weighted power variations of iterated Brownian motion[END_REF] which may be seen as a follow-up of [START_REF] Khoshnevisan | Stochastic calculus for Brownian motion on a Brownian fracture[END_REF]. The formula obtained in [START_REF] Khoshnevisan | Stochastic calculus for Brownian motion on a Brownian fracture[END_REF][START_REF] Khoshnevisan | Iterated Brownian motion and its intrinsic skeletal structure[END_REF] reads, unsurprisingly (due to (1.7) and the similarities between Z and B 1 4 ) and losely speaking, as follows:

f (Z t ) = f (0) + t 0 f (Z s )d • Z s , t 0.
(1.10)

The change-of-variable formula (1.10) is of the same kind than (1.7). In view of what has been done so far for the fractional Brownian motion B H , aiming to provide an answer to the following problem is somehow natural: can we also reinterpret (1.10) in a dynamical way, in the spirit of (1.7), (1.8) and (1.9)? To this end, we first need to introduce a family of processes that contains the iterated Brownian motion Z as a particular element. The family consisting in the so-called fractional Brownian motions in Brownian time, studied in [START_REF] Zeineddine | Fluctuations of the power variation of fractional Brownian motion in Brownian time[END_REF] by the second-named author, does the job. More specifically, it is the family {Z H } H∈(0,1) defined as follows:

Z H t = X H (Y t ), t 0,
where X H is a two-sided fractional Brownian motion of index H and Y is a standard (onesided) Brownian motion independent of X. Roughly speaking, in the present paper we are going to show the following three assertions (see Theorem 2.1 for a precise statement): for any smooth real function f : R → R,

1. If H > 1 6 then f (Z H t ) = f (0) + t 0 f ′ (Z H s )d • Z H s , t 0.
2. If H = 1 6 then, with W a standard Brownian motion independent of the pair (X 1 6 , Y ) and κ 3 ≃ 2.322, f (Z

1 6 t ) = f (0) + t 0 f ′ (Z 1 6 s )d • Z 1 6 s - κ 3 12 Yt 0 f ′′′ (X 1 6 s )dW s , t 0. (1.11) 3. If H < 1 6 , then t 0 (Z H s ) 2 d • Z H s does not exist.
The formula (1.11) is related to a recent line of research in which, by means of Malliavin calculus, one aims to exhibit change-of-variable formulas in law with a correction term which is an Itô integral with respect to martingale independent of the underlying Gaussian processes. Papers dealing with this problem and which are prior to our work include [START_REF] Burdzy | A change of variable formula with Itô correction term[END_REF][START_REF] Harnett | Central limit theorem for a Stratonovich integral with Malliavin calculus[END_REF][START_REF] Harnett | Weak convergence of the Stratonovich integral with respect to a class of Gaussian processes[END_REF][START_REF] Harnett | On Simpson's rule and fractional Brownian motion with H = 1/10[END_REF][START_REF] Nourdin | A change of variable formula for the 2D fractional Brownian motion of Hurst index bigger or equal to 1/4[END_REF][START_REF] Nourdin | Asymptotic behavior of weighted quadratic variations of fractional Brownian motion: the critical case H = 1/4[END_REF][START_REF] Nourdin | The weak Stratonovich integral with respect to fractional Brownian motion with Hurst parameter 1/6[END_REF]; however, it is worthwhile noting that all these mentioned references only deal with Gaussian processes, not with iterated processes (which are arguably more difficult to handle).

A brief outline of the paper is as follows. In Section 2, we introduce the framework in which our study takes place and we provide an exact statement of our result, namely Theorem 2.1. Finally, Section 3 contains the proof of Theorem 2.1, which is divided into several steps.

Framework and exact statement of our results

For simplicity, throughout the paper we remove the superscript H, that is, we write Z (resp. X) instead of Z H (resp. X H ).

Let Z be a fractional Brownian motion in Brownian time of Hurst parameter H ∈ (0, 1), defined as

Z t = X(Y t ), t 0, (2.12)
where X is a two-sided fractional Brownian motion of parameter H and Y is a standard (one-sided) Brownian motion independent of X.

The paths of Z being very irregular (precisely: Hölder continuous of order α if and only if α is strictly less than H/2), we will not be able to define a stochastic integral with respect to it as the limit of Riemann sums with respect to a deterministic partition of the time axis. However, a winning idea borrowed from Khoshnevisan and Lewis [START_REF] Khoshnevisan | Stochastic calculus for Brownian motion on a Brownian fracture[END_REF][START_REF] Khoshnevisan | Iterated Brownian motion and its intrinsic skeletal structure[END_REF] is to approach deterministic partitions by means of random partitions defined in terms of hitting times of the underlying Brownian motion Y . As such, one can bypass the random "time-deformation" forced by (2.12), and perform asymptotic procedures by separating the roles of X and Y in the overall definition of Z.

Following Khoshnevisan and Lewis [START_REF] Khoshnevisan | Stochastic calculus for Brownian motion on a Brownian fracture[END_REF][START_REF] Khoshnevisan | Iterated Brownian motion and its intrinsic skeletal structure[END_REF], we start by introducing the so-called intrinsic skeletal structure of Z H . This structure is defined through a sequence of collections of stopping times (with respect to the natural filtration of Y ), noted

T n = {T k,n : k 0}, n 1, (2.13)
which are in turn expressed in terms of the subsequent hitting times of a dyadic grid cast on the real axis. More precisely, let D n = {j2 -n/2 : j ∈ Z}, n 1, be the dyadic partition (of R) of order n/2. For every n 1, the stopping times T k,n , appearing in (2.13), are given by the following recursive definition: T 0,n = 0, and

T k,n = inf s > T k-1,n : Y (s) ∈ D n \ {Y (T k-1,n )} , k 1.
Note that the definition of T k,n , and therefore of T n , only involves the one-sided Brownian motion Y , and that, for every n 1, the discrete stochastic process (2.14)

Y n = {Y (T k,n ) : k 0}
Based on this fact, one may introduce the counterpart of (1.5) based on T n , namely,

V n (f, t) = ⌊2 n t⌋-1 k=0 1 2 f (Z T k,n ) + f (Z T k+1,n ) (Z T k+1,n -Z T k,n ). (2.15)
Let C ∞ b denote the class of those functions f : R → R that are C ∞ and bounded together with their derivatives. We then have the following result.

Theorem 2.1 Let f ∈ C ∞ b and t > 0. 1. If H > 1 6 then f (Z t ) -f (0) = t 0 f ′ (Z s )d • Z s , (2.16)
where

t 0 f ′ (Z s )d • Z s is the limit in probability of V n (f ′ , t) defined in (2.15) as n → ∞. 2. If H = 1 6 then, with κ 3 ≃ 2.322, f (Z t ) -f (0) + κ 3 12 Yt 0 f ′′′ (X s )dW s (law) = t 0 f ′ (Z s )d • Z s .
(2.17) (2.15) as n → ∞ (its existence is part of the conclusion), W is a two-sided Brownian motion independent of the pair (X, Y ) defining Z, and the integral with respect to W is understood in the Wiener-Itô sense.

Here, t 0 f ′ (Z s )d • Z s denotes the limit in law of V n (f ′ , t) defined in
3. If H < 1 6 then

V n (x → x 2 , t) does not converge, even stably in law.

(

2.18)

This means that there is no way to get a change-of-variable formula for f (x) = x 3 .

3 Proof of Theorem 2.1

Elements of Malliavin calculus

In this section, we gather some elements of Malliavin calculus we shall need thoughout the proof of Theorem 2.1. The reader already familiar with this topic may skip this section. We continue to denote by X = (X t ) t∈R a two-sided fractional Brownian motion with Hurst parameter H ∈ (0, 1). That is, X is a zero mean Gaussian process, defined on a complete probability space (Ω, A , P ), with the covariance function

C H (t, s) = E(X t X s ) = 1 2 (|s| 2H + |t| 2H -|t -s| 2H ), s, t ∈ R.
We suppose that A is the σ-field generated by X. For all n ∈ N * , we let E n be the set of step functions on [-n, n], and

E := ∪ n E n . Set ξ t = 1 [0,t] (resp. 1 [t,0] ) if t 0 (resp. t < 0).
Let H be the Hilbert space defined as the closure of E with respect to the inner product

ξ t , ξ s H = C H (t, s), s, t ∈ R.
The mapping ξ t → X t can be extended to an isometry between H and the Gaussian space H 1 associated with X. We will denote this isometry by ϕ → X(ϕ).

Let F be the set of all smooth cylindrical random variables, i.e. of the form

F = φ(X t 1 , ..., X t l ),
where l ∈ N * , φ : R l → R is C ∞ b and t 1 < ... < t l are some real numbers. The derivative of F with respect to X is the element of L 2 (Ω, H ) defined by

D s F = l i=1 ∂φ ∂x i (X t 1 , ..., X t l )ξ t i (s), s ∈ R.
In particular D s X t = ξ t (s). For any integer k 1, we denote by D k,2 the closure of the set of smooth random variables with respect to the norm

F 2 k,2 = E(F 2 ) + k j=1 E[ D j F 2 H ⊗j ].
The Malliavin derivative D satisfies the chain rule.

If ϕ : R n → R is C 1 b and if F 1 , . . . , F n are in D 1,2 , then ϕ(F 1 , ..., F n ) ∈ D 1,2 and we have Dϕ(F 1 , ..., F n ) = n i=1 ∂ϕ ∂x i (F 1 , ..., F n )DF i .
We have the following Leibniz formula, whose proof is straightforward by induction on q. Let ϕ, ψ ∈ C q b (q 1), and fix 0 u < v and 0 s < t. Then ϕ(X t -X s )ψ(X v -X u ) ∈ D q,2 and

D q ϕ(X t -X s )ψ(X v -X u ) = 2q a=0 2q a ϕ (a) (X t -X s )ψ (2q-a) (X v -X u )1 ⊗a [s,t] ⊗1 ⊗(2q-a) [u,v] , (3.19)
where ⊗ stands for the symmetric tensor product. A similar statement holds fo u < v 0 and s < t 0.

If a random element u ∈ L 2 (Ω, H ) belongs to the domain of the divergence operator, that is, if it satisfies

|E DF, u H | c u E(F 2 ) for any F ∈ F , then I(u) is defined by the duality relationship E F I(u) = E DF, u H , for every F ∈ D 1,2 .
For every n 1, let H n be the nth Wiener chaos of X, that is, the closed linear subspace of L 2 (Ω, A , P ) generated by the random variables {H n (B(h)), h ∈ H , h H = 1}, where H n is the nth Hermite polynomial. The mapping I n (h ⊗n ) = H n (B(h)) provides a linear isometry between the symmetric tensor product H ⊙n and H n . For H = 1 2 , I n coincides with the multiple Wiener-Itô integral of order n. The following duality formula holds

E F I n (h) = E D n F, h H ⊗n , (3.20)
for any element h ∈ H ⊙n and any random variable F ∈ D n,2 . Let {e k , k 1} be a complete orthonormal system in H . Given f ∈ H ⊙n and g ∈ H ⊙m , for every r = 0, ..., n ∧ m, the contraction of f and g of order r is the element of H ⊗(n+m-2r) defined by

f ⊗ r g = ∞ k 1 ,...,kr=1
f, e k 1 ⊗ ... ⊗ e kr H ⊗r ⊗ g, e k 1 ⊗ ... ⊗ e kr H ⊗r .

Note that f ⊗ r g is not necessarily symmetric: we denote its symmetrization by f ⊗r g ∈ H ⊙(n+m-2r) . Finally, we recall the following product formula: if f ∈ H ⊙n and g ∈ H ⊙m then

I n (f )I m (g) =

Notation and reduction of the problem

Throughout all the proof, we shall use the following notation. For all k, n ∈ N we write

ξ k2 -n/2 = 1 [0,k2 -n/2 ] , ξ - k2 -n/2 = 1 [-k2 -n/2 ,0] , δ k2 -n/2 = 1 [(k-1)2 -n/2 ,k2 -n/2 ] , δ - k2 -n/2 = 1 [-k2 -n/2 ,(-k+1)2 -n/2 ]
. Also, •, • ( • , respectively) will always stand for inner product (the norm, respectively) in an appropriate tensor product H ⊗s .

In the sequel, we only consider the case H < 1 2 . The proof of (2.16) in the case H > 1 2 is easier and left to the reader, whereas the proof when H = 1 2 was already done in [START_REF] Khoshnevisan | Stochastic calculus for Brownian motion on a Brownian fracture[END_REF][START_REF] Khoshnevisan | Iterated Brownian motion and its intrinsic skeletal structure[END_REF] by Khoshnevisan and Lewis.

That said, we now divide the proof of Theorem 2.1 in several steps.

Step 1: A key algebraic lemma

For each integer n 1, k ∈ Z and real number t 0, let U j,n (t) (resp. D j,n (t)) denote the number of upcrossings (resp. downcrossings) of the interval [j2 -n/2 , (j + 1)2 -n/2 ] within the first ⌊2 n t⌋ steps of the random walk 2 and Y (T k+1,n ) = (j + 1)2 -n/2 ; D j,n (t) = ♯ k = 0, . . . , ⌊2 n t⌋ -1 :

{Y (T k,n )} k 1 , that is, U j,n (t) = ♯ k = 0, . . . , ⌊2 n t⌋ -1 : Y (T k,n ) = j2 -n/
Y (T k,n ) = (j + 1)2 -n/2 and Y (T k+1,n ) = j2 -n/2 .
While easy, the following lemma taken from [START_REF] Khoshnevisan | Stochastic calculus for Brownian motion on a Brownian fracture[END_REF]Lemma 2.4] is going to be the key when studying the asymptotic behavior of the weighted power variation V (2r-1) n (f, t) of odd order 2r -1, defined as:

V (2r-1) n (f, t) = ⌊2 n t⌋-1 k=0 1 2 f (Z T k,n ) + f (Z T k+1,n ) (Z T k+1,n -Z T k,n ) 2r-1 , t 0.
Its main feature is to separate X from Y , thus providing a representation of V (2r-1) n (f, t) which is amenable to analysis.

Lemma 3.1 Fix f ∈ C ∞ b , t 0 and r ∈ N * . Then V (2r-1) n (f, t) = j∈Z 1 2 f (X j2 -n 2 ) + f (X (j+1)2 -n 2 ) X (j+1)2 -n 2 -X j2 -n 2 2r-1 U j,n (t) -D j,n (t) .
Observe that V

(1)

n (f, t) = V n (f, t), see (2.15).

Step 2: Transforming the weighted power variations of odd order

By [10, Lemma 2.5], one has

U j,n (t) -D j,n (t) =    1 {0 j<j * (n,t)} if j * (n, t) > 0 0 if j * = 0 -1 {j * (n,t) j<0} if j * (n, t) < 0 , where j * (n, t) = 2 n/2 Y T ⌊2 n t⌋,n . As a consequence, V (2r-1) n (f, t) is equal to      2 -nH(r-1 2 ) j * (n,t)-1 j=0
1 2 f (X + j2 -n/2 ) + f (X + (j+1)2 -n/2 ) X n,+ j+1 -X n,+ j 2r-1 if j * (n, t) > 0 0 if j * = 0 2 -nH(r-1 2 ) |j * (n,t)|-1 j=0 1 2 f (X - j2 -n/2 ) + f (X - (j+1)2 -n/2 ) X n,- j+1 -X n,- j 2r-1 if j * (n, t) < 0
, where X + t := X t for t 0, X - -t := X t for t < 0, X n,+ (t) := 2 nH/2 X + 2 -n/2 t for t 0 and X n,-(-t) := 2 nH/2 X - 2 -n/2 (-t) for t < 0. Let us now introduce the following sequence of processes

W (2r-1) ±,n
, in which H p stands for the pth Hermite polynomial:

W (2r-1) ±,n (f, t) = ⌊2 n/2 t⌋-1 j=0 1 2 f (X ± j2 -n 2 ) + f (X ± (j+1)2 -n 2 ) H 2r-1 (X n,± j+1 -X n,± j ), t 0
W (2r-1) n (f, t) = W (2r-1) +,n (f, t) if t 0 W (2r-1) -,n (f, -t) if t < 0 .
We then have, using the decomposition x 2r-1 = r l=1 a r,l H 2l-1 (x) (with a r,r = 1, which is the only explicit value of a l,r we will need in the sequel),

V (2r-1) n (f, t) = 2 -nH(r-1 2 ) r l=1 a r,l W (2l-1) n (f, Y T ⌊2 n t⌋,n ).
(3.22)

Step 3: Known results for fractional Brownian motion

We recall the following result taken * from [START_REF] Nourdin | Central and non-central limit theorems for weighted power variations of the fractional Brownian motion[END_REF] . If m 2 and H ∈

1 4m-2 , 1 2 then, for any f ∈ C ∞ b and as n → ∞, X t , 2 -n/4 W (2m-1) ±,n (f, t) t 0 fdd -→ X t , κ 2m-1 t 0 f (X ± s )dW ± s t 0 , (3.23
) * More precisely: a careful inspection would show that there is no additional difficulty to prove (3.23) by following the same route than the one used to show [START_REF] Nourdin | Central and non-central limit theorems for weighted power variations of the fractional Brownian motion[END_REF]Theorem 1,(1.15)]. The only difference is that the definition of W (r) ±,n is of symmetric type, whereas all the quantities of interest studied in [START_REF] Nourdin | Central and non-central limit theorems for weighted power variations of the fractional Brownian motion[END_REF] are of forward type.

where W + t = W t if t > 0 and W - t = W -t if t < 0, with W a two-sided Brownian motion independent of X, and where t 0 f (X ± s )dW ± s must be understood in the Wiener-Itô sense. Note that in the boundary case m = 2 and H = 1 6 , (3.23) continues to hold, as was shown in [START_REF] Nourdin | The weak Stratonovich integral with respect to fractional Brownian motion with Hurst parameter 1/6[END_REF]Theorem 3.1].

In the case m = 1, it was shown in [START_REF] Nourdin | Central and non-central limit theorems for weighted power variations of the fractional Brownian motion[END_REF]Theorem 4] (case H > 1 6 ) and [16, Theorem 2.13] (case H = 1 6 ) that, for any fixed t > 0, the sequence W

±,n (f, t) converges in probability (when H > 1 6 ) or only in law (when H = 1 6 ) to a non degenerate limit as n → ∞.

3.6

Step 4: Moment bounds for

W (2r-1) n (f, •)
Fix an integer r 1 as well as a function f ∈ C ∞ b . We claim the existence of c > 0 such that, for all real numbers s < t and all n ∈ N,

E W (2r-1) n (f, t) -W (2r-1) n (f, s) 2 c max(|s| 2H , |t| 2H ) |t -s|2 n/2 + 1 . (3.24)
In order to prove (3.24), we will need the following lemma.

Lemma 3.2 If s, t, u > 0 or if s, t, u < 0 then |E X u (X t -X s ) | |t -s| 2H . (3.25)
Proof. When s, t, u > 0 we have

E X u (X t -X s ) = 1 2 t 2H -s 2H + 1 2 |s -u| 2H -|t -u| 2H .
Since |b 2H -a 2H | |b -a| 2H for any a, b ∈ R + , we immediately deduce (3.25). The proof when s, t, u < 0 is similar.

We are now ready to show (3.24). We distinguish two cases according to the signs of s, t ∈ R (and reducing the problem by symmetry):

(1) if 0 s < t (the case s < t 0 being similar), then

E[(W (2r-1) n (f, t) -W (2r-1) n (f, s)) 2 ] = E[(W (2r-1) +,n (f, t) -W (2r-1) +,n (f, s)) 2 ] = 1 4 ⌊2 n/2 t⌋-1 j,j ′ =⌊2 n/2 s⌋ E f (X + j2 -n 2 ) + f (X + (j+1)2 -n 2 ) × f (X + j ′ 2 -n 2 ) + f (X + (j ′ +1)2 -n 2 ) H 2r-1 (X n,+ j+1 -X n,+ j )H 2r-1 (X n,+ j ′ +1 -X n,+ j ′ ) = 1 4 2 nH(2r-1) ⌊2 n/2 t⌋-1 j,j ′ =⌊2 n/2 s⌋ E Θ n j f (X + )Θ n j ′ f (X + )I 2r-1 (δ ⊗(2r-1) (j+1)2 -n/2 )I 2r-1 (δ ⊗(2r-1) (j ′ +1)2 -n/2 ) ,
with obvious notation. Relying to the product formula (3.21), we deduce that this latter quantity is less than or equal to

1 4 2 nH(2r-1) ⌊2 n/2 t⌋-1 j,j ′ =⌊2 n/2 s⌋ 2r-1 l=0 l! 2r -1 l 2 δ (j+1)2 -n/2 ; δ (j ′ +1)2 -n/2 l × E Θ n j f (X + )Θ n j ′ f (X + )I 4r-2-2l (δ ⊗(2r-1-l) (j+1)2 -n/2 ⊗δ ⊗(2r-1-l) (j ′ +1)2 -n/2 ) =: 1 4 2r-1 l=0 l! 2r -1 l 2 Q (+,
d (+,r,l) n (j, j ′ ) := E Θ n j f (X + )Θ n j ′ f (X + )I 4r-2-2l (δ ⊗(2r-1-l) (j+1)2 -n/2 ⊗δ ⊗(2r-1-l) (j ′ +1)2 -n/2 ) = E D 4r-2-2l (Θ n j f (X + )Θ n j ′ f (X + )) ; δ ⊗(2r-1-l) (j+1)2 -n/2 ⊗δ ⊗(2r-1-l) (j ′ +1)2 -n/2 = 4r-2-2l a=0 4r -2 -2l a E f (a) (X + j2 -n/2 )ξ ⊗a j2 -n/2 + f (a) (X + (j+1)2 -n/2 )ξ ⊗a (j+1)2 -n/2 ⊗ f (4r-2-2l-a) (X + j ′ 2 -n/2 )ξ ⊗(4r-2-2l-a) j ′ 2 -n/2 + f (4r-2-2l-a) (X + (j ′ +1)2 -n/2 )ξ ⊗(4r-2-2l-a) (j ′ +1)2 -n/2 ; δ ⊗(2r-1-l) (j+1)2 -n/2 ⊗δ ⊗(2r-1-l) (j ′ +1)2 -n/2
.

Let now c denote a generic constant that may differ from one line to another and recall that f ∈ C ∞ b . We then have the following estimates.

• Case l = 2r -1

Q (+,r,2r-1) n (s, t) c 2 nH(2r-1) ⌊2 n/2 t⌋-1 j,j ′ =⌊2 n/2 s⌋ δ (j+1)2 -n/2 ; δ (j ′ +1)2 -n/2 2r-1 = c ⌊2 n/2 t⌋-1 j,j ′ =⌊2 n/2 s⌋ 1 2 (|j -j ′ + 1| 2H + |j -j ′ -1| 2H -2|j -j ′ | 2H ) 2r-1 = c ⌊2 n/2 t⌋-1 j=⌊2 n/2 s⌋ j-⌊2 n/2 s⌋ q=j-⌊2 n/2 t⌋+1 ρ(q) 2r-1 , with ρ(q) := 1 2 (|q + 1| 2H + |q -1| 2H -2|q| 2H ). By a Fubini argument, it comes Q (+,r,2r-1) n (s, t) c ⌊2 n/2 t⌋-⌊2 n/2 s⌋-1 q=⌊2 n/2 s⌋-⌊2 n/2 t⌋+1 |ρ(q)| 2r-1 (q + ⌊2 n/2 t⌋) ∧ ⌊2 n/2 t⌋ -(q + ⌊2 n/2 s⌋) ∨ ⌊2 n/2 s⌋ c ⌊2 n/2 t⌋-⌊2 n/2 s⌋-1 q=⌊2 n/2 s⌋-⌊2 n/2 t⌋+1 |ρ(q)| 2r-1 ⌊2 n/2 t⌋ -⌊2 n/2 s⌋ c q∈Z |ρ(q)| 2r-1 ⌊2 n/2 t⌋ -⌊2 n/2 s⌋ = c ⌊2 n/2 t⌋ -⌊2 n/2 s⌋ c ⌊2 n/2 t⌋ -2 n/2 t + 2 n/2 t -s + ⌊2 n/2 s⌋ -2 n/2 s c(1 + 2 n/2 |t -s|). (3.27) Note that q∈Z |ρ(q)| 2r-1 < ∞ since H < 1 2 1 -1 4r-2 .
• Preparation to the cases where 0 l 2r -2

In order to handle the terms Q (+,r,l) n (s, t) whenever 0 l 2r -2, we will make use of the following decomposition:

|d (+,r,l) n (j, j ′ )| 1 u,v=0 Ω (u,v,r,l) n (j, j ′ ), (3.28)
where

Ω (u,v,r,l) n (j, j ′ ) = 4r-2-2l a=0 4r -2 -2l a E[f (a) (X + (j+u)2 -n/2 )f (4r-2-2l-a) (X + (j ′ +v)2 -n/2 )] × ξ ⊗a (j+u)2 -n/2 ⊗ξ ⊗(4r-2-2l-a) (j ′ +v)2 -n/2 ; δ ⊗(2r-1-l) (j+1)2 -n/2 ⊗δ ⊗(2r-1-l) (j ′ +1)2 -n/2 . • Case 1 l 2r -2 (only when r 2)
Since f belongs to C ∞ b and since, by (3.25), we have | ξ t ; δ (j+1)2 -n/2 | 2 -nH for all t 0 and all j ∈ N, we deduce that

|d (+,r,l) n (j, j ′ )| c 2 -nH(4r-2-2l) .
As a consequence, and relying to the same arguments that have been used previously in the case l = 2r -1, we get

Q (+,r,l) n (s, t) c 2 -nH(4r-2-2l) 2 nH(2r-1) ⌊2 n/2 t⌋-1 j,j ′ =⌊2 n/2 s⌋ δ (j+1)2 -n/2 ; δ (j ′ +1)2 -n/2 l c 2 -nH(4r-2-2l) 2 nH(2r-1) 2 -nHl q∈Z |ρ(q)| l (1 + 2 n/2 |t -s|) = c 2 -nH(2r-1-l) (1 + 2 n/2 |t -s|) c (1 + 2 n/2 |t -s|). (3.29)
• Case l = 0

Relying to the decomposition (3.28), we get

Q (+,r,0) n (s, t) 2 nH(2r-1) ⌊2 n/2 t⌋-1 j,j ′ =⌊2 n/2 s⌋ 1 u,v=0
Ω (u,v,r,0) n (j, j ′ ).

(3.30)

We will study only the term corresponding to Ω (0,1,r,0) n (j, j ′ ) in (3.30), which is representative of the difficulty. It is given by

2 nH(2r-1) ⌊2 n/2 t⌋-1 j,j ′ =⌊2 n/2 s⌋ 4r-2 a=0 4r -2 a E[f (a) (X + j2 -n/2 )f (4r-2-a) (X + (j ′ +1)2 -n/2 )] × ξ ⊗a j2 -n/2 ⊗ξ ⊗(4r-2-a) (j ′ +1)2 -n/2 ; δ ⊗(2r-1) (j+1)2 -n/2 ⊗δ ⊗(2r-1) (j ′ +1)2 -n/2 c 2 nH(2r-1) ⌊2 n/2 t⌋-1 j,j ′ =⌊2 n/2 s⌋ 4r-2 a=0 ξ ⊗a j2 -n/2 ⊗ξ ⊗(4r-2-a) (j ′ +1)2 -n/2 ; δ ⊗(2r-1) (j+1)2 -n/2 ⊗δ ⊗(2r-1) (j ′ +1)2 -n/2 .
We define

E (a,r) n (j, j ′ ) := ξ ⊗a j2 -n/2 ⊗ξ ⊗(4r-2-a) (j ′ +1)2 -n/2 ; δ ⊗(2r-1) (j+1)2 -n/2 ⊗δ ⊗(2r-1) (j ′ +1)2 -n/2 . By (3.25), recall that | ξ t ; δ (j+1)2 -n/2 |
2 -nH for all t 0 and all j ∈ N. We thus get, with ca some combinatorial constants,

E (a,r) n (j, j ′ ) ca 2 -nH(4r-3) | ξ j2 -n/2 ; δ (j+1)2 -n/2 | + | ξ j2 -n/2 ; δ (j ′ +1)2 -n/2 | +| ξ (j ′ +1)2 -n/2 ; δ (j+1)2 -n/2 | + | ξ (j ′ +1)2 -n/2 ; δ (j ′ +1)2 -n/2 | .
For instance, we can write

⌊2 n/2 t⌋-1 j,j ′ =⌊2 n/2 s⌋ | ξ (j ′ +1)2 -n/2 ; δ (j+1)2 -n/2 | = 2 -nH-1 ⌊2 n/2 t⌋-1 j,j ′ =⌊2 n/2 s⌋ (j + 1) 2H -j 2H + |j ′ -j + 1| 2H -|j ′ -j| 2H 2 -nH-1
⌊2 n/2 t⌋-1 j,j ′ =⌊2 n/2 s⌋ (j + 1) 2H -j 2H +2 -nH-1

⌊2 n/2 s⌋ j j ′ ⌊2 n/2 t⌋-1

(j ′ -j + 1) 2H -(j ′ -j) 2H +2 -nH-1
⌊2 n/2 s⌋ j ′ <j ⌊2 n/2 t⌋-1

(j -j ′ ) 2H -(j -j ′ -1) 2H with W = (W t ) t∈R a two-sided Brownian motion independent of the pair (X, Y ). Indeed, using the decomposition (3.22), the conclusion of Step 4 (to pass from Y T ⌊2 n t⌋,n to Y t ) and the convergence (2.14), we deduce that the limit of V (2r-1) n (f, t) is the same as that of

2 -nH(r-1 2 ) r l=1 a r,l W (2l-1) n (f, Y t ).
Thus, the proofs of (3.32), (3.33) and (3.34) then follow directly from the results recalled in Step 3, as well as the fact that X and Y are independent.

3.8

Step 6: Proving (2.16) and (2.17)

We assume H ∈ [ 1 6 , 1 2 ). We will make use of the following Taylor's type formula.

f (Z T ⌊2 n t⌋,n ) -f (0) = ⌊2 n t⌋-1 k=0 f (Z T k+1,n ) -f (Z T k,n ) = ⌊2 n t⌋-1 k=0 1 2 f (Z T k,n ) + f (Z T k+1,n ) (Z T k+1,n -Z T k,n ) (3.35) - 1 12 V (3) n (f, t) + 7 r=3 2c r V (2r-1) n (f, t) + ⌊2 n t⌋-1 k=0 O((Z T k+1,n -Z T k,n ) 14 ).
As far as the big O in (3.35) is concerned, we have, with G ∼ N(0, 1),

E ⌊2 n t⌋-1 k=0 O((Z T k+1,n -Z T k,n ) 14 ) C f ⌊2 n t⌋-1 k=0 E (Z T k+1,n -Z T k,n ) 14 = C f ⌊2 n t⌋-1 k=0 2 -7nH E[G 14 ] C f E[G 14 ]t 2 n(

  defines a simple random walk over D n . As shown in [10, Lemma 2.2], as n tends to infinity the collection {T k,n : 1 k 2 n t} approximates the common dyadic partition {k2 -n : 1 k 2 n t} of order n of the time interval [0, t]. More precisely, sup 0 s t T ⌊2 n s⌋,n -s → 0 almost surely and in L 2 (Ω).

3 + 7 r=3c

 37 Fix f ∈ C ∞ b .For any a, b ∈ R and for some constants c r whose explicit values are immaterial here,f (b) -f (a) = 1 2 f ′ (a) + f ′ (b) (b -a) -1 24 f ′′′ (a) + f ′′′ (b) (b -a) r f (2r-1) (a) + f (2r-1) (b) (b -a) 2r-1 + O(|b -a| 14 ),where |O(|b -a|14 )| C f |b -a| 14 , C f being a constant depending only on f . One can thus write

  -nH ⌊2 n/2 t⌋ -⌊2 n/2 s⌋ ⌊2 n/2 t⌋ 2H 3t 2H 2 1 + 2 n/2 |t -s| . whereas, if H = 1 6 then, as n → ∞, X t , Y t , V (3) n (f, t) t 0 fdd → X t , Y t , κ 3

	3 2	0 2 Yt	f (X s )dW s	t 0	,	(3.34)

  1-7H) → n→∞ 0 since H

	1 6 .	(3.36)

; IN is supported in part by the (french) ANR grant 'Malliavin, Stein and

Similarly,

⌊2 n/2 t⌋-1 j,j ′ =⌊2 n/2 s⌋ ξ j2 -n/2 ; δ (j+1)2 -n/2 | 3t 2H 2 1 + 2 n/2 |t -s| ;

⌊2 n/2 t⌋-1 j,j ′ =⌊2 n/2 s⌋ ξ j2 -n/2 ; δ (j ′ +1)2 -n/2 | 3t 2H 2 1 + 2 n/2 |t -s| ;

As a consequence, we deduce (

By (1) with s = 0, one can write

.

Similarly

We deduce that

That is, (3.24) also holds true in this case.

Step 5: Limits of the weighted power variations of odd order

Fix f ∈ C ∞ b and t 0. We claim that, if H ∈ 1 6 , 1 2 and r 3 then, as n → ∞,

On the other hand, by continuity of f • Z and due to (2.14), one has, almost surely and as n → ∞,

(3.37)

Finally, when H > 1 6 the desired conclusion (2.16) follows from (3.36), (3.37), (3.32) and (3.33) plugged into (3.35). The proof of (2.17) when H = 1 6 is similar, the only difference being that one has (3.34) instead of (3.33), thus leading to the bracket term

3.9 Step 7: Proving (2.18) 3 , one can write, with 1 denoting the function constantly equal to 1,

As a result, and thank to (2.14), one deduces that if V n (x → x 2 , t) converges stably in law, then V

n (1, t) must converge as well. But it is shown in [17, Corollary 1.2] that 2 -n(1-6H)/4 V (3) n (1, t) converges in law to a non degenerate limit. This being clearly in contradiction with the convergence of V n (1, t), we deduce that (2.18) holds.


Yuedong Xu
email: yuedong.xu@gmail.com

Salah Eddine Elayoubi
email: salaheddine.elayoubi@orange.com

Eitan Altman
email: eitan.altman@inria.fr

Rachid El-Azouzi
email: rachid.elazouzi@univ-avignon.fr

Impact of Flow-level Dynamics on QoE of Video Streaming in Wireless Networks

The Quality of Experience (QoE) of streaming service is often degraded by frequent playback interruptions. To mitigate the interruptions, the media player prefetches streaming contents before starting playback, at a cost of delay. We study the QoE of streaming from the perspective of flow dynamics. First, a framework is developed for QoE when streaming users join the network randomly and leave after downloading completion. We compute the distribution of prefetching delay using partial differential equations (PDEs), and the probability generating function of playout buffer starvations using ordinary differential equations (ODEs). Second, we extend our framework to characterize the throughput variation caused by opportunistic scheduling at the base station in the presence of fast fading. Our study reveals that the flow dynamics is the fundamental reason of playback starvation. The QoE of streaming service is dominated by the average throughput of opportunistic scheduling, while the variance of throughput has very limited impact on starvation behavior.

I. INTRODUCTION

Streaming services are witnessing a rapid growth in mobile networks. According to Allot Communications [START_REF]Allot MobileTrends[END_REF], HTTP streaming service made up 37 percent of mobile broadband traffic during the second half of 2010. This presents new challenges for operators that are used to classify services into real-time (voice-like) and elastic (data-like) services. Indeed, classical QoS metrics in mobile networks are blocking rates for real-time traffic and average user throughput for elastic one, and operators dimension their networks for satisfying targets on those metrics [START_REF] Borst | Integration of Streaming and Elastic Traffic in Wireless Networks[END_REF]. However, the particular nature of streaming applications, halfway between real-time and elastic services, is raising the following difficult questions in wireless environments. First, which QoS metrics best represent the QoE perceived by users. Second, how to predict these QoE metrics for a given traffic intensity and to dimension the network.

The first step towards defining QoE and predicting it is to understand how streaming is played. In general, media players at the devices are equipped with a playout buffer that stores arriving packets. As long as there are packets in the buffer, the video is played smoothly. Once the buffer empties, the spacing between packets does not follow the original one. These starvations cause large jitters and are particularly annoying for end users that see frozen images. One feasible way to avoid starvations is to introduce a startup (also called prefetching) delay before playing the stream, and a rebuffering delay after each starvation event. Then after a number of media frames accumulate in the buffer, the media player starts to work. This leads to two important sets of QoE metrics: starvation properties (probability, frequency, etc.) and startup/re-buffering delays.

Once the behavior of media streaming service is understood, the particularity of offering it over wireless networks is considered. Indeed, the wireless channel is subject to a large variability due to fading, mobility, etc. On top of this, it is a shared channel where multiple users are served simultaneously and cell capacity is divided among them. This introduces two variability time scales: flow level (tens of seconds) driven by the departures/arrivals of calls and wireless channel variability time scale (milliseconds) driven by the fast fading.

Related Literature

Video delivery over wireless networks has been studied for more than 15 years. Large attention was given to enhance coding in order to combat errors due to wireless channel variability. Recent work considered flow level performance in cellular networks delivering real time video. [START_REF] Borst | Integration of Streaming and Elastic Traffic in Wireless Networks[END_REF] studied QoS when integrating elastic and video traffic in cellular networks; video QoS was expressed by a blocking rate, while average throughputs and delays represent QoS for elastic traffic. [START_REF] Rong | Performance Evaluation of Cellular Networks Offering TV Services[END_REF] derived the Erlang-like capacity region for a traffic mix including real time video, the aim being to dimension the network for ensuring a target QoS. [START_REF] Karray | Analytical evaluation of QoS in the downlink of OFDMA wireless cellular networks serving streaming and elastic traffic[END_REF] derived the stability region of the network and showed how it is impacted by real-time video traffic. Recently, QoE for new streaming services was studied, taking into account the initial buffering period. [START_REF] Hao Luan | Impact of network dynamics on users' video quality: analytical framework and QoS provision[END_REF] considered a G/G/1 queue where the arrival and service rates are characterized by their first two moments. [START_REF] Liang | Effect of delay and buffering on jitter-free streaming over random VBR channels[END_REF] considered a channel that oscillates between good and bad states following the extended Gilbert model. [START_REF] Parandehgheibi | Avoiding Interruptions a QoE Reliability Function for Streaming Media Applications[END_REF] considered a P2P video streaming based on random linear network coding; this neglects chunk scheduling and peer selection inherent in P2P systems. An M/M/1 queue model has been adopted in [START_REF] Xu | Probabilistic Analysis of Buffer Starvation in Markovian Queues[END_REF], allowing to derive explicit formula for QoE metrics.

As of the tools used for deriving QoE metrics: [START_REF] Hao Luan | Impact of network dynamics on users' video quality: analytical framework and QoS provision[END_REF] adopted a diffusion approximation where the discrete buffer size is replaced with a Brownian motion whose drift and diffusion coefficients are calculated. [START_REF] Liang | Effect of delay and buffering on jitter-free streaming over random VBR channels[END_REF] presented a probabilistic analysis based on an a priori knowledge of the playback and arrival curves. [START_REF] Parandehgheibi | Avoiding Interruptions a QoE Reliability Function for Streaming Media Applications[END_REF] calculated bounds on the playback interruption probability based on the adopted M/D/1 buffer model. Explicit formula of the distribution of the number of starvations has been obtained in [START_REF] Xu | Probabilistic Analysis of Buffer Starvation in Markovian Queues[END_REF], [START_REF] Xu | QoE Analysis of Media Streaming in Wireless Data Networks[END_REF] based on the Ballot theorem.

The above mentioned works on QoE estimation are very useful for catching the impact of variability of the wireless channel due to fast fading or even user's mobility. However, the underlying models fail to capture the large variations due to flow dynamics. For instance, the diffusion approximation in [START_REF] Hao Luan | Impact of network dynamics on users' video quality: analytical framework and QoS provision[END_REF] supposes that the drift and diffusion coefficients are constant over time, which is not true when the number of concurrent flows changes during playback process. The assumption of Poisson packet arrivals in [START_REF] Parandehgheibi | Avoiding Interruptions a QoE Reliability Function for Streaming Media Applications[END_REF], [START_REF] Xu | Probabilistic Analysis of Buffer Starvation in Markovian Queues[END_REF] also fails to take into account these flow dynamics. Note that the analysis of [START_REF] Parandehgheibi | Avoiding Interruptions a QoE Reliability Function for Streaming Media Applications[END_REF] has been generalized to a two-state Markovian arrival process, but this corresponds more to a bursty traffic due to a Gilbert channel model than to flow dynamics.

Our Work

To the best of our knowledge, this paper is the first attempt to assess the impact of flow dynamics on the QoE of streaming. We model the system as two queues in tandem. The first queue, representing the scheduler of the base station, is modeled as a processor sharing queue, while the second represents the playout buffer whose arrival rates are governed by the output process of the base station queue. We consider a static channel (no fast fading) with Constant Bit Rate (CBR) streaming, and derive the prefetching delay distribution and the starvation probability generation function using Partial Differential Equations (PDEs) as well as Ordinary Differential Equations (ODEs) constructed over the Markov process describing the flow dynamics. We next extend the model to include a fast fading channel and show that the impact of flow dynamics is preponderant over the channel variability due to fast fading. Extensive simulations show that our models are accurate enough to be used in QoE prediction. Our analysis also sheds light on the novel QoE enhancement strategies. The results presented here can be used by the base station to "recommend" the prefetching parameters to the media player, and to guide the admission control and the scheduling algorithms.

The remainder of this paper is organized as follows. Section II describes the system model and the QoE metrics. Section III presents the analytical framework of QoE with flow dynamics. Section IV extends the mathematical framework to include the impact of fast fading. The accuracy of our models is validated in Section V. Section VI eventually concludes the paper.

II. PROBLEM DESCRIPTION AND MODEL

In this section, we first describe our motivation, then define the metrics of QoE for media streaming service, and present a queueing model for the playout buffer at a user.

A. Motivation and Network Description

We consider a wireless data network that supports a number of flows.When a new flow "joins" the network, it requests the streaming service from a media server. After the connection has been built, the streaming packets are transmitted through the base station (BS). The streaming flows have finite sizes, which means that a flow "leaves" the network once the transmission completes. Note that each active user cannot watch more than one streams at the mobile device simultaneously. Hence, we use the terms "flow" and "user" interchangeably.

In wireless data networks, a streaming flow may traverse both wired and wireless links, whereas the BS is the bottleneck for the sake of limited channel capacity. In other words, the queue of an active flow is always backlogged at the BS. This assumption holds because most of Internet streaming servers use TCP/HTTP protocols to deliver streaming packets. The TCP protocol in the transport layer exploits the available bandwidth by pumping as more packets as possible to the BS. The BS can easily perform per-flow congestion control to limit TCP sending rate to avoid buffer overflow (a small number of concurrent flows in total). The adaptive coding and modulation in the physical layer, and ARQ scheme at the MAC layer can effectively avoid TCP packet loss. Due to these reasons, we do not consider TCP packet losses in our system.

Streaming flows may experience fast fading and normalized signal-to-noise ratio (NSNR) scheduling is usually adopted to achieve multiuser diversity with the consideration of fairness [START_REF] Choi | Cell-throughput analysis of the proportional fair scheduler in the single-cell environment[END_REF], [START_REF] Song | Asymptotic throughput analysis for channel-aware scheduling[END_REF]. The scheduling duration is commonly around 2ms [START_REF] Bonald | A Queueing Analysis of Data Networks[END_REF]. NSNR selects the user that has the largest ratio of SNR compared with its mean SNR. It is similar to the well-known proportional fair (PF) scheduler in that they both attempt to achieve channel access-time fairness. We consider NSNR instead of PF for two reasons. First, the moments of throughput of PF do not have explicit results, even asymptotic ones (see [START_REF] Song | Asymptotic throughput analysis for channel-aware scheduling[END_REF] and references therein) when the channel capacity is computed according to the Shannon theorem. Second, NSNR needs the knowledge of the average SNR that can be obtained from the history information. The PF scheduler may cause the slow throughput convergence problem if the initial average throughput is not configured wisely.

At the user side, incoming bits are reassembled into video frames step by step. These video frames are played with a deterministic rate, e.g. 25 frames per second (fps) in the TV and movie-making businesses. The size of a frame is determined by the video codec, i.e. a high definition video streaming or a complex video scenario require more bits to render each frame. We consider constant bit-rate (CBR) streaming in which a codec's output data is consumed with the constant rate (i.e. the same size of frames). We highlight the properties of the streaming system briefly to facilitate the mathematical modeling. In our system, there exist three time scales shown in Fig. 1: i) the scheduling duration (e.g. 2ms); ii) playback interval (e.g. 40ms for a video frame rate of 25fps), and iii) duration of flow dynamics (lasting about tens of seconds). The scheduler and the media player do not work at the same granularity of time scale and job size.

B. QoE Metrics

There exist five industry-standard video quality metrics. Authors in [START_REF] Dobrian | Understanding the Impact of Video Quality on User Engagement[END_REF] summarize them into five terms: join time, buffering ratio, rate of buffering events, average bitrate and rendering quality. The first three metrics reflect the fundamental tradeoff in designing the prefetching process. The last two metrics are concerned with source coding. For analytical convenience, we redefine the QoE metrics regarding "prefetching" process.

-Start-up delay: The start-up delay denotes the duration (measured in seconds) between the time that a user initiates a session and the time that the media player starts playing video frames. In the initial prefetching phase, the player starts until the duration of received video reaches the start-up threshold measured in seconds of video segment. The users are impatient to wait for a very long start-up delay. Once the starvation event happens, the player pauses and resumes until the rebuffered video duration reaches the rebuffering threshold. We use the term rebuffering delay to differentiate the rebuffering time from the initial start-up delay.

-Starvation probabilities: When the playout buffer of a user becomes empty before the video has been completely played, we call this event a starvation. The starvation is very annoying to users. We adopt the starvation probability to evaluate the influence of the start-up threshold. In addition, if the rebuffering process is taken into account, we analyze the probabilities of having a certain number of starvations.

Note that the start-up delay and the starvation probabilities can be used to compute the QoE metrics in [START_REF] Dobrian | Understanding the Impact of Video Quality on User Engagement[END_REF]. The expected number of starvations is the sum of the products of the number of starvations and its probability. The expected buffering time equals to the product of the start-up delay in each rebuffering and the mean number of starvation events (including the initial prefetching).

C. Basic Queueing Model of Playout Buffer

We consider a wireless cellular network that supports up to K simultaneous flows. The purpose of admission control is to avoid the overloading of the cell. We make the following assumptions: -Single user type and static channel: We begin with the case where streaming users coexist in a static channel, as this provides an easier route to understand the developed QoE evaluation model. The impact of fast fading is added in section IV. We also consider that all the flows have the same SNR, and hence, in a static channel case, identical throughput. The extension to multiple user classes is presented in Section VI.

-Exponentially distributed video duration: The video duration, measured in seconds, is exponentially distributed with mean 1/θ. Though the exponential distribution is not the most realistic way to describe video duration, it reveals the essential features of the system, and is the first step for more general distributions (i.e. hyper-exponential distribution in Section VI).

-Processor sharing at the BS: The scheduling slot is very small (e.g. ≤2ms in 3G LTE) compared with the service interval between two video frames (e.g. 40ms at 25fps) in the playout buffer. This property enables us to treat the BS as an egalitarian processor sharing queue where all the flows are served simultaneously. Hence, the per-flow throughput, depicted in continuous time, is a deterministic step-wise function of the number of active users in the static channel (e.g. [START_REF] Borst | User-Level Performance of Channel-Aware Scheduling Algorithms in Wireless Data Networks[END_REF]).

-Continuous time playback: The service of video contents is regarded as a continuous process, instead of a discrete rendering of adjacent video frames spaced by a fixed interval. This assumption is commonly used (see [START_REF] Wang | Multimedia Streaming via TCP: An Analytic Performance Study[END_REF]) and is validated by simulations in this work.

We denote by λ the arrival rate of new video streams. Let Bitrate be the playback speed of video streams in bits per-second, and C (in bps) be the capacity of the static wireless channel. Given the exponential distribution of video duration, the file size F (measured in bits) is also exponentially distribution with mean 1/θ F = Bitrate/θ. Therefore, the dynamics of coexisting flows in the cell can be depicted as a continuous time Markov chain with a finite state space.

We concentrate on one "tagged" flow in order to gain the insight of dynamics of the playout buffer. At any time t, the tagged flow sees i other flows in a finite space S := {0, 1, • • • , K -1}. We denote by {I(t); t ≥ 0} the external environment process that influences the throughput of the tagged flow. The environmental change refers to the join of a new flow, or the departure of an existing flow. From the assumption of Poisson flow arrival and exponentially distributed flow size, we can see that {I(t); t ≥ 0} is a homogeneous, irreducible and recurrent Markov process. Let {π i ; i ∈ S} be the stationary distribution of environmental states that will be computed in the following sections. The throughput of the tagger user is b i := C Bitrate•(i+1) in seconds of video contents at state i. Let N e (t) be the number of changes in the environment by time t. Denote by A l the time that the l th environmental change takes place with A 0 = 0 and by I l := I(A l ) the state to which the environment changes after time A l .

We denote by Q(t) the length of playout buffer measured in seconds of video contents at time t. In the prefetching phase,

Q(t) is expressed as Q a (t) = Ne(t) l=1 b I l (A l -A l-1 ) + b I Ne(t) (t -A Ne(t)
). Denote by q a the start-up threshold. The start-up delay T a is defined as

T a = inf{t ≥ 0|Q a (t) ≥ q a }. The cumulative distribution of T a is expressed as Ψ i (t; q a ) = P{T a < t|I(0) = i} if the tagged flow is in state i upon arrival.
When the media player starts the playback, the queueing process {Q(t); t ≥ 0} is given by

Q b (t) = q -t+ Ne(t) l=1 b I l (A l -A l-1 )+b I Ne (t) (t-A Ne(t)
), if the time axis starts at the instant of playing. Define c i := b i -1 for all i ∈ S. Define

T b := inf{t ≥ 0|Q b (t) < 0} (1)
to be the time of observing empty buffer. Denote by T e (T e < ∞) the completion time of downloading of the tagged flow. If T b is less than T e , a starvation event happens at the playout buffer. Then, the ultimate starvation probability is computed as

W i (q a ) = P{T b < T e |I(0) = i, Q b (0) = q a }.
When the playback begins at state i, and stops upon an empty queue.

The ultimate starvation probability is the weighted sum of starvation probabilities at all the ergodic entry states. Next, we proceed the computations using an approach that is inspired by the ruin analysis in actuarial theory [START_REF] Asmussen | Ruin Probabilities[END_REF], [START_REF] Lu | On the probability of ruin in a Markov-modulated risk model[END_REF].

III. IMPACT OF FLOW LEVEL DYNAMICS ON QOE

In this section, we model the starvation probability and the prefetching delay in a static channel where the media flows join and leave the system dynamically. The key idea is to investigate the queueing process of one "tagged" flow on the basis of differential equations.

A. Markov models of flow dynamics

Our purpose here is to construct two Markov chains to characterize the dynamics of the number of active flows. The first one models flow dynamics before the "tagged" flow joins in the network. Based on this Markov process, we can compute the stationary distribution of the number of active flows observed by the "tagged" flow at the instant when it is admitted. The second one describes the flow dynamics after the tagged flow is admitted. This Markov process enables us to investigate how the playout buffer of the tagged user changes.

We first look into the flow dynamics before the tagged flow joins. When the normalized SNR scheduler is used, the perflow throughput is proportional to the reciprocal of flow population. Given the Poisson arrival rate and the exponentially distributed service time, we can model the dynamics of flow population as a finite-state Markov chain Z a := {0, 1, • • • , K} (illustrated in our technical report [START_REF] Xu | Flow Level QoE of Media Streaming in Wireless Networks[END_REF]). The transition rate from i to i -1 is µ i := Cθ F . Note that the network capacity is a constant in the static channel. Hence, we let µ i = µ for i = 1, 2, • • • , K and µ 0 = 0. Define ρ := λ µ to be the load of the channel. Let z a i be the stationary probability that there exist i flows. We give the expression of z a i (i ∈ S ∪ {K}) directly because it is easy to compute: z a 0 = 1-ρ 1-ρ K+1 , and

z a i = ρ i (1-ρ) 1-ρ K+1 ∀i = 1, • • • , K.
The tagged user cannot be admitted at state K due to the admission control at the BS. Therefore, if it joins in the network successfully, it will observe i other flows with the probability π i := (i+1) so that there has ν i := iCθ F (i+1) = i i+1 µ for all i ∈ S.For the simplicity of notations, we denote by λ i the transition rate from state i to i + 1. It is obvious to have λ i = λ for all i = K-1 and λ K-1 = 0.

z a i 1-z a K = ρ i (1-ρ) 1-ρ K ,

B. Modeling prefetching delay distribution

We want to know how long time the tagged user needs to wait in the prefetching phase. Recall that q a is the start-up threshold. The prefetching time is only meaningful when the video duration is longer than q a . In the prefetching phase, because the playout buffer does not serve video frames, the queue length of the tagged flow evolves in an infinitesimal time interval [0, h] with h(> 0) Q(t + h) = Q(t) + b i h. The distribution of the prefetching time is difficult to solve directly. We resort to the following duality problem: DUALITY PROBLEM: What is the starvation probability by time t if the queue is depleted with rate b i (i ∈ S) and the duration of prefetched contents is q a ? In the duality problem, the queue dynamics in [0, h] is modified as

Q(t+h) = Q(t)-b i h. We define U i (q, t) (∀i ∈ S)
to be the probability of starvation before time t, conditioned on the entry state i and the initially prefetched content q. We use differential equations to obtain U i (q, t). In the infinitesimal time interval [0, h], there are four possible events -no change of the concurrent flows; -arrival of one flow; -departure of one flow (not the tagged one); -occurrence of more than one events. Conditioned on the events occurred in [0, h], we have

U i (q, t) = (1 -λ i h -ν i h)U i (q -b i h, t -h) +λ i hU i+1 (q -b i h, t -h) +ν i hU i-1 (q -b i h, t -h) + o(h), ∀i ∈ S. (2)
The above equations yield for i ∈ S

1/h • (U i (q, t) -U i (q -b i h, t -h)) = λ i U i+1 (q-b i h, t-h) -(λ i + ν i )U i (q -b i h, t -h) +ν i U i-1 (q-b i h, t-h) + o(h)/h. (3)
When h → 0, the left side of eq.( 3) is the partial differentials of U i (q, t) over q and t. In other words, eq.( 3) yields a set of linear partial differential equations (PDEs)

(∂U i /∂t) = -b i • (∂U i /∂q) -(λ i + ν i )U i (q, t) +λ i U i+1 (q, t) + ν i U i-1 (q, t), ∀i ∈ S, (4)
with the initial condition: U i (q, 0) = 0, ∀ q > 0; and the boundary conditions at both sides ∀ t ≥ 0: U i (0, t) = 1, and lim q→∞ U i (q, t) = 0, ∀t ≥ 0. The initial condition means that the starvation cannot happen at time 0 for q > 0.

The right-side boundary condition says that the starvation will not happen before t if the initial prefetching is large enough. Hence, the c.d.f. of start-up delay is the solution of linear PDEs by letting q be q a . We next analyze the probability that the prefetching process starts at state i and ends at state j, for all i, j ∈ S. Define V i,j (q; q a ) := P{I(T a ) = j|I(0

) = i, Q(0) = q}. ( 5
)
We can use the approach of obtaining U i (q, t) to solve V i,j (q; q a ). Note that we now back to the queueing dynamics

Q(t + h) = Q(t) + b i h.
In the time interval [0, h], there exists for all i, j ∈ S V i,j (q; q a ) = (1 -λ i h -ν i h)V i,j (q + b i h; q a ) +λ i hV i+1,j (q+b i h; q a ) + ν i hV i-1,j (q+b i h; q a ). [START_REF] Xu | QoE Analysis of Media Streaming in Wireless Data Networks[END_REF] It is easy to see that V i,j (q; q a ) is the solution of the following differential equation b i Vi,j (q; q a ) = (λ i + ν i )V i,j (q; q a ) -λ i V i+1,j (q; q a ) -ν i V i-1,j (q; q a ), ∀i, j ∈ S,

with the boundary conditions for all i, j ∈ S: V i,j (q a ; q a ) = 1 if i = j; V i,j (q a ; q a ) = 0, otherwise. We interpret the boundary condition in the following way. If there exist I(0) = i and Q(0) = q a , the prefetching duration is 0 and the prefetching process ends at state i. Hence, V i,j (q a ; q a ) is 1 iff i equals to j. We use Laplace Transform (LT) to solve eq.( 7). Define a matrix M V as below:

     -λ0+ν0 b0 λ0 b0 0 • • • 0 0 ν1 b1 -λ1+ν1 b1 λ1 b1 • • • 0 0 • • • • • • • • • • • • • • • • • • 0 0 • • • • • • ν N -1 b N -1 -ν N -1 +λ N -1 b N -1     
Define 1 j to be a column vector in which the j th element is 1 and all other elements are 0. Taking LT of eq.( 6), we obtain

b i sV i,j (s) = (λ i + ν i )V i,j (s) -λ i V i+1,j (s) -ν i V i-1,j (s) + b i V i,j (0; q a ), ∀i, j ∈ S, (8)
where V i,j (s) is the laplace transform of V i,j (q; q a ). Thus, we obtain the vector of probabilities that the prefetching process ends at state j

V j (0; q a ) = [L -1 {(sI + M V ) -1 }(q a )] -1 • 1 j , ∀j ∈ S, (9)
where L -1 {} denotes the inverse LT. For the computation of eq. ( 9), interested readers can refer to the technical report [START_REF] Xu | Flow Level QoE of Media Streaming in Wireless Networks[END_REF].

C. Modeling starvation probability

The modeling of starvation probabilities should take into account the departure of the tagged flow. Recall that the CTMC in Fig. 2 assumes the persistent tagged flow, which is not suitable for the playback process. Before solving the starvation probabilities, we first modify the original CTMC by adding an absorbing state A (illustrated in our technical report [START_REF] Xu | Flow Level QoE of Media Streaming in Wireless Networks[END_REF]). The state A denotes the event that the tagged flow completes its downloading. Because of the exponentially distributed video duration, the transition from state i to state A is Poisson. Denote by Z c := {0, 1, • • • , K-1, A} this absorbing Markov chain. Denote by ϕ i the transition rate from state i to A. At state i, the bandwidth of a flow is C i+1 , resulting in ϕ i := µ i+1 . Define c i := b i -1. The queue length of the tagged flow changes in an infinitesimal interval h according to the rule

Q(t + h) = Q(t) + c i h. If c i > 0,
the bandwidth is sufficient for continuous playback of the tagged flow and i other flows. For mathematical convenience, we suppose that q is 0 -if buffer starvation happens.When the tagged flow enters the absorbing state, it has downloaded the whole file with a nonempty playout buffer. Thus, the starvation probability at state A is 0 for any q ≥ 0. Let W i (q) be the starvation probability with q seconds of contents in the playout buffer at state i. We derive a system of ordinary differential equations for W i (q). In an infinitesimal interval [0, h], there are five possible events: -all four events the same as those in the prefetching stage; -the tagged flow entering the absorbing state.

The above events give rise to the following bunch of equations

W i (q) = (1 -(λ i + µ i )h)W i (q + c i h) (10) +λ i W i+1 (q + c i h) + ν i W i-1 (q + c i h) + o(h).
When h → 0, we obtain

c i Ẇi (q) = (λ i + µ i )W i (q) -λ i W i+1 (q) -ν i W i-1 (q). (11)
Similarly, we solve the above ODEs using LT. Define Ŵi (s) as the LT of W i (q), and Ŵ(s) as the vector

{ Ŵ0 (s), Ŵ1 (s), • • • , ŴK-1 (s)}. They yield for all i ∈ S c i (s Ŵi (s) -W i (0))=(λ i + µ i ) Ŵi (s)-λ i Ŵi+1 (s)-ν i Ŵi-1 (s).
We write the above equations in matrix form:

(sI + M W ) • Ŵ(s) = W(0) ( 12
)
where M W is expressed as the following:

     -λ0+µ0 c0 λ0 c0 0 • • • 0 0 ν1 c1 -λ1+µ1 c1 λ1 c1 • • • 0 0 • • • • • • • • • • • • • • • • • • 0 0 • • • • • • ν N -1 c N -1 -µ N -1 +λ N -1 c N -1     
and W(0) is the vector of initial state. In order to compute Ŵi (s), we first need to find out W i (0) (∀i ∈ S), the starvation probability without prefetching. However, eq.( 12) does not give explicit expressions for W i (0). Theorem 1: Let adj(sI + M W ) be the adjugate matrix of sI + M W . {W i (0)} are solved by

-W i (0) = 1 if c i < 0, ∀i ∈ S ; -W i (0) = 0, ∀i ∈ S if c K-1 ≥ 0; -adj(sI+M W )•W(0)| s=ξ l =0, ∀ l = 0, 1, • • • , m, m < K-1, if {ξ l } are positive roots of the determinant |sI+M W |.
Proof: Please refer to the technical report [START_REF] Xu | Flow Level QoE of Media Streaming in Wireless Networks[END_REF].

Once the vector W(0) is solved, the LTs of starvation probabilities can be computed by Ŵ(s) = [sI + M W (s)] -1 • W(0). Using inverse LT, we obtain the starvation probability W(q) := L -1 {[M W (s)] -1 }(q) • W(0). Substituting q by q a , we derive the starvation probabilities W(q a ).

Next, we build a bridge to interconnect the prefetching threshold and the starvation probability function W i (q). For a given prefetching threshold q a , the starvation event takes place only when the video duration T video is longer than q a . Thus a flow with T video > q a can be regarded as a tagged flow. When the prefetching process is finished, the tagged flow enters the playback process. Conditioned on the distribution of entry states {π}, the distribution of the states that the playback process begins (or the prefetching process ends) is computed by π•V(0; q a ). Then, the starvation probability with the prefetching threshold q a is obtained by [START_REF] Dobrian | Understanding the Impact of Video Quality on User Engagement[END_REF] strength is no longer a constant at different scheduling slots. To explore the multiuser diversity gain, the base station adopts the normalized SNR scheduling algorithm for allocating time slots to coexisting flows.

P s (q a ) = P{T video > q a } • π • V(0; q a ) • W(q a ) = exp -θq a • π • V(0; q a ) • W(q a ).
We begin with the scenario with a fixed population of i users (or flows) served by a single base station. In each slot, the users measure their channel qualities and feedback them to the BS. Based on the channel quality indications, the BS transmits to only one of the users every slot. Denote by γ j,n the instantaneous signal to noise ratio (SNR) of user j, (1 ≤ j ≤ i), at slot n. As stated in most of previous work, we assume that all the users experience Rayleigh fast-fading. Denote by γj the average SNR of user j. Then, the received SNR of user j is an exponentially distributed random variable with the following probability density function g j (γ) = 1 γj exp(-γ γj ). The NSNR scheduler selects the user that has the highest relative SNR for transmission, j * n = max j {γ j,n /γ j , j = 1, 2, • • • , i}, where j * is the scheduled user at slot n. In this section, we consider the case of homogeneous average SNRs (i.e. γj = γ for all j). Therefore, the NSNR scheduler is equivalent to the maximum sum rate (MSR) scheduler that gives the largest per-user throughput.Since the SNRs of different users are independently distributed, the scheduled SNR, denoted by γ * , has the following probability density function [START_REF] Chang | Cross-layer QoS Analysis of Opportunistic OFDM-TDMA and OFDMA Networks[END_REF] g

* (γ) = i γ exp(-γ γ ) 1 -exp(-γ γ ) i-1 . Denote by f (γ) the
data rate of a user with the SNR γ. Here, f (•) can be a linear function in the low-SNR regime and a logarithmic function in the high SNR regime if the modulation scheme is continuous. For discrete modulations, f (•) is a step function of γ. Without loss of generality, we let f (γ) = log 2 (1 + γ).

Analysis of throughput process.

The fast fading along with NSNR scheduling brings variation of bit arrivals to the receiver. Here, our purpose is to show that the throughput variance may have negligible influence on the computation of starvation behaviors. In what follows, we calculate the mean throughput and its variance measured in video duration. To achieve this goal, we shall obtain the mean throughput and its variance measured in bits first.

Denote by r * i the transmission rate of the user with the best SNR at a slot in each Hz when there are i active flows in the cell. Denote by r i the transmission rate to one particular flow at a slot per Hz. Given the assumption that all the flows have the same average SNR, each flow has the equal probability of being scheduled. Hence, we can see

r i := r * i w.p. 1 i ; 0 w.p. i-1 i .
For the r.v. r * i , its mean and variance are computed by

E[r * i ] = ∞ 0 f (γ) • g * (γ)dγ, Var[r * i ] = ∞ 0 f (γ) 2 • g * (γ)dγ -(E[r * i ]) 2 .
The above equations yield:

E[r i ] = 1 i E[r * i ], Var[r i ] = E[r 2 i ] -(E[r i ]) 2 = 1 i E[(r * i ) 2 ] -1 i 2 (E[r * i ]) 2 = 1 i Var[r * i ] + (E[r * i ]) 2 ( 1 i -1 i 2 )
. Denote by D s the duration of scheduling slot (usually 2ms), and by B the width of wireless spectrum in Hz. Then, the mean and the variance of per-flow throughput measured in the duration of video contents are

B•Ds•E[ri] Bitrate and ( B•Ds Bitrate ) 2 •
Var[r i ] respectively in one slot. Let R i be the r.v. of per-flow throughput in one second that is measured by the duration of video contents. In one second, the total throughput of a flow at one Hz is the sum of throughput in 1 Ds slots. Therefore, the r.v. R i is the sum of 1 Ds i.i.d. r.v.s corresponding to the per-slot throughput. We can express the mean and the variance of R i as follows:

E[R i ] = 1 Ds • B•Ds•E[ri] Bitrate = B•E[r * i ] i•Bitrate , Var[R i ] = 1 Ds • ( B•Ds Bitrate ) 2 • Var[r i ] = 1 i Var[r * i ]+(E[r * i ]) 2 ( 1 i -1 i 2 ) • B 2 •Ds Bitrate 2 .
In general, the frequency width B is 1∼5 MHz, the bit-rate is usually greater than 200 Kbps, and D s equals to 0.002s. Then, Var[R i ] is usually at the order of 10 -2 . If starvation happens at state i, E[R i ] is usually less than 1, which means that B Bitrate needs to be small. However, the small B Bitrate results in the small variance Var[R i ]. This is to say, if the variance of bit arrival process is large, there might not exist starvations. On the contrary, if the starvations appear, the variance is usually small so that its impact on the starvation is negligible. For this reason, we directly use the framework without diffusion approximation to model the streaming QoE in a fast fading channel.

Markov model of flow dynamics

To analyze the interaction between NSNR scheduling and the flow dynamics, a fluid-level capacity model is required. When the average SNR of all active users are the same, the per-flow throughput in each slot is i.i.d. and only depends on the quantity of flows (see section VI.A in [START_REF] Xu | Flow Level QoE of Media Streaming in Wireless Networks[END_REF]). Given the exponentially distributed video size, we can model the flow dynamics as a Markov process.

The Markov processes {Z a }, {Z b }, and {Z c }, contain transitions rates such as µ i , ν i and ϕ i . However, it is not direct to feed the parameters of this section into the above Markov processes. In {Z a }, state i refers to the number of flows in the system. The departure rate is computed by µ i = iθE[R i ] for i ∈ S ∪ {K}, recalling that E[R i ] is average per-user throughput in video duration per second. It is easy to obtain the stationary distribution of having i flows by

z a i = λ i i l=1 µ l 1 + K j=1 λ j j l=1 µ l -1 , ∀i = 0, • • • , K,
(with the convention that over an empty set is 1). When a tagged user joins in the system and is also admitted, it observes i other flows with the following stationary distribution {π} :

π i = z a i 1 -z a K = λ i i l=1 µ l 1 + K-1 j=1 λ j j l=1 µ l , ∀i ∈ S.
The Markov processes {Z b } and {Z c } are conditioned on the existence of the tagged flow. At state i, the per-user throughput is E[R i+1 ] because there are i flows plus the tagged one. Hence, the transition rate ν i is computed by

ν i := iθ • E[R i+1 ] for all i ∈ S. The transition rate ϕ i is expressed as ϕ i := θ•E[R i+1 ].
Define μi as the total departure rate at state i that has

μi := ϕ i + ν i = (i + 1)θE[R i+1 ] = µ i+1 , (20)
in the presence of the tagged flow. The constants b i and c i are obtained by

b i = E[R i+1 ] and c i = b i -1, ∀i ∈ S. (21)
Substituting the above parameters to the framework in section III, we can derive the approximated QoE metrics in a fast fading channel with flow dynamics.

V. SIMULATION

In this section, we compare the numerical experiments with the developed framework using MATLAB. Our models are shown to match the simulations accurately.

A. Flow dynamics without fast fading

We consider a network with maximum number of ten simultaneous streaming flows and the capacity of 2.5Mbps. Flows arrive to the network with a Poisson rate λ = 0.12. Let the video duration be exponentially distributed with the mean 60 seconds. Then, there have µ = 0.1302 and ρ = 0.9216 at the playback rate 360Kbps, and µ = 0.0868 and ρ = 1.3824 at the playback rate 480Kbps. The simulation lasts 5 × 10 5 s. Starvation probabilities: In this set of experiments, we will illustrate the overall starvation probability, the starvation probabilities when the playback process begins at different states, as well as the p.g.f. of starvation events. Fig. 4 shows the overall starvation probabilities with different settings of the start-up threshold. When it increases from 0 to 20s of video contents, the starvation probability decreases. The higher playback rate (e.g. 480Kbps) incurs larger starvation probabilities in comparison with the lower playback rate (e.g. 360Kbps). Our mathematical models match the simulations very well. Figure 5 compares the starvation probabilities when the playback process begins at different states. A higher state refers to more coexisting flows (or congestions), and hence causing a larger starvation probability. Note that the arrival rates at state 7 and 9 are less than 360Kbps. Without prefetching, the starvation event happens for sure. We evaluate the probabilities of having one or two starvations in the whole procedure in Fig. 6. For clarity, we choose the same value for startup and re-buffering thresholds. The starvation probabilities increase in the beginning and decrease afterwards when q a (or q b ) increases from 0 to 30s. This is because there are many starvations with very small start-up threshold and few starvations with very large start-up threshold. Our analytical model predicts the starvation probabilities accurately. Start-up delay: We illustrate the distribution of start-up delays in Fig. 7. The start-up threshold is set to 10s. We highlight the c.d.f. curves when the tagged flow sees 3, 5, and 7 other flows respectively after entering the network. We use MATLAB PDE function pdepe to compute the model in eq.( 4) numerically. Fig. 7 demonstrates accurate estimation of start-up delay in the simulation. When the cumulative probability is close to 1, the PDE model oscillates slightly. This is because U i (0, 0) is discontinuous in the boundary and initial conditions. The elimination of this oscillation demands an iterative solution of linear PDEs, which will be our future work.

B. Rayleigh Fading Channel

Consider a wireless channel with frequency width of 1MHz. The average SNRs of users is 5dB. The base station allows at most 10 flows simultaneously, and schedules the transmission to one of them in every slot of duration 0.002s. The video duration is exponentially distributed with the mean of 90 seconds and the video bit rate is chosen to be 480Kbps. Then, the mean throughput are {3.5749, 2.3702, 1.7844, 1.4369, 1.2061, 1.0412, 0.9174, 0.8207, 0.7432, 0.6794} times the playback rate at states from 0 to 9. In other words, the mean throughput at states 6∼9 are insufficient to support the continuous playback. The variances at all states are {0.0083, 0.0144, 0.0144, 0.0134, 0.0124, 0.0114, 0.0105, 0.0098, 0.0091, 0.0086}, which are small enough. We consider two flow arrival rates, λ = 0.07 and λ = 0.09. For λ = 0.07, the traffic load ρ is greater than 1 at states 0∼5 and less than 1 at states 6∼9. For the latter case, there have ρ > 1 at all the states. Each set of simulation lasts 2 × 10 7 time slots.

In Fig. 8 we compare the starvation probabilities measured from a Rayleigh fading channel, and those computed from the model without considering throughput variation. The simulation matches the model quite well, which means that the flow-level dynamics have a dominant impact on the playback interruption, while the impact of throughput variation due to Rayleigh fading is negligible. In Fig. 9 we examine the starvation probabilities when the playback process begins at different states. We test two start-up thresholds, q a = {5, 10}, and two flow arrival rates, λ = {0.07, 0.09}. One can observe that the starvation probabilities do not differ much in high states (e.g. 8 and 9). However, the starvation probabilities in the states with mean throughput around 1 are distinguishable, in which state 6 is an example. With λ = 0.09, a tagged flow sees the congested network (more other flows) with a higher probability, and also encounters a higher probability of starvation afterwards.

VI. CONCLUSION

In this work, we developed an analytical framework to compute the QoE metrics of media streaming service in wireless data networks. Our framework takes into account the dynamics of playout buffer at three time scales, the scheduling duration, the video playback variation, as well as the flow arrivals or departures. We show that the proposed models can accurately predict the distribution of prefetching delay and the probability generating function of buffer starvations. The analytical results demonstrate that the flow dynamics have dominant influence on QoE metrics compared to the jittering in the throughput due to the fast fading. Our result can be used by the streaming publishers to tune the video prefetching, and by the network operator to configure admission control and scheduling algorithms.

As of future works, we aim at studying the QoE with the variation of playback process. The more heterogeneous case where users of different radio conditions or different service types share the network resources is also left as an interesting future extension.

Fig. 1 .

 1 Fig. 1. Illustration of 3 different time scales

  ∀i ∈ S. After the tagged flow joins in the network, the Markov process Z a has been altered. The states are the number of flows observed by the tagged user, and the transition rates are conditioned on the presence of the tagged flow. Therefore, we model the flow dynamics observed by the tagged flow through a finite-state Markov chain Z b := {0, 1, • • • , K-1} in Fig.2. Denoted by ν i the transition rate from state i to i-1. The per-flow throughput at state i is C

Fig. 2 .

 2 Fig. 2. Flow dynamics observed by tagged flow

Fig. 4 .Fig. 6 .Fig. 7 .Fig. 8 .Fig. 9 .

 46789 Fig. 4. Overall starvation probability VS startup threshold

D. Modeling P.G.F. of starvation events

When a starvation event happens, the media player pauses until q b seconds of video contents are re-buffered. A more interesting but challenging problem is how many starvations may happen in a streaming session. In this section, we come up with an approach to derive the probability generating function of starvation events.

We define a path as a sequence of prefetching and starvation events, as well as the event of completing the downloading. Obviously, the probability of a path depends on the number of starvations. We illustrate a typical path with L starvations in figure 3 that starts from a prefetching process and ends at a playback process. We denote by I A l the beginning state of the l th prefetching, by I B l the beginning state of the l th playback, and by I e the end of downloading. The end of a prefetching process is exactly the beginning of a playback process. The end of a playback process is also the beginning of a subsequent prefetching process if the video has not been downloaded completely. This path contains a sequence of events happening at the states

, I e }. The process between I A l and I B l is the l th prefetching process, while that between I B l and I A l+1 is the l th playback process, (1 ≤ l ≤ L). The first starvation takes place at the instant that the second prefetching process begins. The starvation event (e.g. I B l , 1 ≤ l ≤ L) cannot happen at the state i that has c i ≥ 0. The sample path in figure 3 demonstrates a roadmap to find the p.g.f. of starvation events. We need to compute the transition probability along the path with all possible states. Recall that the transition probabilities from state I A l to I B l have been computed in section III-B. The only missing part is the transition probabilities from state I B l to I A l+1 . Denote by X i,j (q) the probability that a playback process starts at state i and meets with the empty buffer at state j with the prefetching threshold q. Define a matrix X(q) := {X i,j (q); i, j ∈ S}. Denote by X j (q) the vector of probabilities that the starvation takes place at state j with the prefetching threshold q, i.e. X j (q

0)] T be the vector of those probabilities without the prefetching. Using the same argument, we get the differential equation of X i,j (q), ∀i, j ∈ S,

The solution of eq.( 14) is directly given by

The computation of X j (q) requires the knowledge of the boundary condition X j (0).

Theorem 2:

The probability that the playback process starts at state i and ends at state j without prefetching is computed by

• The third item of Theorem 1. Proof: The proof is the same as that of Theorem 1.

When replacing q by q a , we obtain the probability X ij (q a ) that the first starvation happens at state j with i other flows observed by the tagged flow at the beginning of the playback process. The starvation probability in a rebuffering process is calculated by X ij (q b ), given the rebuffering threshold q b .

The probability of having L starvations can be expressed as the product of the probabilities from the first prefetching to the last playback. The probability vector from I A 1 to I B 1 is obtained by

The probability vector from

Recall that the starvation happens at state I A 2 , and the rebuffering process ends at state I B 2 with the prefetched video duration q b . We next compute the probability of having only one starvation denoted by P 1starv . The possible paths include

The first part of P 1starv refers to the case that the remaining video duration is less than the rebuffering threshold q b . The second part refers to the case that the remaining video duration is longer than q b and there is no starvation after the rebuffering process. P 1starv

Here, the expression (1 -W(q b )) is the probability I A 2 → I e in the first path and the expression (1 -W(q b )) is that of I B 2 → I e in the second path. Similarly, we can deduce the probability of having L(L > 1) starvations recursively by

Though the expression in eq.( 19) looks complicated, it only involves duplicated products of matrices with dimension K that can be calculated easily.

IV. IMPACT OF FAST FADING ON QOE

This section models the starvation behavior of CBR streaming when users experience fast channel fading. We compute the first two moments of bit arrival process and show how these parameters can be fed into our analytical framework.

Network description. Due to the change of radio condition (e.g. user mobility, or a car passing by the user), the signal