Speed of convergence to the quasi-stationary distribution for L\'evy input fluid queues

In this note we prove that the speed of convergence of the workload of a L\'evy-driven queue to the quasi-stationary distribution is of order $1/t$. We identify also the Laplace transform of the measure giving this speed and provide some examples.


Introduction
In this paper, we consider a storage system with Lévy netput. In other words, the workload process {Q(t), t ≥ 0} is a spectrally one-sided Lévy process X(t) that is reflected at 0: We assume that the drift of the process X(t) is negative; that is, we have EX(1) < 0. This stability condition guarantees the existence of a stationary distribution π of Q, which by virtue of Reichs identity' can be expressed in term of the all-time supremum: (2) π(x) = P sup t≥0 X(t) ≤ x .
In the sequel, we consider the initial distribution Q(0) sampled from this steady-state distribution which is indicated by adding the subscript π to the probability measure P and to the associated expectation E. Now, let T denote the busy period; that is T = inf{t ≥ 0 : Q(t) = 0}.
We will further consider the Yaglom limit lim t→∞ P π (Q(0) ∈ dx, Q(t) ∈ dy | T > t) := µ(dx, dy), where the convergence is to be understood in the weak sense. Yaglom limits are a probability measure and a particular case of the quasi-stationary (QS) distribution, which is and invariant distribution for the process conditioned on non-extinction; that is, we condition on the event that the process survives some killing event (e.g. related with exiting from some subset of possible values). Yaglom [30] was the first to explicitly identify QS distributions for the subcritical Bienaymé-Galton-Watson branching process. This result has been generalized in the context of the continuous-time branching process and the Fleming-Viot process; see [1,9,18]. Similar results were also derived for Markov chains on positive integers with an absorbing state at the origin; see Seneta and Vere-Jones [27], Tweedie [29], Jacka and Roberts [13] and the bibliographic database of Pollet [25]. Recently, Foley and McDonald showed that Yaglom limits may depend on the starting state [10].
Work on QS distributions has been very extensive. Martinez and San Martin [22] analyze the Brownian motion with drift exiting from the positive half-line, complementing the result for random walks obtained by Iglehart [12]. Later, QS laws have been studied for various Lévy processes. Kyprianou [15] found the Laplace transform of the QS distribution for the workload process of the stable M/G/1 queue with service times that have a rational moment generating function. Kyprianou and Palmowski [17] identified the quasi-stationary distribution associated with a general light-tailed Lévy process. Haas and Rivero [26] found (after appropriate scaling) the QS distribution when the Lévy process under study has a jump measure with a regularly varying tail. The speed of convergence (in total variation) to the quasi-stationary distribution for population processes has been studied in [5]. Finally, Mandjes et al. [21] derived the QS distribution of the workload process Q(t). This paper builds upon [21].
A contribution of this paper lies in proving that the speed of convergence to the quasi-stationary distribution is surprisingly slow (of order 1/t). We also identify a measure ξ(dx, dy) (which we call second-order quasi-stationary measure), such that We hence prove the conjecture posed in Polak and Rolski [24], which proved the above statement for a birth-death process by using an asymptotic expansion of a transition function and certain properties of Bessel functions. In this paper, we suggest new method, which relies on a refined Tauberian-type expansion of the Laplace transform. We also analyze in detail the M/M/1 queue and a Browniandriven queue.
If we want to simulate the quasi-stationary distribution directly from definition (1), then the result stated in (3) shows that the speed of such a simulation is very slow. Still, in the literature there are papers giving other efficient algorithms of simulation of quasi-stationary measures; see e.g. Blanchet et al. [4] and references therein.
The main result in (3) contrasts the typical results derived for the regular stationary distribution of Markov processes where, in most of the cases, the rate is exponential. More precisely, for many models the distance between the distribution of the stochastic process at time t and its stationary distribution decays exponentially fast in t. The typical distances used are the total variation distance, the separation distance, and the L 2 distance. The classical results concern mainly Markov chains and use Perron and Frobenius theory, renewal equations or the coupling method; see e.g. [6,8,14,19,20] and references therein. Another method concerns Harris recurrent Markov processes and it is based on the construction of a special Lyapunov function and then the application of Foster-Lyapunov criteria; see e.g. [2,23,28]. All the above-mentioned methods though are different from the one used in this paper, which is based on expansions of Laplace transforms.
The paper is organized as follows. In next section, we introduce the notation and basics facts that are used later. In Section 3, we present the main results. The central step for the proof of the main results is given in Section 4. Finally, the last section provides some examples.

Preliminaries
We follow [16] for definitions, notations and basic facts on Lévy processes. Let X ≡ (X(t)) t be a spectrally negative Lévy process, which is defined on the filtered space (Ω, F, {F t } t≥0 , P) with the natural filtration that satisfies the usual assumptions of right continuity and completion. We define P x as P x (X(0) = x) = 1 and P 0 = P; similarly, E x is the expectation with respect to P x . We denote by Π(·) the jump measure of X, which is supported by spectral negativity in the non-positive half-line; in other words, jumps are non-positive. We define the Laplace exponent ψ(η) by (4) Ee ηX(t) = e tψ(η) , for η ∈ R such that the left hand side of (4) is well-defined (which holds at least for η ≥ 0). We denote by Φ(s) := inf{η ≥ 0 : ψ(η) > s} the right inverse of ψ; see [16] for details.
Dual process. We also consider the dual processX t = −X t with jump measureΠ (0, y) = Π (−y, 0). Note thatX(t) is a spectrally positive process having only non-negative jumps. Characteristics ofX are indicated by using the same symbols as for X, but with a 'ˆ' added. In particular, We skip the symbol 'ˆ' for Q(t) and hence for T and all quantities related to Q as it will be clear from the context if a statement concerns the spectrally negative or the spectrally positive case.

Asymptotic expansions. Consider a function
for z ∈ C be its Laplace transform. Consider singularities off (z); among these, let a 0 < 0 be the one with the largest real part. Notice that this yields the integrability of for some (and then any) a > a 0 . In this paper, we need the following Tauberian theorem, found in Doetsch [7, Theorem 37.1], where the behaviour of the Laplace transform around the singularity a 0 < 0 plays a crucial role. First recall the concept of the W-contour, centered at a 0 , with a half-angle of opening π/2 < ψ ≤ π, as depicted on [7, Fig. 30, p. 240] ; also, for the purposes of our problem, G a 0 (ψ) is the region between the contour W and the line ℜ(z) = 0. More precisely, where arg z is the principal part of the argument of the complex number z. In the following theorem, conditions are identified that provide an asymptotic expansion of the Laplace transform.

Definition 2 (Semiexponentiality)
. A function f is said to be semiexponential if for some 0 < φ ≤ π/2, there exists a finite and strictly negative function γ(ϑ), called the indicator function, defined as the infimum of all a ∈ R such that f (e iϑ r) < e ar for all sufficiently large r; here −φ ≤ ϑ ≤ φ and sup γ(ϑ) < 0.
Quasi-stationary distribution. We define a new probability measure P η x , needed in the sequel, by the exponential change of measure The quasi-stationary distribution of the workload process Q(t) in the stationary regime was identified in [21].
Theorem 3 (Mandjes et. al [21]). (i) If X is a spectrally negative Lévy process satisfying conditions (SN), then is a renewal function related with the downward ladder height process of X considered on the shifted P ϑ * x probability measure. (ii) If X is a spectrally positive Lévy process satisfying conditions (SP), then is a renewal function related with the upward ladder height process of X considered on the shifted P −ϑ * x probability measure.
In addition, when X is spectrally negative the bivariate Laplace transform of the quasi-stationary measure. From Proposition 3 stated below (see also [21]), we obtain the following equivalent expressions.
The key component of the proof of this corollary is based on Wiener-Hopf factorization, from which master formulas can be derived (given below), and on either some expansion theorems (see [3] and [17]) or some Tauberian-type theorems.
Master formulas. Recall that Q(t) given in (1) is a workload process with stationary distribution (2) and busy period T . We define now the double Laplace-Stieltjes transform: In [21], the following representations of L were derived.

Main results
We state now the main results of this paper. Define the constants We start with the following expansion for the double Laplace-Stieltjes transform.

Now, straightforward calculations give
which completes the proof due to the definition of the measure ξ.
Remark 8. Formally starting from Proposition 5, then using generalization of Proposition 6 in the next step combined with Theorem 1, will give more terms of the expansion of E π [e −αQ(0)−βQ(t) |T > t], and hence a longer expansion of the measure P π (Q(0) ∈ dx, Q(t) ∈ dy|T > t) as t → ∞.
In the second step, we plug the expansion (12) into Proposition 5 and order the outcome according to powers of s − ζ * . This will complete the proof of spectrally positive case (SP).
Similarly, when (SN) holds then as s ↓ ζ * . Then, using Proposition 5 in the same way as before gives the required assertion after some simple manipulations.

Examples
In this section, we illustrate our theory through a few examples.