The correlation function of a queue with Levy and Markov additive input

TL;DR

This paper investigates the correlation function of queue workloads driven by Levy and Markov additive processes, proving conjectures, providing counterexamples, and analyzing decay rates and tail behaviors in different process settings.

Contribution

It extends the analysis of correlation functions to spectrally two-sided Levy processes and Markov additive processes, including proofs, counterexamples, and decay rate identification.

Findings

Correlation function properties extend to spectrally two-sided Levy processes with reflection.

Counterexamples show correlation can be negative, decreasing, and concave for Markov additive processes.

Decay rate analysis reveals tail behavior differences between Levy and Markov additive inputs.

Abstract

Let $(Q_t)$ be a stationary workload process, and $r(t)$ the correlation coefficient of $Q_0$ and $Q_t$. In a series of previous papers (i) the transform of $r(\cdot)$ has been derived for the case that the driving process is spectrally-positive (sp) or spectrally-negative (sn) Levy, (ii) it has been shown that for sp-Levy and sn-Levy input $r(\cdot)$ is positive, decreasing, and convex, (iii) in case the driving Levy process is light-tailed (a condition that is automatically fulfilled in the sn case), the decay of the decay rate agrees with that of the tail of the busy period distribution. In the present paper we first prove the conjecture that property (ii) carries over to spectrally two-sided Levy processes; we do so for the case the Levy process is reflected at 0, and the case it is reflected at 0 and $K > 0$. Then we focus on queues fed by Markov additive processes (maps). We start by the establishing the counterpart of (i) for sp- and sn-maps. Then we refute property (ii) for maps: we construct examples in which the correlation coefficient can be (locally) negative, decreasing, and concave. Finally, in relation to (iii), we point out how to identify the decay rate of $r(\cdot)$ in the light-tailed map case, thus showing that the tail behavior of $r(\cdot)$ does not necessarily match that of the busy-period tail; singularities related to the transition rate matrix of the background Markov chain turn out to play a crucial role here.

Equations

$$r(t):=\frac{\mathbb{C}\mathrm{ov}(Q_0,Q_t)}{\sqrt{\mathbb{V}\mathrm{ar}{Q}_0\cdot \mathbb{V}\mathrm{ar}{Q}_t}},$$

$$Q_t:=X_t+\max \left\{Q_0,-\underset{0⩽u⩽t}{\inf }{X}_u\right\}.$$

$$Q_t=\underset{s⩽t}{\sup }(X_t-X_s);$$

$$X_{s,t}:=X_t-X_s,L_{s,t}:=-\underset{s⩽u⩽t}{\inf }{X}_{s,t}.$$

$$Q_t=X_{s,t}+\max \{Q_s,L_{s,t}\}=X_{s,t}+Q_s+\max \{0,L_{s,t}-Q_s\};$$

$$Q_u-Q_t=$$

$$=$$

$$r(t):=\frac{\mathbb{C}\mathrm{ov}(Q_0,Q_t)}{\sqrt{\mathbb{V}\mathrm{ar}{Q}_0\cdot \mathbb{V}\mathrm{ar}{Q}_t}}=\frac{\mathbb{E}{Q}_0{Q}_t-(\mathbb{E}{Q}_0{)}^2}{\mathbb{V}\mathrm{ar}{Q}_0}.$$

$$\mathbb{E}(Af(A))⩾\mathbb{E}(Bf(A)),\mathbb{E}(Ag(A))⩽\mathbb{E}(Bg(A)).$$

$$\mathbb{C}\mathrm{ov}(Q_0,Q_t)$$

$$\displaystyle\>={\mathbb{E}}\big{(}{\mathbb{E}}((Q_{0}-\bar{Q})(X_{0,t}+\max\{Q_{0},L_{0,t}\}))\,|\,{\mathscr{G}}_{0}\big{)},$$

$$f(a):=X_{0,t}+\max \{a,L_{0,t}\}.$$

$${\mathbb{C}}{\rm ov}(Q_{0},Q_{t})={\mathbb{E}}\big{(}{\mathbb{E}}((Q_{0}-\bar{Q})(X_{0,t}+\max\{Q_{0},L_{0,t}\}))\,|\,{\mathscr{G}}_{0}\big{)}\geqslant 0,$$

$$f(a):=X_{s,t}+\max \{a,L_{s,t}\}.$$

$${\mathbb{E}}\big{(}{\mathbb{E}}(Q_{0}(X_{s,t}+\max\{Q_{s},L_{s,t}\}))|\,{\mathscr{G}}_{s}\big{)}\leqslant{\mathbb{E}}\big{(}{\mathbb{E}}(Q_{s}(X_{s,t}+\max\{Q_{s},L_{s,t}\}))|\,{\mathscr{G}}_{s}\big{)}.$$

$$\mathbb{E}(Q_0{Q}_t)$$

$$\displaystyle\>\leqslant{\mathbb{E}}\big{(}{\mathbb{E}}(Q_{s}(X_{s,t}+\max\{Q_{s},L_{s,t}\}))|\,{\mathscr{G}}_{s}\big{)}={\mathbb{E}}(Q_{s}Q_{t})={\mathbb{E}}(Q_{0}Q_{t-s}),$$

$$\mathbb{E}(Q_0(Q_u-Q_t)⩾\mathbb{E}(Q_s(Q_u-Q_t))=\mathbb{E}(Q_0(Q_{u-s}-Q_{t-s})),$$

$$f(a):=X_{t,u}+\max \{0,L_{t,u}-X_{s,t}-\max \{a,L_{s,t}\}\}.$$

$$\mathbb{E}(Q_0(Q_u-Q_t))$$

$$\displaystyle\>={\mathbb{E}}\big{(}{\mathbb{E}}(Q_{0}(X_{t,u}+\max\{0,L_{t,u}-X_{s,t}-\max\{Q_{s},L_{s,t}\}\}))|\,{\mathscr{G}}_{s}\big{)}$$

$$\displaystyle\>\geqslant{\mathbb{E}}\big{(}{\mathbb{E}}(Q_{s}(X_{t,u}+\max\{0,L_{t,u}-X_{s,t}-\max\{Q_{s},L_{s,t}\}\}))|\,{\mathscr{G}}_{s}\big{)}$$

$$=\mathbb{E}(Q_s(Q_u-Q_t)).$$

$$W_{n+1}=\max \{W_n+U_n,0\}.$$

$$W_{n+1}=\max \left\{\underset{i=-\infty ,\dots ,n}{\max }\sum _{j=i}^n{U}_j,0\right\},$$

$$Q_t=Q_0+X_t+R_t^{-}-R_t^{+},$$

$${\int }_0^{\infty }{X}_t\mathrm{d}{R}_t^{-}={\int }_0^{\infty }(K-X_t)\mathrm{d}{R}_t^{+}=0.$$

$$Q_t=Q_s+X_{s,t}+R_{s,t}^{-}-R_{s,t}^{+}.$$

$$Q_u-Q_t=X_{t,u}+R_{t,u}^{-}-R_{t,u}^{+}.$$

$${\Gamma }^{-}[Y{]}_t:=Y_t+\max \left\{0,-\underset{0⩽u⩽t}{\inf }{Y}_u\right\};$$

$${\Gamma }^{+}[Y{]}_t:=Y_t-\underset{0⩽s⩽t}{\sup }\left(\max \{Y_s-K,0\},\underset{s⩽u⩽t}{\inf }{Y}_u\right).$$

$$Q_t=({\Gamma }^{+}\circ {\Gamma }^{-})[Q_0+X{]}_t={\Gamma }^{+}[{\Gamma }^{-}[Q_0+X]{]}_t.$$

$${\Gamma }^{-}[Y^{(2)}{]}_t-Y_t^{(2)}=$$

Full text

In the present paper we first prove the conjecture that property (ii) carries over to spectrally two-sided Lévy processes; we do so for the case the Lévy process is reflected at 0, and the case it is reflected at [math] and $K>0$. Then we focus on queues fed by Markov additive processes (maps). We start by the establishing the counterpart of (i) for sp- and sn-maps. Then we refute property (ii) for maps: we construct examples in which the correlation coefficient can be (locally) negative, decreasing, and concave. Finally, in relation to (iii), we point out how to identify the decay rate of $r(\cdot)$ in the light-tailed map case, thus showing that the tail behavior of $r(\cdot)$ does not necessarily match that of the busy-period tail; singularities related to the transition rate matrix of the background Markov chain turn out to play a crucial role here.

Keywords. Lévy processes $\circ$ reflection $\circ$ workload $\circ$ Markov additive processes

Affiliations. Wouter Berkelmans and Agata Cichocka are with Department of Mathematics, Vrije Universiteit Amsterdam, De Boelelaan 1081a, 1081 HV Amsterdam, the Netherlands. Michel Mandjes is with Korteweg-de Vries Institute for Mathematics, University of Amsterdam, Science Park 904, 1098 XH Amsterdam, the Netherlands. He is also with Eurandom, Eindhoven University of Technology, Eindhoven, the Netherlands, and Amsterdam Business School, Faculty of Economics and Business, University of Amsterdam, Amsterdam, the Netherlands. His research is partly funded by NWO Gravitation project Networks, grant number 024.002.003.

Acknowledgments. The third author wishes to thank Lars Nørvang Andersen (Aarhus) and François Baccelli (Austin) for helpful interactions regarding the spectrally two-sided Lévy case, Jevgenijs Ivanovs (Aarhus) and Zbigniew Palmowski (Wrocław) for providing useful references, and Peter Glynn (Stanford) for discussions regarding the tail behavior examples in Section 6.

1. Introduction

Consider a queueing resource with workload process $(Q_{t})_{t\in{\mathbb{R}}}$. Where much effort has been spent on characterizing transient and stationary distributional properties pertaining to this workload process, this paper focuses on the workload correlation

[TABLE]

assuming the workload is in stationarity at time [math]. Insight into $r(\cdot)$ is useful in various ways, as it tells us from what time horizon $\tau$ on the random variables $Q_{0}$ and $Q_{\tau}$ can safely be assumed independent (in that $r(\tau)$ drops below some predefined level $\varepsilon$).

In a series of articles, the function $r(\cdot)$ has been analyzed for the important class of Lévy-driven queues. These can be seen as queues fed by Lévy input [6], and are also known as reflected Lévy processes [15]. Important subclasses of Lévy processes are spectrally positive (in this paper abbreviated to sp) Lévy processes, with jumps only in the upward direction, and spectrally negative (sn) Lévy processes, with jumps only in the downward direction. The results on the workload correlation function $r(\cdot)$ that were derived so far can be divided into three types: (i) expressions for the transform $\widehat{r}(\cdot)$ of $r(\cdot)$, (ii) structural properties concerning the shape of $r(\cdot)$, and (iii) characterizations of $r(t)$ for $t$ large.

We proceed by giving a brief account of the literature. Ott [18] considered a subclass of the sp Lévy input case, viz. queues with compound Poisson input (i.e., queues of the M/G/1 type). He succeeded in finding $\widehat{r}(\cdot)$ in closed-form. In addition, he used this transform to prove, relying on the machinery of completely-monotone functions, that $r(\cdot)$ is non-negative, non-increasing, and convex. These results were extended by Es-Saghouani and Mandjes [8] to the full class of sp Lévy inputs (thus also covering Lévy processes with a Brownian component, as well as components with so-called ‘small jumps’ in the upward direction). In addition, they identified the tail behavior of $r(\cdot)$; in particular, when the sp Lévy input is light-tailed, it was shown that its (exponential) decay rate matches that of the tail of the busy-period distribution. Glynn and Mandjes [9] proved the counterparts of the results of [8] for the case of sn Lévy input, and in addition devised a coupling-based simulation technique to efficiently estimate $r(t)$. The findings of [8, 9] show that the derivation of the Laplace transform $\widehat{r}(\cdot)$ requires the availability of an explicit expression for the Laplace transform of the stationary workload, as well as for the Laplace transform of the transient workload (i.e., after an exponentially distributed interval $T$, conditional on the workload at time [math]).

Extensions in various directions could be thought of. Two of them are:

• $\circ$

In the first place, one could consider queues with general Lévy input (queues fed by a spectrally two-sided Lévy process, that is). A long outstanding conjecture is that the fact that $r(\cdot)$ is non-negative, non-increasing, and convex, which was proven for sp and sn Lévy input, carries over to the spectrally two-sided Lévy setting. There is a similar conjecture for the correlation of the waiting times in the GI/G/1 queue.

It is clear that a proof for the spectrally two-sided case requires a different approach than the one used in the one-sided cases. The major complication is that, unlike in the sp and sn cases, no explicit results for the Laplace transforms of the stationary and transient workload are available — results are typically in terms of ‘Wiener-Hopf expressions’; see e.g. [6, Thm. 4.3]. As a result, identifying the transform $\widehat{r}(\cdot)$ seems to be out of reach. As in the sp and sn cases this Laplace transform was crucial in establishing that $r(\cdot)$ is non-negative, non-increasing, and convex (using the concept of completely-monotone functions), one may wonder whether these structural properties can be proven without knowing $\widehat{r}(\cdot).$

In addition, one may wonder whether the property that $r(\cdot)$ is non-negative, non-increasing, and convex remains valid if the Lévy process is doubly reflected, i.e., reflected both from below at level [math] and from above at level $K>0$ (corresponding with the workload process of a queue with finite buffer capacity $K$).

• $\circ$

Secondly, one could consider queues fed by spectrally one-sided Markov additive processes (maps). These processes, dating back to [5, 17], can be seen as the Markov-modulated of Lévy processes. Informally, a map can be seen as a Lévy process, but with its local behavior depending on the state of an external finite-state Markov process, usually referred to as the background process. Only for very specific cases, partial results on the workload correlation function are known; see e.g. the findings in [12].

We now describe the contributions of our paper. In the first place, we prove the conjecture that in the case of two-sided Lévy input the workload correlation function is non-negative, non-increasing, and convex. We do so without relying on explicit expressions for the transform $\widehat{r}(\cdot)$; instead we provide a clean and insightful proof that uses first principles only. The argument carries over to the correlation of waiting times in the GI/G/1 queue. We also show that in case we add reflection at $K>0$, the property that $r(\cdot)$ is non-negative, non-increasing, and convex remains true.

We then consider map-driven queues, focusing on the cases with sp-map and sn-map input. For both classes we succeed in deriving an explicit expression for the Laplace transform $\widehat{r}(\cdot)$, making intensive use of earlier fluctuation-theoretic results for spectrally one-sided maps. One may wonder whether also for map input the workload correlation function is non-negative, non-increasing, and convex; by a series of insightful examples, we show that this is not the case. We also study the tail behavior of $r(\cdot)$, with the main conclusion that in the light-tailed setting the decay rate does not match with the busy-period decay rate (unlike in the Lévy case). More specifically, as it turns out, the rate of decay of $r(\cdot)$ follows by comparing the busy-period decay rate with singularities related to the transition rate matrix of the background process.

2. Lévy input: structural properties of correlation function

In this section we consider the workload process pertaining to a queue with Lévy input. In the first subsection we consider the case of reflection from below at level 0, whereas the second studies its finite-buffer counterpart in which reflection from above at level $K>0$ is added.

2.1. Single-sided reflection

Let $(X_{t})_{t\in{\mathbb{R}}}$ be a single-dimensional Lévy process [15] such that ${\mathbb{E}}\,X_{1}<0.$ Then we define $(Q_{t})_{t\in{\mathbb{R}}}$ as the associated workload process (or: the process $(X_{t})_{t\in{\mathbb{R}}}$ reflected at 0) by

[TABLE]

We assume the workload process is in stationarity, meaning that (under our stability condition ${\mathbb{E}}\,X_{1}<0$) we can alternatively write

[TABLE]

due to the stationarity, we in particular have that, for any $t$, $Q_{t}$ is distributed as $Q_{0}.$ For future reference, we also introduce for $s\leqslant t$,

[TABLE]

In addition, we let ${\mathscr{G}}_{s}$ be $\sigma\{(X_{s,t})_{t\geqslant s}\}$, i.e., the sigma-algebra generated by the increments of $(X_{t})_{t\in{\mathbb{R}}}$ after time $s$ (relative to $X_{s}$, that is). In the sequel, we also work with $\bar{Q}$, which has the stationary workload distribution (and is therefore distributed as $Q_{0}$), but is independent of the process $(X_{t})_{t\in{\mathbb{R}}}$ (and hence also of the process $(Q_{t})_{t\in{\mathbb{R}}}$).

For $s\leqslant t$,

[TABLE]

observe that (i) $Q_{t}$ is non-decreasing in $Q_{s}$, and (ii) $L_{s,t}$ is independent of $Q_{s}$. This relation implies that, for $s\leqslant t\leqslant u$,

[TABLE]

observe that (i) $Q_{u}-Q_{t}$ is non-increasing in $Q_{s}$, and (ii) $L_{s,t}$, $X_{s,t}$, $L_{t,u}$, and $X_{t,u}$ are independent of $Q_{s}$.

We throughout assume that ${\mathbb{V}}{\rm ar}\,Q_{0}<\infty$, implying that ${\mathbb{C}}{\rm ov}(Q_{0},Q_{t})<\infty$. The object of study is

[TABLE]

Our main objective is to prove, using first principles, that $r(\cdot)$ is non-negative, non-increasing, and convex (which was shown for sp-Lévy and sn-Lévy input in [8, 9]).

The following lemma will play a pivotal role in our analysis; its proof is provided in the appendix.

Lemma 1**.**

Let $A$ and $B$ be identically distributed non-negative random variables. Let $f(\cdot)$ be non-decreasing, and $g(\cdot)$ be non-increasing. Then,

[TABLE]

The main result of this section is the following.

Theorem 1**.**

In the system with the Lévy process $(X_{t})_{t\in{\mathbb{R}}}$ being reflected from below at 0, $r(\cdot)$ is non-negative, non-increasing, and convex.

Proof: (A) We first show that ${\mathbb{C}}{\rm ov}(Q_{0},Q_{t})$ is non-negative. To this end, observe that

[TABLE]

where the last equality is due to (3). Now apply Lemma 1 with $A:=Q_{0}$, $B:=\bar{Q}$, and

[TABLE]

To verify that the conditions of the lemma are met, first note that $f(\cdot)$ is non-decreasing. In the second place, $Q_{0}$ is independent of ${\mathscr{G}}_{0}$ (i.e., the Lévy process after time 0), and has therefore the same distribution as $\bar{Q}$. The lemma thus yields

[TABLE]

as desired.

(B) We now prove that ${\mathbb{C}}{\rm ov}(Q_{0},Q_{t})$ is non-increasing, for which it is sufficient that ${\mathbb{E}}(Q_{0}Q_{t})$ is non-increasing. Again we use Lemma 1, but now with $A:=Q_{s}$, $B:=Q_{0}$ (with $0\leqslant s\leqslant t$), and

[TABLE]

Evidently, $f(\cdot)$ is non-decreasing and conditional on ${\mathscr{G}}_{s}$ the random variables $Q_{0}$ and $Q_{s}$ have the same (conditional) distributions. The lemma thus yields

[TABLE]

We obtain

[TABLE]

with the last equality being valid due to stationarity.

(C) In the third part of the proof, we show that ${\mathbb{C}}{\rm ov}(Q_{0},Q_{t})$ is convex. We do so by proving the sufficient property, for $s\leqslant t\leqslant u$,

[TABLE]

where the last equality is due to the stationarity; combined with the monotonicity of ${\mathbb{E}}(Q_{0}Q_{t})$ this implies convexity. In this case we define

[TABLE]

Again we wish to apply Lemma 1 with $A:=Q_{s}$ and $B:=Q_{0}$. Now $f(\cdot)$ is non-increasing. In addition, conditional on ${\mathscr{G}}_{s}$, we again have that $Q_{0}$ and $Q_{s}$ are identically distributed. It thus follows that, using (3),

[TABLE]

This completes the proof. $\Box$

Remark 1**.**

Waiting times in the GI/G/1 queue are known to satisfy the Lindley recursion. Indeed, with $W_{n}$ the waiting time of the $n$-th customer, $A_{n}$ the time between the arrivals of the $n$-th and $(n+1)$-st arrivals, $S_{n}$ the service time of the $n$-th customer, and $U_{n}:=S_{n}-A_{n}$,

[TABLE]

Iterating this relation, we obtain that, under the stability condition ${\mathbb{E}}(S_{n})<{\mathbb{E}}(A_{n})$,

[TABLE]

which is the discrete-time counterpart of (2). Mimicking all steps that we used in the proof of Thm. 1, it follows that the correlation function of $(W_{n})_{n\in{\mathbb{Z}}}$ is non-negative, non-increasing and convex as well. $\hfill\diamond$

2.2. Double-sided reflection

The objective of this section is to show that Thm. 1 carries over to finite-buffer systems. We use an argumentation very similar to the one used for the case of single-sided reflection.

With $K>0$ denoting the system’s buffer capacity, one can write [2, Section 1] the workload under double-sided reflection as

[TABLE]

with the ‘regulators’ $(R_{t}^{-})_{t\geqslant 0}$ and $(R_{t}^{+})_{t\geqslant 0}$ non-decreasing processes satisfying

[TABLE]

In self-evident notation we thus have, for $s\leqslant t$, with $R_{s,t}^{-}$ representing the local time in 0 in $[s,t]$, and $R_{s,t}^{+}$ the local time in $K$ in $[s,t]$,

[TABLE]

Therefore, for $s\leqslant t\leqslant u$,

[TABLE]

Like in the system with single-sided reflection (presented in the previous subsection), $Q_{t}$ is non-decreasing in $Q_{s}$, and in addition $Q_{u}-Q_{t}$ is non-increasing in $Q_{s}$. The former statement follows by comparing two instances of the workload process, one being at $q_{-}\in[0,K)$ at time $s$, and the other at $q_{+}\in(q_{-},K]$; then the claim can be concluded by observing that, when using the same driving Lévy process after time $s$, the trajectories of both workload processes do not cross. The latter statement is effectively a consequence of $R_{t,u}^{-}$ (recalling that it represents the local time in 0 in $[t,u]$) being non-increasing in $Q_{s}$ and $R_{t,u}^{+}$ (i.e., the local time in $K$ in $[t,u]$) being non-decreasing in $Q_{s}$; formally, this property can be established by using the following lemma.

Lemma 2**.**

Let $0\leqslant t\leqslant u$. Then $Q_{u}-Q_{t}$ is non-increasing in $Q_{0}.$

Proof: We first define, for a given initial level $Q_{0}\in[0,K]$ and a path $X$ of the Lévy process, a reflection from below at level 0. To this end, we introduce, as in (1), the functional $\Gamma^{-}[\,\cdot\,]$, given by

[TABLE]

The single-sided reflection (from below, at level [math]) we obtain by imposing $\Gamma^{-}[\,\cdot\,]$ on $Y$, where $Y_{s}=Q_{0}+X_{s}$, in line with (1).

Now introduce the functional $\Gamma^{+}[\,\cdot\,]$ by

[TABLE]

Then, according to [14], for our model with double reflection,

[TABLE]

In other words, the process $(Q_{t})_{t\geqslant 0}$ is obtained by imposing $\Gamma^{+}\circ\Gamma^{-}$ on the process $Q_{0}+(X_{t})_{t\geqslant 0}.$

Now consider paths $Y^{(1)}$ and $Y^{(2)}$ such that $Y^{(1)}_{t}\leqslant Y^{(2)}_{t}$ for all $t\geqslant 0.$ Then, considering the functional $\Gamma^{-}[\,\cdot\,]$,

[TABLE]

implying that

[TABLE]

A similar reasoning applies for the functional $\Gamma^{+}[\,\cdot\,]$: for trajectories $Y^{(1)}$ and $Y^{(2)}$ such that $Y^{(1)}_{t}\leqslant Y^{(2)}_{t}$ for all $t\geqslant 0$,

[TABLE]

such that

[TABLE]

The next step is to combine the bounds (6) and (7). Our goal is to show that $Q_{t}-Q_{0}$ is non-increasing in $Q_{0}$ (conditional on the path of the driving Lévy process, say $X$). To this end, consider again two instances of the workload process, one being at $Q_{0}:=q_{-}\in[0,K]$ at time $s$, and the other at $Q_{0}:=q_{+}\in[q_{-},K]$; we thus have that the workloads (at time $t\geqslant 0$) corresponding with both instances are given by

[TABLE]

First recall that $\Gamma^{-}[q_{-}+X]_{t}\leqslant\Gamma^{-}[q_{+}+X]_{t}$. Hence applying (7) yields, for any $t\geqslant 0$,

[TABLE]

Now using (6), the quantity in the right-hand side of the previous display is no larger than $q_{+}-q_{-}.$ Combining the above, we conclude

[TABLE]

recalling that $q_{+}\geqslant q_{-}$, this is equivalent to $Q_{t}-Q_{0}$ being non-increasing in $Q_{0}$ (for a given path $X$, that is). By shifting time, this entails that, for $0\leqslant t\leqslant u$, $Q_{u}-Q_{t}$ being non-increasing in $Q_{t}$. Combining this with the fact that $Q_{t}$ is non-decreasing in $Q_{0}$, we have proven the claim. $\Box$

Note that the above lemma immediately implies that (for $s\leqslant t\leqslant u$) $Q_{u}-Q_{t}$ is non-increasing in $Q_{s}$. This means that from this point on, the argumentation we have developed for the infinite-buffer case can be applied in the finite-buffer context as well. We conclude that Thm. 1 extends to the finite-buffer case.

Theorem 2**.**

In the system with the Lévy process $(X_{t})_{t\in{\mathbb{R}}}$ being reflected from below at 0 and from above at $K>0$, $r(\cdot)$ is non-negative, non-increasing, and convex.

In the setting above there is reflection from below at [math] and from above at $K>0$, but evidently this strip can be shifted. In this way we can cover the setting of reflection from below at $K_{1}$ and from above at $K_{2}$ (with $K_{2}>K_{2}$).

3. Markov additive input: model and preliminaries

In the following sections, our objective is to analyze the workload correlation function of a queue fed by spectrally one-sided Markov additive input. In this section, we first provide a formal definition of maps, and introduce the notation needed. Then we present a number of useful preliminary results.

3.1. Model and notation

As mentioned, we study the workload of a queue fed by a map. A map is a bivariate Markovian process $(X_{t},J_{t})$, defined as follows.

• $\circ$

Let $(J_{t})_{t}$ be an irreducible continuous-time Markov chain with finite state space $E=\{1,\ldots,d\}$. Define by $\mathscr{Q}:=(q_{ij})_{i,j=1}^{d}$ the $(d\times d)$ transition rate matrix of $(J_{t})_{t}$. Let ${\boldsymbol{\pi}}$ be the (unique) stationary distribution; recall that ${\boldsymbol{\pi}}^{\rm T}\mathscr{Q}={\boldsymbol{0}}^{\rm T}$. In addition, $q_{ij}\geqslant 0$ if $i\not=j$ and $\mathscr{Q}{\boldsymbol{1}}={\boldsymbol{0}}$ (i.e., the row sums are 0). Define $\widehat{q}_{i}:=-q_{ii}\geqslant 0.$

• $\circ$

For each state $i$ that $J_{t}$ can attain, let the process $(X^{(i)}_{t})_{t}$ be a Lévy process. We either assume that either all of these are sp or that all of them are sn. In the former case they have Laplace exponents

[TABLE]

whereas in the latter case they have cumulant generating functions

[TABLE]

for $i=1,\ldots,d$

• $\circ$

Letting $T_{n}$ and $T_{n+1}$ be two successive transition epochs of $J_{t}$, and given that $J_{t}$ jumps from state $i$ to state $j$ at $t=T_{n}$, we define the additive process $X_{t}$ for $t$ in the time interval $[T_{n},T_{n+1})$ through

[TABLE]

where the $(U_{ij}^{n})_{n}$ constitute a sequence of i.i.d.  random variables (each of which is distributed as a generic random variable $U_{ij}$), describing the jumps at transition epochs. In the sp case we work with the Laplace-Stieltjes transform (with $i\in\{1,\ldots,d\}$)

[TABLE]

where $U_{ij}\geqslant 0$ almost surely. Likewise, in the sn case we work with the moment generating functions

[TABLE]

where $U_{ij}\leqslant 0$ almost surely.

Define, in the sp case,

[TABLE]

and in the sn case,

[TABLE]

Then, as pointed out in e.g. [6, Prop. 11.3], the matrices with entries $\sigma_{ij}(\alpha,t)$ and $\Sigma_{ij}(\beta,t)$ can be written in terms of matrix exponentials. More concretely,

[TABLE]

for a matrix $s(\alpha)$ with entries $s_{ij}(\alpha)$, and a matrix $S(\beta)$ with entries $S_{ij}(\beta)$.

Throughout we write ${\mathbb{P}}_{i}(A):={\mathbb{P}}(A\,|\,J_{0}=i)$ and ${\mathbb{E}}_{i}(Y):={\mathbb{E}}(Y\,|\,J_{0}=i)$ for events $A$ and random variables $Y$. Also, $\langle{\boldsymbol{a}},{\boldsymbol{b}}\rangle:={\boldsymbol{a}}^{\rm T}{\boldsymbol{b}}$ for ${\boldsymbol{a}},{\boldsymbol{b}}\in{\mathbb{R}}^{d}.$

The queueing processes are assumed to be stable, meaning that we impose, for the sp and sn case, respectively,

[TABLE]

3.2. Objectives and preliminary results

We now define the workload covariance function through its transform with respect to time, as follows. Assume the queueing process is in stationarity at time [math]. The goal is to study, for $\vartheta>0$,

[TABLE]

where $T\equiv T(\vartheta)$ is exponentially distributed with mean $\vartheta^{-1}$, independently of the driving map. Withe $r(t)$ the correlation between $Q_{0}$ and $Q_{t}$ and $\widehat{r}(\cdot)$ the Laplace transform of $r(\cdot)$, we thus have the evident relation

[TABLE]

In the next two sections we focus on deriving an expression for $\gamma(\vartheta)$ (for the sp and sn cases, respectively), whereas the examples in Section 6 provide insight in the shape of $r(\cdot)$ and its asymptotics for $t$ large.

If the underlying queueing system would have been Markovian, this expression could have been rewritten as

[TABLE]

with $p(x)\,{\rm d}x:=\,{\mathbb{P}}(Q_{0}\in{\rm d}x)$ and $m(x):={\mathbb{E}}(Q_{T}\,|\,Q_{0}=x)$. In our case, however $(Q_{t})_{t\geqslant 0}$ (alone) is not Markovian, but $(Q_{t},J_{t})_{t\geqslant 0}$ jointly are. Define $p_{i}(x)={\mathbb{P}}(Q_{0}\in{\rm d}x,J_{0}=i)$, and $m_{i}(x):={\mathbb{E}}(Q_{T}\,|\,Q_{0}=x,J_{0}=i).$ In evident vector notation, we obtain that in our setting with map input

[TABLE]

The following result will be used several times. Define $\tau(x):=\inf\{t\geqslant 0:X_{t}\leqslant-x\}$ for $x\geqslant 0$.

Lemma 3**.**

For all $x\geqslant 0$,

[TABLE]

Proof: The quantity ${\boldsymbol{m}}(x)$ can be determined by investigating the effect of slightly perturbing the initial workload level $x$. More precisely, we compare $m_{i}(x+\delta)$ and $m_{i}(x)$ for $\delta$ small. Observe that if before $T$ the queues have drained, the queueing processes have coupled, and as a consequence the workloads are the same; if, on the contrary, the busy period has not ended before $T$, the workloads still differ by $\delta.$ This entails that, as $\delta\downarrow 0$, for an $f\in(0,1)$,

[TABLE]

with $\underline{X}_{t}$ denoting the minimum of $X_{s}$ over $s\in[0,t]$. Dividing both sides of the previous display by $\delta$ and letting $\delta\downarrow 0$, we thus obtain the differential equation

[TABLE]

By conditioning on the value of $\tau(x)$,

[TABLE]

We thus find that, for constants $C$ and $K$,

[TABLE]

The constant $K$ can be determined as follows. As $x\to\infty$, we have that $m_{i}(x)-x\downarrow{\mathbb{E}}_{i}X_{T}$ (the underlying reason being that as the initial level $x$ increases to $\infty$, the probability of reflection before $T$ vanishes). It thus follows that $K={\mathbb{E}}_{i}X_{T}.$ $\hfill\Box$

Later in this paper we use a result for the transform of ${\mathbb{E}}_{i}\,{\rm e}^{-\vartheta\,\tau(y)}$ (with respect to $y$, that is). To this end, recall that we have observed above that

[TABLE]

Recall that ${\mathbb{P}}_{i}(\tau(y)<t)={\mathbb{P}}_{i}(\underline{X}_{t}<-y).$ Upon combining these two facts,

[TABLE]

Now using integration by parts, we end up with the following result. We remark that its validity is not restricted to spectrally one-sided Lévy processes.

Lemma 4**.**

For $\eta\geqslant 0$,

[TABLE]

4. Markov additive input: spectrally positive case

The main objective of this section is to identify the transform $\gamma(\cdot)$ for the case of sp-map input. The main idea is to use representation (9). To do so, we first collect a number of useful expressions.

In the first place, in this sp-map case the Laplace transform of $\kappa(\alpha)$ of $Q_{0}$ is given by

[TABLE]

see for instance [6, Eqn. (11.1)] and [7, Corollary 3.1]. The row vector ${\boldsymbol{\ell}}$ can be evaluated using [7, Lemmas 3.1 and 3.2]. This result can be seen as the extension of the classical Pollaczek-Khinchine formula [6, Thm. 3.2] for spectrally positive Lévy input [19] to the spectrally positive map case. It follows that ${\mathbb{E}}\,Q_{0}=-\kappa^{\prime}(0).$ In addition, ${\mathbb{V}}{\rm ar}\,Q_{0}=\kappa^{\prime\prime}(0)-(\kappa^{\prime}(0))^{2}.$

A somewhat more precise form of the above result for $\kappa(\alpha)$ concerns the distribution of $Q_{0}$ jointly with the background state $J_{0}$:

[TABLE]

Bearing in mind Eqn. (9) and Lemma 3, the next question is: how to evaluate the quantities $\xi_{i}(y,\vartheta):={\mathbb{E}}_{i}\,{\rm e}^{-\vartheta\,\tau(y)}$ and $\nu_{i}\equiv\nu_{i}(\vartheta):={\mathbb{E}}_{i}X_{T}$; if we could evaluate these, we have an expression for $\gamma(\vartheta).$

• $\circ$

We can determine $\nu_{i}\equiv\nu_{i}(\vartheta)$ as follows. As we have seen before

[TABLE]

so that

[TABLE]

There is, however, an easier way to compute $\nu_{i}$. Observe that, using the memoryless property,

[TABLE]

this (diagonally dominant) linear system has a unique solution, and can be solved by standard methods.

• $\circ$

We now point out how to evaluate ${\mathbb{E}}_{i}\,{\rm e}^{-\vartheta\,\tau(y)}$. The following standard reasoning can be used. In the first place, observe that $(\tau(y),J_{\tau(y)})_{y\geqslant 0}$ is a map; $\tau_{x}$ is increasing (and therefore sp), and (due to the drift condition) $\tau_{x}\uparrow\infty$ as $x\to\infty.$ Now define the matrix $A(\vartheta)$ whose $(i,j)$-th component is given by

[TABLE]

Because $(\tau(x),J_{\tau(x)})_{x\geqslant 0}$ is an sp-map, there is matrix-valued function $\Psi(\cdot)$ such that

[TABLE]

for characterizations of $\Psi(\cdot)$, see e.g. [6, 16]. We thus conclude that there are matrices $A(\vartheta)$ and $\Psi(\vartheta)$ such that

[TABLE]

We refer for more background to e.g. [16], and, for the case of phase-type jumps, [4]. Depending on the linear-algebraic properties of the matrix $\Psi(\vartheta)$, this expression can be further evaluated. For instance if all eigenvalues of $\Psi(\vartheta)$ are distinct (say $\psi_{1}(\vartheta)$ up to $\psi_{d}(\vartheta)$), we can write

[TABLE]

for appropriate $c_{ij}(\vartheta)$. For the moment we assume that this expansion applies; if it does not apply (when eigenvalues have multiplicities larger than 1), then straightforward adaptations can be made, as explained in Remark 3.

We thus obtain that

[TABLE]

We are now in a position to evaluate $\gamma(\vartheta).$ We have

[TABLE]

Here $\pi_{i}\,{\mathbb{E}}_{i}Q_{0}=-\kappa_{i}^{\prime}(0),$ and, as was mentioned earlier, ${\mathbb{V}}{\rm ar}\,Q_{0}=\kappa^{\prime\prime}(0)-(\kappa^{\prime}(0))^{2}.$ We have thus found the following representation for $\gamma(\vartheta)$.

Theorem 3**.**

If the expansion $(\ref{expa})$ applies, then, for any $\vartheta\geqslant 0$,

[TABLE]

Remark 2**.**

The above theorem can be considered as the map counterpart of the corresponding result for a queue driven by a spectrally positive Lévy process [8, Thm. 2.1]. Observe that, so as to translate the results into one another, the expression for $\gamma(\vartheta)$ has to be divided by $\vartheta$ as well as by ${\mathbb{V}}{\rm ar}\,Q_{0}$, cf. (8).* $\diamond$*

Remark 3**.**

Condition (10) is imposed for convenience. If there are mixed exponential-polynomial terms in the expression for $\xi_{i}(y,\vartheta)$ (due to eigenvalues with multiplicities higher than 1), one could proceed as follows. Consider for instance the case that the eigenvalue $\psi_{j}(\vartheta)$ having multiplicity two leads to a term proportional to $y\,{\rm e}^{-\psi_{j}(\vartheta)\,y}.$ Due to the identity

[TABLE]

we find that for this $j$ (and all $i\in\{1,\ldots,d\}$) the $-\kappa^{\prime}_{i}(\psi_{j}(\vartheta))/\psi_{j}(\vartheta)$ appearing in Thm.  3 has to be replaced by

[TABLE]

Analogously, if a term $y^{k}\,{\rm e}^{-\psi_{j}(\vartheta)\,y}$ appears for some $k\in\{0,\ldots,d-1\}$, then the $(k+1)$-st derivative of $\kappa(\cdot)$ will end up in the expression for $\gamma(\vartheta).$

• $\diamond$*

5. Markov additive input: spectrally negative case

In this section we provide an expression for $\gamma(\vartheta)$ for the case of sn-map input. As it turns out, now the quantities $\xi_{i}(y,\vartheta)={\mathbb{E}}_{i}\,{\rm e}^{-\vartheta\,\tau(y)}$ can not be given (unlike in the sp-map case), but their transforms with respect to $y$ can be evaluated. Combining this with the fact that $Q_{0}$ is of phase type facilitates the computation of $\gamma(\vartheta).$ Below we detail this reasoning.

Lemma 5**.**

In the sn-map case, $Q_{0}$ (jointly with the state $J_{0}$ being $i$) has a phase-type distribution. In particular, we can write $Q_{0}$ (jointly with $J_{0}=i$) is given by the initial distribution ${\boldsymbol{t}}_{i}$ and transition rate matrix $E$.

Proof: This property follows directly using (for instance) the results from [6]. The reasoning consists of the following steps. (i) From [6, Prop. 4.2] it follows how to translate (by using time-reversal) the distribution of $Q_{0}$ (jointly with $J_{0}=i$) into the distribution of the all-time maximum of the non-reflected process (given that the background state at time [math] is $i$). (ii) Then we use the material of [6, Section 3] on the all-time minimum for an sp-map; clearly, one obtains results for the all-time maximum of an sn-map by looking at the negative of the process. (iii) This means that we can make use of [6, Thm. 3.2] that provides the transform of the all-time minimum $\underline{Y}$ for an sp-map $(Y_{t})_{t\geqslant 0}$ (or, alternatively, we could use [10, Section 6]). Note that this transform also covers the epoch the minimum is attained, but (as we are interested in $\underline{Y}$ only) one has to set the corresponding parameter in the transform equal to [math]. (iv) We thus obtain that, with $\underline{J}$ the state of the background process at the epoch the maximum is attained, the matrix of which the $(i,j)$-th element is given by

[TABLE]

can be written as

[TABLE]

for a matrix $F$ with non-negative elements and row-sums equal to 1, and a defective transition rate matrix $E$ (and $I$ defining the $d$-dimensional identity matrix). (v) Observe that, by e.g. [3, Prop. III.4.1(iii)], this defines a phase-type distribution. $\Box$

Due to Lemma 5, by [3, Prop. III.4.1(i)],

[TABLE]

For the moment we assume that, for vectors ${\boldsymbol{\zeta}}$ and ${\boldsymbol{\eta}}$, the above expression allows the expansion

[TABLE]

a sufficient condition being that the eigenvalues of $E$ be distinct; in Remark 4 we comment on adaptations when (11) does not apply. In a way the above assumption can be considered as the sn-counterpart of the assumption (10) that we imposed in the sp-case. It follows that

[TABLE]

so that

[TABLE]

Recall from Lemma 4 that the transform of ${\mathbb{E}}_{i}{\rm e}^{-\vartheta\,\tau(y)}$ (with respect to $y$, that is) can be expressed in terms of the transform of $\underline{X}_{T}$. In addition, we remark that an expression for this transform of $\underline{X}_{T}$, which we denote by

[TABLE]

is given in e.g. [13, Thm. 1.(ii), Eqn. (17)] and [10, Prop. 6.1]. Now, similarly to what we have done in the sp-case, $\gamma(\vartheta)$ can be obtained appealing to (9) and Lemma 3. To this end, we have that

[TABLE]

with $\nu_{i}(\vartheta)$ as defined (and evaluated) in the sp-case. With (recalling Lemma 4)

[TABLE]

we have (by swapping the order of the integrals)

[TABLE]

where $\bar{\Omega}^{\prime}_{i}(\eta,\vartheta)$ is understood as the derivative of $\bar{\Omega}_{i}(\eta,\vartheta)$ with respect to $\eta$. Subtracting $({\mathbb{E}}Q_{0})^{2}$, we thus arrive at the following result.

Theorem 4**.**

If the expansion $(\ref{expand})$ applies, then, for any $\vartheta\geqslant 0$,

[TABLE]

Remark 4**.**

Regarding condition $(\ref{expand})$, we remark that one may have situations in which mixed exponential-polynomial terms show up in $p_{i}(x)$ (due to non-simple eigenvalues). These can be dealt with as in Remark 3.* $\diamond$*

6. Markov additive input: examples

In this section we treat a number of illustrative examples. In the first example we present explicit computations for one of the most frequently used maps, viz. Markov-modulated fluid. The second example is such that the structural properties, that have been derived for the Lévy case ($r(\cdot)$ being non-negative, non-increasing, convex), still apply, as well as the tail behavior agreeing with the tail of the busy-period distribution. The third example then shows that in general $r(\cdot)$ is not non-negative, non-increasing, convex; in addition the rate of decay of $r(t)$ for $t$ large is shown to potentially differ from that of the busy-period tail, depending on properties of the modulating process $(J_{t})_{t\in{\mathbb{R}}}$.

6.1. Markov-modulated fluid with a general number of states

We define the Markov-modulated fluid model as follows. In any of the states of the background process, the process behaves as a (state-specific) drift. Let $\mu_{i}$ be the constant drift when the background state is $i$, where $i\in\{1,\ldots,d\}$; assume that the $\mu_{i}$ s do not equal [math]. Let ${\mathscr{E}}_{-}$ the subset of $\{1,\ldots,d\}$ such that $\mu_{i}<0$ and let ${\mathscr{E}}_{+}$ the subset of $\{1,\ldots,d\}$ such that $\mu_{i}>0$. We start by considering the case of a general dimension $d\in\{2,3,\ldots\}.$ We first explain how ${\boldsymbol{\nu}}(\vartheta)$ and ${\boldsymbol{\xi}}(y,\vartheta)$ can be found.

• $\circ$

It is straightforward to determine the mean rate $\nu_{i}(\vartheta).$ It takes some elementary algebra to verify that ${\boldsymbol{\nu}}(\vartheta)=-(\mathscr{Q}-\vartheta I)^{-1}{\boldsymbol{\mu}}.$

• $\circ$

We set up a system of linear differential equations for the $\xi_{i}(y,\vartheta)$. Using standard arguments, as $\delta\downarrow 0$,

[TABLE]

It thus follows that, where $\xi^{\prime}(y,\vartheta)$ corresponds to differentiation with respect to $y$,

[TABLE]

In self evident notation, this is, in vector/matrix form,

[TABLE]

with $\Psi(\vartheta):={\mathscr{M}}^{-1}(\mathscr{Q}-\vartheta I)$ and ${\mathscr{M}}:={\rm diag}\{{\boldsymbol{\mu}}\}.$ By solving the system of linear differential equations, it thus follows that ${\boldsymbol{\xi}}(y,\vartheta)={\rm e}^{-\Psi(\vartheta)\,y}{\boldsymbol{\xi}}(0+,\vartheta),$ where

[TABLE]

Observe that ${\xi}_{i}(0+,\vartheta)=1$ for $i\in{\mathscr{E}}_{-}$, whereas ${\xi}_{i}(0+,\vartheta)<1$ for $i\in{\mathscr{E}}_{+}$; in the latter case ${\xi}_{i}(0+,\vartheta)$ is the transform of the queue’s busy period given the background process is in state $i$ at the start of this busy period.

6.2. Markov-modulated fluid with two states

In this case the steady-state distribution follows from the results of [11]. Write $q_{1}:=\widehat{q}_{1}=q_{12}$, $q_{2}:=\widehat{q}_{2}=q_{21}$, and $\bar{q}:=q_{1}+q_{2}$. The steady-state distribution of the background process is ${\boldsymbol{\pi}}=(q_{2}/\bar{q},q_{1}/\bar{q}).$ Following [11], we assume that $\mu_{1}>0$ and $\mu_{2}<0$. The average drift is

[TABLE]

which we throughout assume to be negative. We define $\lambda:=q_{2}/\mu_{1}+q_{1}/\mu_{2}>0.$ As argued in [11],

[TABLE]

This entails that

[TABLE]

Hence,

[TABLE]

and

[TABLE]

We now evaluate ${\boldsymbol{\nu}}(\vartheta)$:

[TABLE]

We now concentrate on determining $\xi_{1}(y,\vartheta)$ and $\xi_{2}(y,\vartheta).$ In this special case there is a direct approach; we start by analyzing $\xi_{2}(y,\vartheta).$ Let $N(x)$ be Poisson distributed with mean $x>0$. Let $\tau$ denote a busy period (starting at workload level 0, with the background process just having turned to state 1), and $\tau_{1},\tau_{2},\ldots$ i.i.d. copies of $\tau$; also, $\pi(\vartheta):={\mathbb{E}}\,{\rm e}^{-\vartheta\tau}$. Then, using basic self-similarity properties, given that $J_{0}=2$,

[TABLE]

In other words, with $\varrho(\vartheta):=-(\vartheta+q_{2}(1-\pi(\vartheta)))/\mu_{2}$,

[TABLE]

It is not hard to see that

[TABLE]

So we are left with determining $\pi(\vartheta)$. This we do by relating $\tau$ to the busy period $\tau^{\circ}$ in an M/M/1 queue with arrival rate $\alpha$, service rate $\beta$, and depletion rate $R$. With $\pi^{\circ}(\vartheta)\equiv\pi^{\circ}(\vartheta\,|\,\alpha,\beta,R):={\mathbb{E}}\,{\rm e}^{-\vartheta\tau^{\circ}}$, we have

[TABLE]

Now $\pi^{\circ}(\vartheta)$ can be found be solving the above second-order equation and picking the correct root. A straightforward geometric argument then yields that

[TABLE]

It then follows that

[TABLE]

Now $\gamma(\vartheta)$ can be explicitly evaluated. Define $c_{1}(\vartheta):=\pi(\vartheta)$, $c_{2}(\vartheta):=1$, and

[TABLE]

It thus follows that

[TABLE]

For ease we renormalize time and space such that $\mu_{1}=q_{2}=1$; we simply put $\mu:=-\mu_{2}>0$ and $q:=q_{1}.$ Hence,

[TABLE]

which equals

[TABLE]

The asymptotic behavior of $c(t):={\mathbb{C}}{\rm ov}(Q_{0},Q_{t})$ for $t$ large can be found by inspecting the singularities of $\gamma(\vartheta).$ More specifically, as described in full detail in [1, Section 3], the rightmost singularity of $\gamma(\vartheta)$ in the left half-plane determines the tail behavior of $c(t)$. In our case there is a pole at $\vartheta_{\rm p}=-\bar{q}=-q-1$. In addition, there is a branching point where $\pi(\vartheta)$ has a branching point, which is at

[TABLE]

Noting that $2\sqrt{q\mu}\geqslant-\mu-q$,

[TABLE]

We conclude that in this example the tail behavior of $c(t)$ is determined by $\vartheta_{\rm b}.$ More specifically,

[TABLE]

6.3. Example illustrating tail behavior of correlation

In this example the background process has a cyclic structure: $q_{12}=q_{23}=\cdots=q_{d-1,d}=q_{d1}=q$, $\mu_{1}>0$, and $\mu_{i}=-\infty$ for $i\in\{2,\ldots,d\}.$ To make the notation more compact, we again renormalize time and space by putting $\mu_{1}=q=1.$ Strictly speaking, this example is not a map (due to $\mu_{i}=-\infty$ for $i\in\{2,\ldots,d\}$), but it can be approximated arbitrarily closely by a map. By this example we show that the tail behavior of $c(\cdot)$ does not necessarily match that of the tail of the busy period distribution, depending on specific features of the matrix $\mathscr{Q}.$

As it turns out, in this setting we can compute ${\boldsymbol{p}}(x)$ and ${\boldsymbol{m}}(x)$ explicitly. It is easily seen that busy periods are exponentially distributed with mean $1$ and that during such busy periods the background state is 1. When the background process is in $i\in\{2,\ldots,d\}$, the workload level is [math]. Hence $p_{1}(x)={\rm e}^{-x}/d$, whereas ${\mathbb{P}}(Q_{0}=0,J_{0}=i)=1/d$ for $i\in\{2,\ldots,d\}$. Also,

[TABLE]

where $E_{d-1}$ denotes an Erlang random variable with $d-1$ phases, each phase being exponentially distributed with parameter 1. As

[TABLE]

we conclude that

[TABLE]

It thus follows that

[TABLE]

which simplifies to

[TABLE]

This expression has poles at $\vartheta=-1$ and at $\vartheta=-1+{\rm e}^{2\pi{\rm i}\,k/d}$, for $k\in\{1,2,\ldots,d-1\}$; observe that, in the expression between the brackets in (12), the zero in the denominator for $\vartheta=0$ is compensated by a zero in the corresponding numerator, entailing that $\vartheta=0$ is not a pole. We now show that increasing $d$ drastically affects the nature of the asymptotics of $c(t)$.

• $\circ$

We start with $d=2.$ Then

[TABLE]

This we rewrite to

[TABLE]

Explicit inversion yields

[TABLE]

In this case $\gamma(\vartheta)$ has poles at $\vartheta=-1$ and $\vartheta=-2.$ As mentioned, the right-most pole dominates the asymptotic behavior. We have $c(t)\,{\rm e}^{t}\to 1$ as $t\to\infty.$

• $\circ$

We continue with $d=3$. Now, besides the usual pole at $\vartheta=-1$, there are two more poles at $\vartheta={\rm e}^{2\pi{\rm i}/3}=(-3+{\rm i}\,\sqrt{3})/2$ and ${\rm e}^{4\pi{\rm i}/3}=(-3-{\rm i}\,\sqrt{3})/2$, both of them having real part $-3/2$. By comparing these poles, it follows that the pole with the largest real part remains $\vartheta=-1$. We have

[TABLE]

eventually leading to

[TABLE]

and thus $c(t)\,{\rm e}^{t}\to 2/3$ as $t\to\infty.$ Conclude that we see essentially the same behavior as for $d=2.$

• $\circ$

For $d=4$ we have

[TABLE]

In this case the poles are, besides $\vartheta=-1$, also $\vartheta=-2,$ $-1+{\rm i}$, and $-1-{\rm i}.$ After tedious but straightforward calculus,

[TABLE]

so that $\lim_{t\to\infty}c(t)\,{\rm e}^{-t}$ does not exist. More specifically,

[TABLE]

Conclude that in this case there are three poles with the same ‘magnitude’ (in terms of the size of their real parts), which consequently all appear in the asymptotic behavior.

• $\circ$

We now consider $d\in\{5,6,\ldots\}$. The crucial observation is that, for $d$ larger than $5$, the poles $-1+{\rm e}^{2\pi{\rm i}\,/d}$ and $-1+{\rm e}^{2\pi{\rm i}\,/d}$ are dominant, as these have a real part that is strictly larger than $-1$; this real part is $\phi_{d}:=-1+\cos(2\pi/d).$ As in the case of $d=4$, $\lim_{t\to\infty}c(t)\,{\rm e}^{-t}$ does not exist because of oscillating terms. Logarithmic asymptotics can be derived, though: the decay rate $\bar{c}$ equals

[TABLE]

So for $d$ up to 4 the decay rate corresponds with the pole at $-1$, whereas for larger values of $d$ the pole at $\cos(2\pi/d)$ takes over: we have found that $\bar{c}=-1+\max\{0,\cos(2\pi/d)\}.$

Remark 5**.**

The fact that the decay rate of $c(\cdot)$ does not necessarily coincide with that of the tail of the busy period distribution can be understood as follows. In [9] a coupling-based approach was developed for the case of Lévy input. The argumentation relied on writing $c(t)$ as ${\mathbb{E}}(Q_{0}(Q_{t}-Q^{\star}_{t}))$, with the process $(Q^{\star}_{t})_{t\geqslant 0}$ being constructed as follows: sample $Q_{0}^{\star}$ independently from $Q_{0}$, but use the same driving Levy process $(X_{t})_{t\geqslant 0}$ as for $(Q_{t})_{t\geqslant 0}$. It is now easily seen that a condition for $c(t)$ to be non-zero is that at least one of the two busy periods has not ended at time $t$; more precisely, we should have that $\underline{X}_{t}>-Q_{0}$ or $\underline{X}_{s}>-Q_{0}^{\star}$. This explains why in the Lévy case the tail behavior of $c(\cdot)$ is intimately connected with the tail distribution of the busy period. These ideas have been made precise in [9].

In the case of map input however, the states of the background processes corresponding to the sample $Q_{0}$ and $Q_{0}^{\star}$ do not necessarily coincide. As a consequence, we cannot use the same driving Lévy process $(X_{t})_{t\geqslant 0}$ after time 0, and hence the coupling idea does not apply. To remedy this complication, an idea could be to let the Lévy processes $(X_{t})_{t\geqslant 0}$ and $(X^{\star}_{t})_{t\geqslant 0}$ run independently till the background processes reach the same state, and apply the above coupling from that point on. It may, however, take relatively long before the background states meet. It is precisely this insight that explains why the structure of the background process $(J_{t})_{t\in{\mathbb{R}}}$ (represented by its transition rate matrix $\mathscr{Q}$) plays a role. This is nicely illustrated in the above example: when the dimension $d$ grows, the coupling takes longer, entailing that from a certain $d$ on (in the example $d=4$) the tail of $c(\cdot)$ starts to diverge from the tail of the busy-period distribution. * $\diamond$*

Appendix

Proof of Lemma 1: We believe this result has appeared, in various forms, in the literature. Several approaches can be followed to prove the claim; for completeness we include a compact proof here. The starting point is the relation

[TABLE]

for any non-negative $A$ and $B$ for which the above objects are well-defined.

First suppose that $f(\cdot)$ is non-negative; then, by (13),

[TABLE]

Then observe that

[TABLE]

using the fact that $f(\cdot)$ is non-decreasing in the (non-strict) inequality, and

[TABLE]

in the last equality. We thus establish the claim under the proviso that $f(\cdot)$ is non-decreasing and non-negative.

We now lift the condition that $f(\cdot)$ be non-negative. To this end, suppose that $f(0)<0$. Then, realizing that $a\mapsto f(a)-f(0)$ is non-negative, by applying the previous result,

[TABLE]

The result for $g(\cdot)$ follows analogously. $\Box$